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

void Align(){

	// Customize

	const Int_t numSens = 6;
	const Int_t maxEvents = 1000000;

	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);

	TH1F *DX = new TH1F("DX","DX",1000,-2,+2);
	TH1F *DY = new TH1F("DY","DY",1000,-2,+2);
	
	TCanvas *can0 = new TCanvas();
	TCanvas *can1 = new TCanvas();
	TCanvas *can2 = new TCanvas();
	TCanvas *can3 = new TCanvas();
	TCanvas *can4 = new TCanvas();
	TCanvas *can5 = 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;

	int nrHits1=0;
	int nrHits2=0;
	int nrHits3=0;
	int nrHits4=0;

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

	Double_t xa,xb,ya,yb,za_1,za_2,zb_1,zb_2;

	bool OneHit[6];

	for (Int_t j = 0 ; j < t->GetEntries(); j++) // loop on events
	{
		t->GetEvent(j);

		for (Int_t s = 0 ; s < 6 ; s++) OneHit[s] = false;

		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();
			
			if (!OneHit[0]) OneHit[0] = true;
			else
			  OneHit[0] = false;
		      }


		    If (name == SensName[5])
		      {
			x2 = point->GetX();
			y2 = point->GetY();
			z2 = point->GetZ();
		      
			if (!OneHit[5]) OneHit[5] = true;
			else
			  OneHit[5] = false;
		 
		      }

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

			if (!OneHit[1]) OneHit[1] = true;
			else
			  OneHit[1] = false;

		      }


		    if (name == SensName[2]) // first single sided - y
		      {
			ya = point->GetY();
			za_2 = point->GetZ();
		      
			if (!OneHit[2]) OneHit[2] = true;
			else
			  OneHit[2] = false;

		      }


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

			if (!OneHit[3]) OneHit[3] = true;
			else
			  OneHit[3] = false;

		      }


		    if (name == SensName[4]) // second single sided - y
		      {
			yb = point->GetY();
			zb_2 = point->GetZ();
		      
			if (!OneHit[4]) OneHit[4] = true;
			else
			  OneHit[4] = false;
		
		      }

		  }// end loop on hits

		
		if (OneHit[0] && OneHit[1] && OneHit[2] && OneHit[3] && OneHit[4] && OneHit[5])

		  {
		
		    Dx3->Fill(xa - CalcX(x1,z1,x2,z2,za_1) );
		    Dy3->Fill(ya - CalcY(y1,z1,y2,z2,za_2) );
		    Dx4->Fill(xb - CalcX(x1,z1,x2,z2,zb_1) );
		    Dy4->Fill(yb - CalcY(y1,z1,y2,z2,zb_2) );
		    
		    DX->Fill(x2-x1);
		    DY->Fill(y2-y1);

		  }
		
	} // End loop on events
	
	can0->cd();
	Dx3->Draw();
	//	Dx3->Fit("gaus");
	TSpectrum *sp0 = new TSpectrum(1,1);
	sp0->Search(Dx3,2,"",0.05);

	can1->cd();
	Dy3->Draw();
	Dy3->Fit("gaus");

	can2->cd();
	Dx4->Draw();
	Dx4->Fit("gaus");

	can3->cd();
	Dy4->Draw();
	Dy4->Fit("gaus");
	
	cout << "MAX: " << FindMax(Dy4) << endl;

	can4->cd();
	DX->Draw();

	can5->cd();
	DY->Draw();

	//rr->Draw("COLZ");
}


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;
  
}

Double_t FindMax(TH1F *h)
{

  return (h->GetBinCenter(h->GetBin(h->GetMaximum())));

}