/*
anaMomentum_new.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 = 1000;
  bool verbose = false;
  double pi = TMath::Pi();

  int p1 = 5;
  //  int theta1 = 7;

  TString outPath="tmpOutput/p";
  outPath += p1;
  outPath += "/";

  TString outfile=outPath+"Momentum_Resolution.root";                  //root file to store histogram for dtheta
    
  for(int theta1 =5; theta1<6; theta1++)
    {
      
      TString storePath="tmpOutput/p";
      storePath += p1;
      storePath += "/theta";
      storePath += theta1;
      storePath += "/";
      
      TString directoryName="theta_";
      directoryName += theta1;
      
      TFile *dtheta_histo = new TFile(outfile, "update");
      dtheta_histo->mkdir(directoryName);
      
      TString simInput=storePath+"points_";
      simInput += p1;
      simInput += "_";
      simInput += theta1;
      simInput += ".root";
      
      TString storehisto="dtheta_";                                 //name of histogram in file dtheta_histo.root it is like dtheta_i , where i = 1,2... 15
      storehisto += p1;
      storehisto += "_";
      storehisto += theta1;  
 
      TString MCmom="MC_Momentum_";
      MCmom += p1;
      MCmom += "_";
      MCmom += theta1;  
      
      TString Recomom="Reco_Momentum_";
      Recomom += p1;
      Recomom += "_";
      Recomom += theta1;  
      
      TString Resomom="Resolution_Momentum_";
      Resomom += p1;
      Resomom += "_";
      Resomom += theta1;  
      
      TString title = "Theta determined by difference in track after and before dipole for P = ";
      title += p1;
      title += " GeV/c";
      title += " for theta ";
      title += theta1; 
      title += " degrees"; 
      
      cout<<simInput<<endl;
      
      //to open data file for MC Points and tracks
      TFile *mcfile = new TFile(simInput);
      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);
      
      //To define Histograms
      TH1F *mcTheta = new TH1F("MCHit Theta", "Theta fitted by MCHits", 500, 0.00, 0.1);
      TH1F *deltaMcTheta = new TH1F(storehisto,title , 300, 0.001, 0.7);
      
      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", 800, -0.005, 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(MCmom, "Momentum", 1500, 0.0, 15.0);
      TH1F *mom_reco = new TH1F(Recomom, "Reco Momentum", 1500, 0.0, 15.0);
      TH1F *mom_reso = new TH1F(Resomom, "Momentum Resoluition", 1500, -1.0, 1.0);
      
      
      //Define some global variables
      TVector3 vecmc, vecrc;
      double mcX=0, mcY=0, mcZ=0, rcX=0, rcY=0, rcZ=0, mcPx=0, mcPy=0, mcPz=0 ;
      double trackthe=0, mctheta=0, deltamctheta=0, mctheta1=0, mcfirstthe=0, dtheta=0, momentum =0, momentum_reco=0, momentum_reso=0 ;
      double mc_x[16]=0, mc_y[16]=0, mc_z[16]=0, mc_P[48]=0, mc1_x[48]=0, mc1_y[48]=0, mc1_z[48]=0, r_mc[16]=0, r1_mc[48]=0, rc_x[16]=0, rc_y[16]=0, rc_z[16]=0, r_rc[16]=0;
      double mc2_x[16]=0, mc2_y[16]=0, mc2_z[16]=0, r2_mc[16]=0;
      double r_mc_err[16]=0, z_mc_err[16]=0, r_rc_err[16]=0, z_rc_err[16]=0;
      
      double par0 = 0.09029, par1 = -0.6945, par2 = 0.2994; //fts5richfts6
      //      double par0 = 0.1274, par1 = -0.8118, par2 = 0.2956; //fts5fts6rich
      
      //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);
	  
	  cout<<"the Event No. is "<<i<<endl;
	  if(mc_array->GetEntriesFast()==48)
	    {
	      for (int j=0; j<mc_array->GetEntriesFast(); 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);
		      
		      Int_t iTrack = point->GetTrackID();
		      PndMCTrack *track = (PndMCTrack*)mc_track->At(iTrack);
		      
		      if(track->GetMotherID()==-1){
			trackthe = track->GetMomentum().Theta();
			momentum = track->GetMomentum().Mag();
			mcPx = track->GetMomentum().X();
			mcPy = track->GetMomentum().Y();
			mcPz = track->GetMomentum().Z();
			cout<<"the momentum_mag, mcPx, mcPy, mcPz are "<<momentum<< " " <<mcPx<<" "<<mcPy<<" "<<mcPz<<endl;
		      }
		      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]=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;
		      
		      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);
		      
		      Int_t iTrack = point->GetTrackID();
		      PndMCTrack *track = (PndMCTrack*)mc_track->At(iTrack);
		      
		      if(track->GetMotherID()==-1){
			trackthe = track->GetMomentum().Theta();
			momentum = track->GetMomentum().Mag();
			mcPx = track->GetMomentum().X();
			mcPy = track->GetMomentum().Y();
			mcPz = track->GetMomentum().Z();
			cout<<"the momentum_mag, mcPx, mcPy, mcPz are "<<momentum<< " " <<mcPx<<" "<<mcPy<<" "<<mcPz<<endl;
		      }
		      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]=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;
		      
		    }
		  }
		}//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;
	      
	      momentum_reco = par0 + par1*sin(deltamctheta/2) + par2/sin(deltamctheta/2);
	      
	      cout<<"mctheta  is "<<mctheta<<endl;
	      cout<<"mctheta1 is "<<mctheta1<<endl;
	      cout<<"difference in mctheta and mctheta1 is "<<deltamctheta<<endl;
	      
	      
	      momentum_reso = (momentum_reco-momentum)/momentum;
	      
	      dtheta = mcfirstthe-trackthe;
	      mom->Fill(momentum);
	      mom_reco->Fill(momentum_reco);
	      mom_reso->Fill(momentum_reso);
	      deltaMcTheta->Fill(deltamctheta);
	      mcTheta->Fill(mctheta);
	      refTheta->Fill(trackthe);
	      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();
      deltaMcTheta->DrawCopy();
      deltaMcTheta->GetXaxis()->SetTitle("#Delta#theta_{x} (in rad)");
      deltaMcTheta->GetYaxis()->SetTitle("counts");
      deltaMcTheta->GetXaxis()->SetTitleOffset(1.3);
      deltaMcTheta->SetTitleOffset(1.3, "Y");
      dtheta_histo->cd(directoryName);
      deltaMcTheta->Write();
      
      can1->cd(2); //dTheta ->DrawCopy();   dTheta->Fit("total","R");
      mom->DrawCopy();
      dtheta_histo->cd(directoryName);
      mom->Write();
      
      can1->cd(3); //refTheta->DrawCopy();
      mom_reco->DrawCopy();
      dtheta_histo->cd(directoryName);
      mom_reco->Write();
      
      can1->cd(4); //mom->DrawCopy();
      mom_reso->DrawCopy();
      dtheta_histo->cd(directoryName);
      mom_reso->Write();
      
      //      gStyle->SetOptFit();
      dtheta_histo->Close(); 
    }
  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;
  
}