/*
 * TrackHits.C
 *
 *  Created on: Mar  5, 2010
 *      Author: Simone Bianco
 */

void Align(){

	// Customize

	const Int_t numSens = 6;
   const Int_t maxEvents = 10000;

	TString SensName[numSens];
	
	SensName[0] = "/TS_1/TTVol_0/TTDouble_0/StripActiveTD1_0";
	SensName[1] = "/TS_1/TTVol_0/TTSingle_0/StripActiveTS3a_0";
	SensName[2] = "/TS_1/TTVol_0/TTSingle_0/StripActiveTS3b_0";
	SensName[3] = "/TS_1/TTVol_0/TTSingle_0/StripActiveTS4a_0";
	SensName[4] = "/TS_1/TTVol_0/TTSingle_0/StripActiveTS4b_0";
	SensName[5] = "/TS_1/TTVol_0/TTDouble_0/StripActiveTD2_0";


	// load libs

	gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");

	//

	TString directory = gSystem->Getenv("VMCWORKDIR");
	TString geomFile = directory + "/geometry/TrackingStation.root";
	TString HitsFile = "hits.root";


	// Loading the geometry and defining the geo handler
	
	TFile *geo = new TFile(geomFile);

	TGeoManager *myGeo = geo->Get("FAIRGeom");

	fGeoH = new PndMvdGeoHandling(myGeo);


	TH1F *Dx3 = new TH1F("Dx3","Dx3",1000,-0.96,+0.96);
	TH1F *Dy3 = new TH1F("Dy3","Dy3",1000,-0.96,+0.96);
	TH1F *Dx4 = new TH1F("Dx4","Dx4",1000,-0.96,+0.96);
	TH1F *Dy4 = new TH1F("Dy4","Dy4",1000,-0.96,+0.96);

	
	TCanvas *can0 = new TCanvas();
	TCanvas *can1 = new TCanvas();
	TCanvas *can2 = new TCanvas();
	TCanvas *can3 = new TCanvas();

	TString name = "";
	
	// Load the hits
	TFile *f = new TFile(HitsFile);

	TTree *t=(TTree *) f->Get("cbmsim");
	
	TClonesArray* tr_array=new TClonesArray("PndMvdHit");
	t->SetBranchAddress("MVDHitsStrip",&tr_array);//Branch names

    	cout << "Events: " << t->GetEntries() << endl;


	Double_t x1,y1,z1,x2,y2,z2;

	Double_t xa,xb,ya,yb,za_1,za_2,zb_1,zb_2;
 
    
	for (Int_t j = 0 ; j < t->GetEntries() && j<maxEvents; j++) // loop on events
	{
		t->GetEvent(j);

		Int_t nrHits0=0;
		Int_t nrHits1=0;
		Int_t nrHits2=0;
		Int_t nrHits3=0;
		Int_t nrHits4=0;
		Int_t nrHits5=0;
    

		if (j%10000 == 0) cout << "Ev. " << j << endl;

		for (Int_t y = 0 ; y < tr_array->GetEntries() ; y++) // loop on hits
		  {

		    PndMvdHit* point = (PndMvdHit*)tr_array->At(y);
		    //cout<< "index: "<<point->GetClusterIndex() << endl;
          //cout<<" (x,y)=("<<point->GetX()<<","<<point->GetY()<<") charge="<<point->GetEloss()<<endl;
		
		    name =  fGeoH->GetPath(point->GetDetName());

		    if (j < 10) cout << name << endl;
		    

		      if (name == SensName[0])
		      {
			      x1 = point->GetX();
			      y1 = -point->GetY();
			      z1 = point->GetZ();
			      nrHits0++;
		      }
			

		      if (name == SensName[5])
		      {
			      x2 = point->GetX();
			      y2 = -point->GetY();
			      z2 = point->GetZ();
               nrHits1++;
            }

		    
		      if (name == SensName[1]) // first single sided - x
		      {
			      xa = point->GetX();
			      za_1 = point->GetZ();
               nrHits2++;
	         }


		      if (name == SensName[2]) // first single sided - y
		      {
			      ya = point->GetY();
			      za_2 = point->GetZ();
		         nrHits3++;
		      }


		      if (name == SensName[3]) // second single sided - x
		      {
			      xb = point->GetX();
			      zb_1 = point->GetZ();
               nrHits4++;
		      }


		      if (name == SensName[4]) // second single sided - y
		      {
			      yb = point->GetY();
			      zb_2 = point->GetZ();
			      nrHits5++;
		      }

		   }// end loop on hits

		   cout << "***hit:"<<endl;
         cout << " n1="<<nrHits0<<endl;
         cout << " n2="<<nrHits1<<endl;
         cout << " n3="<<nrHits2<<endl;
         cout << " n4="<<nrHits3<<endl;
         cout << " n5="<<nrHits4<<endl;
         cout << " n6="<<nrHits5<<endl;

         if ((nrHits0==1) && (nrHits1==1))
	      {
            if (nrHits2==1)
	      Dx3->Fill(xa - CalcX(x1,z1,x2,z2,za_1) );
            if (nrHits3==1)
	      Dy3->Fill(ya - CalcY(y1,z1,y2,z2,za_2) );
            if (nrHits4==1)
	      Dx4->Fill(xb - CalcX(x1,z1,x2,z2,zb_1) );
            if (nrHits5==1)
	      Dy4->Fill(yb - CalcY(y1,z1,y2,z2,zb_2) );
	      }
		
	} // End loop on events
	
	can0->cd();
	Dx3->Draw();

	can1->cd();
	Dy3->Draw();

	can2->cd();
	Dx4->Draw();

	can3->cd();
	Dy4->Draw();

	cout << fff << endl;
}


Double_t CalcX(Double_t xFirst, Double_t zFirst, Double_t xLast, Double_t zLast, Double_t z)
{
  
  Double_t ratio = z/(zLast - zFirst);

  Double_t pos = xFirst + (xLast-xFirst)*ratio;

  return pos;
  
}

Double_t CalcY(Double_t yFirst, Double_t zFirst, Double_t yLast, Double_t zLast, Double_t z)
{
  
  Double_t ratio = z/(zLast - zFirst);

  Double_t pos = yFirst + (yLast-yFirst)*ratio;

  return pos;
  
}