/*
anaMomentum.C
the macro measures difference between monte carlo momentum and momentum from hit point fitting

Auther: Manoj
*/

{
#include "TCanvas.h"
  std::cout;
  
  //load library by run macro
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
  rootlogon();
  
  //Timer
  TStopwatch timer;
  timer.Start();
  
  //variable defination 
  int nEvents = 10000;
  bool verbose = false;
  double pi = TMath::Pi();

  //to open data file for MC Points and tracks
  TFile *mcfile = new TFile("points.root");
  TTree *tree = (TTree*)mcfile->Get("cbmsim");
  TClonesArray *mc_array = new TClonesArray("PndFtsPoint");
  tree->SetBranchAddress("FTSPoint", &mc_array);

  TClonesArray *mc_track = new TClonesArray("PndMCTrack");
  //  TClonesArray *mc_track = new TClonesArray("PndFtsTrackerIdeal");
  tree->SetBranchAddress("MCTrack", &mc_track);

/*
  TFile *recofile = new TFile("../run/reco_complete.root");
  TTree *tree1 = (TTree*)mcfile->Get("cbmsim");
  TClonesArray *hit_array = new TClonesArray("PndTrack");
  tree->SetBranchAddress("FtsIdealGenTrack", &hit_array);
*/
  //To define Histograms
  TH1F *mcTheta = new TH1F("MCHit Theta", "Theta fitted by MCHits", 300, 0.001, 0.1);
  TH1F *deltaMcTheta = new TH1F("delta MC theta","Theta determined by difference in track after and before dipole",300, 0.001, 0.2);

  TH1F *refTheta = new TH1F("MCTrue Theta", "The track theta set in Generator", 300, 0.001, 0.1);
  TH1F *dTheta = new TH1F("FirstMCHit-MCTrue", "delta theta between FirstMCHit and MCTrue", 300, -0.003, 0.003);

  TH1F *dTheta0 = new TH1F("RecoHit-MCTrue","delta theta between RecoHit and MCTrue theta",300, -0.003,0.003);//theta resolution
  //	TH1F *dTheta1 = new TH1F("FirstRecoHit-MCTrue","",300, -0.003,0.003);
  TH1F *dTheta2 = new TH1F("MCHit-MCTrue","",300, -0.003,0.003);
  TH1F *dTheta3 = new TH1F("RecoHit-MCHit","delta theta between RecoHit and MCHit",300, -0.003,0.003); //digi&reco effect


  TF1 *trackFit = new TF1("trackFit", "pol1",-40, 40);

  TH1F *mom = new TH1F("MC Momentum", "Momentum", 400, 4.0, 6.0);


  //Define some global variables
  TVector3 vecmc, vecrc;
  double mcX=0, mcY=0, mcZ=0, rcX=0, rcY=0, rcZ=0;
  double trackthe, mctheta, deltamctheta, mctheta1, mcfirstthe, dtheta, momentum;
  double mc_x[16], mc_y[16], mc_z[16], mc1_x[48], mc1_y[48], mc1_z[48], r_mc[16], r1_mc[48], rc_x[16], rc_y[16], rc_z[16], r_rc[16];
  double mc2_x[16], mc2_y[16], mc2_z[16], r2_mc[16];
  double r_mc_err[16], z_mc_err[16], r_rc_err[16], z_rc_err[16];

  //Loop the event and get Theta
//  cout<<"No. Of Events"<<tree->GetEntriesFast()<<endl;
  for(int i =0; i<nEvents&&i<tree->GetEntriesFast(); i++)
  {
	  tree->GetEntry(i);
//	  tree1->GetEntry(i);

	  cout<<"the Event No. is "<<i<<endl;
	  if(mc_array->GetEntriesFast()==48)
	  {
		  for (int j=0; j<mc_array->GetEntriesFast(); j++)
		  {
//		  if(mc_array.fChamberID==1 || mc_array.fChamberID==2)
//		  {
			  //			  if(hit_array->GetEntriesFast()!=mc_array->GetEntriesFast()) continue;
			  //			  FtsIdealGenTrack *hit = (FtsIdealGenTrack*)hit_array->At(j);
			  PndFtsPoint *point = (PndFtsPoint*)mc_array->At(j);
			  if(0<=j && j<16){
			  if(point->GetChamberID()==1 || point->GetChamberID()==2){

				  cout<<"chamber ID for stations in FT"<<point->GetChamberID()<<endl;

				  mcX = point.GetX();
				  mcY = point.GetY();
				  mcZ = point.GetZ();

				  cout<<"the mcX, mcY, mcZ are "<<mcX<<" "<<mcY<<" "<<mcZ<<endl;
				  vecmc.SetXYZ(mcX, mcY, mcZ);
/*
			  	  rcX = hit->GetX();
			  	  rcY = hit->GetY();
			  	  rcZ = hit->GetZ();

			  	  cout<<"the rcX, rcY, rcZ are "<<rcX<<","<<rcY<<","<<rcz<<endl;
			  	  vecrc.SetXYZ(rcX, rcY, rcZ);
*/
				  Int_t iTrack = point->GetTrackID();
				  PndMCTrack *track = (PndMCTrack*)mc_track->At(iTrack);

				  if(track->GetMotherID()==-1){
					  trackthe = track->GetMomentum().Theta();
					  momentum = track->GetMomentum().Mag();
				  }
				  cout<<"the track theta is "<<trackthe<<endl;

				  mc_x[j]= vecmc.X();
				  mc_y[j]= vecmc.Y();
				  mc_z[j]= vecmc.Z();

				  r_mc[j]=TMath::Sqrt(mc_x[j]*mc_x[j]);//+mc_y[j]*mc_y[j]);

				  std::cout<<"the r_mc["<<j<<"] is "<<r_mc[j]<<std::endl;
				  std::cout<<"the mc_z["<<j<<"] is "<<mc_z[j]<<std::endl;

				  r_mc_err[j]=0;
				  z_mc_err[j]=0;
/*
			  	  rc_x[j]=vecrc.X();
			  	  rc_y[j]=vecrc.Y();
			  	  rc_z[j]=vecrc.Z();
			  	  r_rc[j]=TMath::Sqrt(rc_x[j]*rc_x[j]+rc_y[j]*rc_y[j]);
			  	  std::cout<<"the r_rc[j] is "<<r_rc[j]<<std::endl;
			  	  std::cout<<"the rc_z[j] is "<<rc_z[j]<<std::endl;
			  	  r_rc_err[j]=0;
			  	  z_rc_err[j]=0;
*/
				  double tempmc = TMath::Sqrt(mc_x[0]*mc_x[0]);//+mc_y[0]*mc_y[0]);
				  mcfirstthe = TMath::ATan(tempmc/(mc_z[0]));
				  cout<<"the first hit theta is "<<mcfirstthe<<endl;
			  	  }
			  }
			  if(32<=j && j<48){
			  if(point->GetChamberID()==5 || point->GetChamberID()==6){

				  cout<<"chamber ID for stations in FT"<<point->GetChamberID()<<endl;

				  mcX = point.GetX();
				  mcY = point.GetY();
				  mcZ = point.GetZ();

				  cout<<"the mcX, mcY, mcZ are "<<mcX<<" "<<mcY<<" "<<mcZ<<endl;
				  vecmc.SetXYZ(mcX, mcY, mcZ);
/*
			  	  rcX = hit->GetX();
			  	  rcY = hit->GetY();
			  	  rcZ = hit->GetZ();

			  	  cout<<"the rcX, rcY, rcZ are "<<rcX<<","<<rcY<<","<<rcz<<endl;
			  	  vecrc.SetXYZ(rcX, rcY, rcZ);
*/
				  Int_t iTrack = point->GetTrackID();
				  PndMCTrack *track = (PndMCTrack*)mc_track->At(iTrack);
/*
				  if(track->GetMotherID()==-1){
					  trackthe = track->GetMomentum().Theta();
					  momentum = track->GetMomentum().Mag();
				  }
				  cout<<"the track theta is "<<trackthe<<endl;
*/
				  mc1_x[j]= vecmc.X();
				  mc1_y[j]= vecmc.Y();
				  mc1_z[j]= vecmc.Z();

				  r1_mc[j]=TMath::Sqrt(mc1_x[j]*mc1_x[j]);//+mc1_y[j]*mc1_y[j]);

				  std::cout<<"the r_mc["<<j<<"] is "<<r1_mc[j]<<std::endl;
				  std::cout<<"the mc_z["<<j<<"] is "<<mc1_z[j]<<std::endl;

//				  r_mc_err[j]=0;
//				  z_mc_err[j]=0;
/*
			  	  rc_x[j]=vecrc.X();
			  	  rc_y[j]=vecrc.Y();
			  	  rc_z[j]=vecrc.Z();
			  	  r_rc[j]=TMath::Sqrt(rc_x[j]*rc_x[j]+rc_y[j]*rc_y[j]);
			  	  std::cout<<"the r_rc[j] is "<<r_rc[j]<<std::endl;
			  	  std::cout<<"the rc_z[j] is "<<rc_z[j]<<std::endl;
			  	  r_rc_err[j]=0;
			  	  z_rc_err[j]=0;

				  double tempmc = TMath::Sqrt(mc_x[0]*mc_x[0]+mc_y[0]*mc_y[0]);
				  mcfirstthe = TMath::ATan(tempmc/(mc_z[0]));
				  cout<<"the first hit theta is "<<mcfirstthe<<endl;
*/}
			  	  }
		  	  }//loop all events containing for 48 hits

		  TGraphErrors *mctrack = new TGraphErrors(16, mc_z, r_mc, z_mc_err, r_mc_err);
		  mctrack->Fit("trackFit", "ONF");
		  mctheta = TMath::ATan(trackFit->GetParameter(1));
		  cout<<"the mc fit theta is "<<mctheta<<endl;

		  for(int k=32; k<48; k++){
			  int j=k-32;
			  mc2_x[j]= mc1_x[k];
			  mc2_y[j]= mc1_y[k];
			  mc2_z[j]= mc1_z[k];

			  r2_mc[j]= r1_mc[k];
		  }

		  TGraphErrors *mctrack1 = new TGraphErrors(16, mc2_z, r2_mc, z_mc_err, r_mc_err);
		  mctrack1->Fit("trackFit", "ONF");
		  mctheta1 = TMath::ATan(trackFit->GetParameter(1));
		  cout<<"the mc fit theta1 is "<<mctheta1<<endl;

		  deltamctheta = mctheta1-mctheta;

		  cout<<"mctheta  is "<<mctheta<<endl;
		  cout<<"mctheta1 is "<<mctheta1<<endl;
		  cout<<"difference in mctheta and mctheta1 is "<<deltamctheta<<endl;

/*
		  TGraphErrors *rctrack = new TGraphErrors(4,rc_z,r_rc,z_rc_err,r_rc_err);
		  rctrack->Fit("trackFit","ONF");
		  rctheta=TMath::ATan(trackFit->GetParameter(1));
		  cout<<"the rc fit theta is "<<rctheta<<endl;

	      dthe0 = rctheta-trackthe;
		  dthe2 = mctheta-trackthe;
	      dthe3 = rctheta-mctheta;
*/
		  dtheta = mcfirstthe-trackthe;
		  mom->Fill(momentum);
		  deltaMcTheta->Fill(deltamctheta);
		  mcTheta->Fill(mctheta);
		  refTheta->Fill(trackthe);
/*
		  dTheta0->Fill(dthe0);
		  dTheta2->Fill(dthe2);
		  dTheta3->Fill(dthe3);
		  */
		  dTheta->Fill(dtheta);
	  }
  }
  //draw the plots

  Double_t par[3];
  TF1 *total = new TF1("total", "gaus", -0.06, 0.06);
  total->SetLineColor(kRed);
  total->SetLineWidth(1);

  TCanvas *can1 = new TCanvas("Theta","MC, Track theata", 0, 0, 800, 800);
  can1->Divide(2,2);
  TPad* mypad = 0;

  can1->cd(1); mcTheta->DrawCopy();
  can1->cd(2); deltaMcTheta->DrawCopy();
  can1->cd(3); refTheta->DrawCopy();
  can1->cd(4); mom->DrawCopy();

  gStyle->SetOptFit();
  TCanvas *can2 = new TCanvas("Dealta Theta", "Track", 0, 0, 800, 800);
//  can2->Divide(2,2);

//  can2->cd(1);    dTheta0->DrawCopy();   dTheta0->Fit("total","R");
//    can2->cd(3);    dTheta1->DrawCopy(); //  dTheta1->Fit("total","R");
//  can2->cd(2);    dTheta2->DrawCopy();   dTheta2->Fit("total","R");
//  can2->cd(3);    dTheta3->DrawCopy();   dTheta3->Fit("total","R");
//  can2->cd(4);
  dTheta ->DrawCopy();   dTheta->Fit("total","R");


  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout<<endl;
  cout<<"the macro finished successfully"<<endl;
  cout<<"the real time is "<<rtime<<" seconds and the CPU time is "<<ctime<<" seconds."<<endl;

}