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

  //to open data file for MC Points and tracks

  TString inRecoFile  = "reco_complete.root";
  TString inDigiFile  = "digi_complete.root";
  TString inSimFile = "sim_complete.root";

  TFile *inFile = TFile::Open(inSimFile,"READ");
  TFile *fiDigiFile = TFile::Open(inDigiFile,"READ");
  TFile *fiReco1File = TFile::Open(inRecoFile,"READ");

  TTree *tree=(TTree *) inFile->Get("cbmsim") ;
  tree->AddFriend("cbmsim",inRecoFile);
  tree->AddFriend("cbmsim",inDigiFile);

  TClonesArray *ftsmcarray = new TClonesArray("PndFtsPoint");
  tree->SetBranchAddress("FTSPoint", &ftsmcarray);

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

  TClonesArray *trk_array = new TClonesArray("PndTrack");
  tree->SetBranchAddress("FtsIdealGenTrack", &trk_array);

  TClonesArray *ftsarray = new TClonesArray("PndFtsHit");
  tree->SetBranchAddress("FTSHit", &ftsarray);


  TList *bl = (TList*) fiDigiFile->Get("BranchList");
  FairRootManager *fairman = new FairRootManager();
  fairman->SetBranchNameList(bl);
  int ftsdetid = fairman->GetBranchId("FTSHit");
  cout<< "detector id for fts is "<<ftsdetid<<endl;


  //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 *rcTheta = new TH1F("reco theta","Theta determined by RecoHits",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, rctheta;
  double mc_x[48], mc_y[48], mc_z[48], mc1_x[48], mc1_y[48], mc1_z[48], r_mc[48], r1_mc[48], rc_x[48], rc_y[48], rc_z[48], r_rc[48];
  double mc2_x[16], mc2_y[16], mc2_z[16], r2_mc[16];
  double r_mc_err[48], z_mc_err[48], r_rc_err[48], z_rc_err[48];

  //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;
	  cout<<"the hit array no is "<<trk_array->GetEntriesFast()<<endl;
	  cout<<"the mc array no is "<<ftsmcarray->GetEntriesFast()<<endl;
	  //	  if(ftsmcarray->GetEntriesFast() > 0)
	  //	  {
		  for (int j=0; j<trk_array->GetEntriesFast(); j++)
		  {
//		  if(mc_array.fChamberID==1 || mc_array.fChamberID==2)
//		  {
//			  if(hit_array->GetEntriesFast()!=mc_array->GetEntriesFast()) continue;
			  PndTrack *trk = (PndTrack*)trk_array->At(j);
			  if(!trk)
			  {
				  cout<<"ERROR track "<< itrk << " does not exit"<<endl;
				  continue;
			  }
			  PndTrackCand *cand = trk->GetTrackCandPtr();
			  if(!cand)
			  {
				  cout<<"ERROR track "<< itrk << " has no candidate association"<<endl;
				  continue;
			  }

			  for(Int_t ihit = 0; ihit < cand->GetNHits(); ihit++){

			  PndTrackCandHit candhit = cand->GetSortedHit(ihit);
			  Int_t hitId = candhit.GetHitId();
			  Int_t detId = candhit.GetDetId();

			  cout<<"HitId: "<< hitId << " DetId: "<<detId<<endl;

			  Float_t x = 0, y = 0, z = 0, mcx = 0, mcy = 0, mcz = 0;
			  FairHit *hit = NULL;
			  FairMCPoint *point = NULL;

			 if(detId == FairRootManager::Instance()->GetBranchId("FTSHit"))
			 {
				 hit = (FairHit*) ftsarray->At(hitId);
				 if(!hit)
				 {
					 cout<<"Error fts " << hitId << " does not exist"<<endl;
				 }
				 else
				 {
					 if(hit->GetRefIndex()==-1) cout <<"Error fts: ref index -1"<<endl;
					 else
					 {
						 point = (PndFtsPoint*)(ftsmcarray->At(hit->GetRefIndex()));
						 mcx = point->GetX(); mcy = point->GetY(); mcz = point->GetZ();
						 
						 cout << "hits for reco X:"<< hit->GetX() <<" Y: " << hit->GetY() << " Z: " << hit->GetZ() << endl; 
					 }
				 }

				 cout<<"the mcX, mcY, mcZ are "<<mcx<<" "<<mcy<<" "<<mcz<<endl;
				vecmc.SetXYZ(mcx, mcy, mcz);

				  Int_t iTrack = point->GetTrackID();
				  PndMCTrack *track = (PndMCTrack*)mc_array->At(iTrack);

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

				  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(48, 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;


		  TGraphErrors *rctrack = new TGraphErrors(48,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;
		  dtheta = mcfirstthe-trackthe;
		  */
		  mom->Fill(momentum);

		  //		  mcTheta->Fill(mctheta);
		  //		  refTheta->Fill(trackthe);

		  //		  dTheta->Fill(dtheta);
		  //	  }
		  cout<<"===================================================================================="<<endl;
  }
  //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); rcTheta->DrawCopy();
  can1->cd(3); refTheta->DrawCopy();
  can1->cd(4); mom->DrawCopy();

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

  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;

}