Merge pull request #27 from eic/tracking_benchmark

Tracking benchmark
eic · Jun 27, 2024 · 0a85ff9 · 0a85ff9
2 parents c93e12d + 0d0c9c6
commit 0a85ff9
Show file tree

Hide file tree

Showing 4 changed files with 381 additions and 0 deletions.
diff --git a/benchmarks/tracking_performances/Script_widebin.sh b/benchmarks/tracking_performances/Script_widebin.sh
@@ -0,0 +1,64 @@
+#!/bin/bash
+# Script; shyam kumar; INFN Bari, Italy
+# [email protected]; [email protected]
+rm -rf truthseed/ realseed/ *.root
+mkdir -p truthseed/pi-/mom realseed/pi-/mom truthseed/pi-/dca realseed/pi-/dca
+mom_array=(0.5 1.0 2.0 5.0 10.0 15.0)
+particle_array=("pi-")
+filename=("tracking_output") 
+etabin_array=(-3.5 -2.5 -1.0 1.0 2.5 3.5)
+nevents=10000
+
+# run the simulation
+source ../epic/install/bin/thisepic.sh 
+dd_web_display --export -o epic_craterlake_tracking_only.root ../epic/install/share/epic/epic_craterlake_tracking_only.xml
+for ((i=0; i<${#mom_array[@]}; i++)); do
+npsim --compactFile ../epic/install/share/epic/epic_craterlake_tracking_only.xml --outputFile sim${mom_array[i]}.edm4hep.root --numberOfEvents $nevents --enableGun --gun.thetaMin 3*deg --gun.thetaMax 177*deg --gun.distribution eta --gun.particle pi- --gun.momentumMin ${mom_array[i]}*GeV --gun.momentumMax ${mom_array[i]}*GeV --gun.multiplicity 1 --random.seed 100000
+done
+# run the reconstruction
+source ../epic/install/setup.sh 
+source ../EICrecon/install/bin/eicrecon-this.sh
+
+for ((i=0; i<${#mom_array[@]}; i++)); do
+eicrecon \
+ -Pnthreads=1 \
+ -Pjana:debug_plugin_loading=1 \
+ -Pjana:nevents=$nevents \
+ -Pacts:MaterialMap=calibrations/materials-map.cbor \
+ -Ppodio:output_file="${filename}"_${mom_array[i]}.edm4eic.root \
+ -Pdd4hep:xml_files=../epic/install/share/epic/epic_craterlake_tracking_only.xml   \
+ -Ppodio:output_collections="MCParticles,CentralCKFTrajectories,CentralCKFTrackParameters,CentralCKFSeededTrackParameters,CentralTrackVertices" \
+ sim${mom_array[i]}.edm4hep.root
+done 
+
+# run the analysis
+for ((iparticle=0; iparticle<${#particle_array[@]}; iparticle++)); do
+# truth seeding
+for ((i=0; i<${#mom_array[@]}; i++)); do
+root -b -l -q Tracking_Performances.C'("'${filename}'","'${particle_array[iparticle]}'",'${mom_array[i]}',0.15,"")'
+done
+
+# real seeding
+for ((i=0; i<${#mom_array[@]}; i++)); do
+root -b -l -q Tracking_Performances.C'("'${filename}'","'${particle_array[iparticle]}'",'${mom_array[i]}',0.15,"Seeded")'
+done
+done
+cd truthseed/pi-/dca
+hadd final_hist_dca_truthseed.root *.root
+cd ../../../realseed/pi-/dca/
+hadd final_hist_dca_realseed.root *.root
+cd ../../../
+
+rm -rf Final_Results/ Debug_Plots/ 
+mkdir -p Final_Results/pi-/mom Final_Results/pi-/dca  Debug_Plots/truth/pi-/mom  Debug_Plots/truth/pi-/dca  Debug_Plots/real/pi-/mom Debug_Plots/real/pi-/dca
+# loop over particles
+for ((iparticle=0; iparticle<${#particle_array[@]}; iparticle++)); do
+for ((i=0; i<${#etabin_array[@]}-1; i++)); do
+xmax_hist=0.3 
+if [ $i == 2 || $i == 1 ]; then 
+xmax_hist=0.01 
+fi
+root -b -l -q doCompare_truth_real_widebins_mom.C'("'${particle_array[iparticle]}'",'${etabin_array[i]}','${etabin_array[i+1]}','$xmax_hist')'
+done
+done
+
diff --git a/benchmarks/tracking_performances/Tracking_Performances.C b/benchmarks/tracking_performances/Tracking_Performances.C
@@ -0,0 +1,118 @@
+// Code to extract the Tracking Performances
+// Shyam Kumar; INFN Bari, Italy
+// [email protected]; [email protected]
+
+#include "TGraphErrors.h"
+#include "TF1.h"
+#include "TRandom.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TMath.h"
+#define mpi 0.139  // 1.864 GeV/c^2
+
+void Tracking_Performances(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.15, TString name = "")
+{
+
+  // style of the plot
+   gStyle->SetPalette(1);
+   gStyle->SetOptTitle(1);
+   gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y");
+   gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y");
+   gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y");
+   gStyle->SetHistLineWidth(2);
+   gStyle->SetOptFit(1);
+   gStyle->SetOptStat(1);
+
+   TString dir = "";
+   TString dist_dir_mom = ""; TString dist_dir_dca = "";
+   if (name=="") {dist_dir_mom = "mom_resol_truth"; dist_dir_dca = "dca_resol_truth"; dir = "truthseed";}
+   else {dist_dir_mom = "mom_resol_realseed"; dist_dir_dca = "dca_resol_realseed"; dir = "realseed";}
+
+   bool debug=true;	  
+  // Tree with reconstructed tracks
+   const int nbins_eta = 5;
+   int theta_val[nbins_eta+1] ={3,50,45,135,130,177};
+   int nfiles = 100; 
+   double eta[nbins_eta+1]={-3.5,-2.5,-1.0,1.0,2.5,3.5};
+   double pt[nbins_eta+1]={0.1,0.5,1.0,2.0,5.0,10.0};
+   TH1D *histp[nbins_eta]; 
+
+    TH3D *h_d0xy_3d= new TH3D("h_d0xy_3d","Transverse Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,200,0.,20.);
+    TH3D *h_d0z_3d= new TH3D("h_d0z_3d","Longitudinal Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,200,0.,20.);
+
+   for (int i=0; i<nbins_eta; i++){
+   histp[i] = new TH1D(Form("hist_etabin%d",i),Form("hist_etabin%d",i),600,-0.3,0.3);
+   histp[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f ",eta[i],eta[i+1],mom));
+   histp[i]->SetName(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
+   }
+
+   TFile* file = TFile::Open(Form("./%s_%1.1f.edm4eic.root",filename.Data(),mom));
+   if (!file) {printf("file not found !!!"); return;}
+   TTreeReader myReader("events", file); // name of tree and file
+   if (debug) cout<<"Filename: "<<file->GetName()<<"\t NEvents: "<<myReader.GetEntries()<<endl;
+
+   // MC and Reco information 
+   TTreeReaderArray<Float_t> charge(myReader, "MCParticles.charge"); 
+   TTreeReaderArray<Double_t> vx_mc(myReader, "MCParticles.vertex.x"); 
+   TTreeReaderArray<Double_t> vy_mc(myReader, "MCParticles.vertex.y"); 
+   TTreeReaderArray<Double_t> vz_mc(myReader, "MCParticles.vertex.z"); 
+   TTreeReaderArray<Float_t> px_mc(myReader, "MCParticles.momentum.x"); 
+   TTreeReaderArray<Float_t> py_mc(myReader, "MCParticles.momentum.y"); 
+   TTreeReaderArray<Float_t> pz_mc(myReader, "MCParticles.momentum.z"); 
+   TTreeReaderArray<Int_t> status(myReader, "MCParticles.generatorStatus"); 
+   TTreeReaderArray<Int_t> pdg(myReader, "MCParticles.PDG"); 
+   TTreeReaderArray<Int_t> match_flag(myReader, Form("CentralCKF%sTrackParameters.type",name.Data()));
+   TTreeReaderArray<Float_t> d0xy(myReader, Form("CentralCKF%sTrackParameters.loc.a",name.Data())); 
+   TTreeReaderArray<Float_t> d0z(myReader, Form("CentralCKF%sTrackParameters.loc.b",name.Data()));    
+   TTreeReaderArray<Float_t> theta(myReader, Form("CentralCKF%sTrackParameters.theta",name.Data()));  
+   TTreeReaderArray<Float_t> phi(myReader, Form("CentralCKF%sTrackParameters.phi",name.Data()));
+   TTreeReaderArray<Float_t> qoverp(myReader, Form("CentralCKF%sTrackParameters.qOverP",name.Data()));  
+
+  int count =0;
+  int matchId = 1; // Always matched track assigned the index 0 
+  while (myReader.Next()) 
+  {
+   if (match_flag.GetSize()==0) continue;  // Events with no reco tracks skip them
+   for (int i = 0; i < matchId; ++i){
+
+     for (int j = 0; j < pdg.GetSize(); ++j){
+
+      if (status[j] !=1 && pdg.GetSize()!=1) continue;
+      Double_t ptmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]); 
+
+      if (fabs(ptmc) < pTcut) continue;
+
+      Double_t pmc = charge[j]*sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); 
+      Double_t prec = 1./qoverp[j]; 
+
+      Double_t pzrec = prec*TMath::Cos(theta[j]);  Double_t pt_rec = sqrt(prec*prec-pzrec*pzrec);  
+      Double_t pzmc = pz_mc[j];  
+
+      Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2));
+      Double_t p_resol = (prec-pmc)/pmc;
+
+      for (int ibin=0; ibin<nbins_eta; ++ibin){ 
+      if(etamc>eta[ibin] && etamc<eta[ibin+1]) histp[ibin]->Fill(p_resol); 
+      }
+      h_d0xy_3d->Fill(d0xy[j]*0.1, etamc, ptmc); // cm
+      h_d0z_3d->Fill(d0z[j]*0.1, etamc, ptmc); // cm
+      } // Generated Tracks  
+     } // Reco Tracks
+
+   }// event loop ends    
+
+   TFile *fout_mom = new TFile(Form("%s/%s/mom/Performances_mom_%1.1f_%s_%s.root",dir.Data(),particle.Data(),mom,dist_dir_mom.Data(),particle.Data()),"recreate");
+   fout_mom->cd();
+   for (int ibin=0; ibin<nbins_eta; ++ibin) histp[ibin]->Write();
+   fout_mom->Close();
+
+   TFile *fout_dca = new TFile(Form("%s/%s/dca/Performances_dca_%1.1f_%s_%s.root",dir.Data(),particle.Data(),mom,dist_dir_dca.Data(),particle.Data()),"recreate");
+   fout_dca->cd();
+    h_d0xy_3d->Write();
+    h_d0z_3d->Write();
+    fout_dca->Close();
+}
+
+
+
+
diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C
@@ -0,0 +1,193 @@
+// Code to compare the tracking performances: Truth seeding vs real seeding
+// Shyam Kumar; [email protected]; [email protected]
+
+#include "TGraphErrors.h"
+#include "TF1.h"
+#include "TRandom.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TMath.h"
+#define mpi 0.139  // 1.864 GeV/c^2
+
+void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.);
+void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, double range =0.3, Bool_t drawreq=1) // name = p, pt for getting p or pt dependence fitted results
+{
+
+//=== style of the plot=========
+   gStyle->SetPalette(1);
+   gStyle->SetOptTitle(1);
+   gStyle->SetTitleOffset(1.0,"XY");
+   gStyle->SetTitleSize(.04,"XY");
+   gStyle->SetLabelSize(.04,"XY");
+   gStyle->SetHistLineWidth(2);
+   gStyle->SetOptFit(1);
+   gStyle->SetOptStat(1);
+
+   const Int_t nfiles = 6;
+   double mom[nfiles] ={0.5,1.0,2.0,5.0,10.0,15.0};
+   std::vector<double> momV_truth, momV_real, momresolV_truth, err_momresolV_truth, momresolV_real, err_momresolV_real;
+   momV_truth.clear(); momV_real.clear(); momresolV_truth.clear(); err_momresolV_truth.clear(); momresolV_real.clear(); err_momresolV_real.clear();
+
+   TF1 *f1=new TF1("f1","FitMomentumResolution",0.,30.0,2);
+   f1->SetParLimits(0,0.,0.1);	
+   f1->SetParLimits(1,0.,5.0);	
+
+   TCanvas *c_mom = new TCanvas("cmom","cmom",1400,1000);
+   c_mom->SetMargin(0.10, 0.05 ,0.1,0.05);
+   c_mom->SetGridy();
+
+ //Reading the root file
+  TFile *fmom_truth[nfiles], *fmom_real[nfiles];
+  TGraphErrors *gr_mom_truth, *gr_mom_real;
+
+  TMultiGraph *mgMom; 
+  TLegend *lmom; 
+  mgMom = new TMultiGraph("mgMom",";p (GeV/c); #sigmap/p %");
+
+  lmom = new TLegend(0.70,0.80,0.90,0.93);
+  lmom->SetTextSize(0.03);
+  lmom->SetBorderSize(0);
+  lmom->SetHeader(Form("Particle (%s): %1.1f < #eta < %1.1f",particle.Data(),etamin,etamax),"C");
+
+  TF1 *func_truth = new TF1("func_truth","gaus",-0.5,0.5);
+  TF1 *func_real = new TF1("func_real","gaus",-0.5,0.5);
+
+  for (int i =0; i<nfiles; ++i){
+
+   TCanvas *cp = new TCanvas("cp","cp",1400,1000);
+   cp->SetMargin(0.10, 0.05 ,0.1,0.07);
+
+	 fmom_truth[i] = TFile::Open(Form("./truthseed/pi-/mom/Performances_mom_%1.1f_mom_resol_truth_%s.root",mom[i],particle.Data()));
+	 fmom_real[i] = TFile::Open(Form("./realseed/pi-/mom/Performances_mom_%1.1f_mom_resol_realseed_%s.root",mom[i],particle.Data()));
+
+	 TH1D *hist_truth = (TH1D*) fmom_truth[i]->Get(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom[i],etamin,etamax));
+	 hist_truth->Rebin(2);
+	 hist_truth->SetName(Form("truthseed_hist_mom_%1.1f_%1.1f_eta_%1.1f_%s",mom[i],etamin,etamax,particle.Data()));
+	 hist_truth->SetTitle(Form("Truth Seed (%s): Momentum = %1.1f && %1.1f<#eta<%1.1f;#Delta p/p; Entries(a.u.)",particle.Data(),mom[i],etamin,etamax));
+
+	 double mu_truth = hist_truth->GetMean(); 
+	 double sigma_truth = hist_truth->GetStdDev();
+         hist_truth->GetXaxis()->SetRangeUser(-1.0*range,1.0*range);
+	 func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range
+	 hist_truth->Fit(func_truth,"NR+");
+	 hist_truth->Fit(func_truth,"R+");
+	 float truth_par2 = func_truth->GetParameter(2)*100;
+	 float truth_par2_err = func_truth->GetParError(2)*100;
+	 momV_truth.push_back(mom[i]);
+	 momresolV_truth.push_back(truth_par2);
+	 err_momresolV_truth.push_back(truth_par2_err);
+
+	 TH1D *hist_real = (TH1D*) fmom_real[i]->Get(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom[i],etamin,etamax));
+	 hist_real->Rebin(2);
+	 hist_real->SetName(Form("realseed_hist_mom_%1.1f_%1.1f_eta_%1.1f_%s",mom[i],etamin,etamax,particle.Data()));
+	 hist_real->SetTitle(Form("Realistic Seed (%s): Momentum = %1.1f && %1.1f<#eta<%1.1f;#Delta p/p; Entries(a.u.)",particle.Data(),mom[i],etamin,etamax));
+
+	 double mu_real = hist_real->GetMean(); 
+	 double sigma_real = hist_real->GetStdDev();
+         hist_real->GetXaxis()->SetRangeUser(-1.0*range,1.0*range);
+         func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real); // fit with in 2 sigma range
+	 hist_real->Fit(func_real,"NR+");
+	 mu_real = func_real->GetParameter(1); 
+	 sigma_real = func_real->GetParameter(2);
+	 func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real);
+	 hist_real->Fit(func_real,"R+");
+	 float real_par2 = func_real->GetParameter(2)*100;
+	 float real_par2_err = func_real->GetParError(2)*100;
+	 momV_real.push_back(mom[i]);
+	 momresolV_real.push_back(real_par2);
+	 err_momresolV_real.push_back(real_par2_err);
+	 cp->cd();
+	 hist_truth->Draw();
+	 cp->SaveAs(Form("Debug_Plots/truth/%s/mom/truth_mom_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),mom[i],etamin,etamax));
+	 cp->Clear();
+	 cp->cd();
+	 hist_real->Draw();
+	 cp->SaveAs(Form("Debug_Plots/real/%s/mom/real_mom_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),mom[i],etamin,etamax));
+ } // all files
+
+	const int size_truth = momV_truth.size();
+	double p_truth[size_truth], err_p_truth[size_truth], sigma_p_truth[size_truth], err_sigma_p_truth[size_truth]; 
+
+	for (int i=0; i<size_truth; i++){
+	p_truth[i] = momV_truth.at(i);
+	sigma_p_truth[i] = momresolV_truth.at(i);
+	err_sigma_p_truth[i] = err_momresolV_truth.at(i);
+	err_p_truth[i] = 0.;
+	}
+
+     const int size_real = momV_real.size();
+	double p_real[size_real], err_p_real[size_real], sigma_p_real[size_real], err_sigma_p_real[size_real]; 
+
+	for (int i=0; i<size_real; i++){
+	p_real[i] = momV_real.at(i);
+  sigma_p_real[i] = momresolV_real.at(i);
+	err_sigma_p_real[i] = err_momresolV_real.at(i);
+	err_p_real[i] = 0.;
+	}
+
+  TFile *fout = new TFile(Form("Final_Results/%s/mom/mom_resol_%1.1f_eta_%1.1f.root",particle.Data(),etamin,etamax),"recreate");
+	TGraphErrors *gr1 = new TGraphErrors(size_truth,p_truth,sigma_p_truth,err_p_truth,err_sigma_p_truth);
+        gr1->SetName("gr_truthseed");
+	gr1->SetMarkerStyle(25);
+	gr1->SetMarkerColor(kBlue);
+	gr1->SetMarkerSize(2.0);
+	gr1->SetTitle(";p (GeV/c);#sigmap/p");
+	gr1->GetXaxis()->CenterTitle();
+	gr1->GetYaxis()->CenterTitle();
+
+       TGraphErrors *gr2 = new TGraphErrors(size_real,p_real,sigma_p_real,err_p_real,err_sigma_p_real);
+       gr2->SetName("gr_realseed");
+	gr2->SetMarkerStyle(34);
+	gr2->SetMarkerColor(kRed);
+	gr2->SetMarkerSize(2.0);
+	gr2->SetTitle(";p (GeV/c);#sigmap/p");
+	gr2->GetXaxis()->CenterTitle();
+	gr2->GetYaxis()->CenterTitle();
+
+	mgMom->Add(gr1);
+	mgMom->Add(gr2);
+	c_mom->cd();
+	mgMom->GetXaxis()->SetRangeUser(0.40,10.2);
+	mgMom->GetYaxis()->SetRangeUser(0.,10.0);
+	mgMom->Draw("AP");
+	lmom->AddEntry(gr1,"Truth Seeding");
+	lmom->AddEntry(gr2,"Realistic Seeding");
+	lmom->Draw("same");
+	draw_req_Mom(etamin,etamax,0.,mgMom->GetXaxis()->GetXmax());
+	c_mom->SaveAs(Form("Final_Results/%s/mom/mom_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax));
+	fout->cd();
+	mgMom->SetName(Form("mom_resol_%1.1f_eta_%1.1f",etamin,etamax));
+	mgMom->Write();
+	fout->Close();
+}
+
+//===Fit Momentum Resolution
+float FitMomentumResolution(Double_t *x, Double_t *par)
+{
+  float func = sqrt(par[0]*par[0]*x[0]*x[0]+par[1]*par[1]);
+  return func;
+}
+
+//From Yellow report from section 11.2.2
+
+void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.)
+{
+
+   TF1 *dd4hep_p;
+   if (etamin >= -3.5 && etamax <= -2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax);
+   else if (etamin >= -2.5 && etamax <= -1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax);
+   else if (etamin >= -1.0 && etamax <= 1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+0.5^2)",xmin,xmax);
+   else if (etamin >= 1.0 && etamax <= 2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax);
+   else if (etamin >= 2.5 && etamax <= 3.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax);
+   else return;
+   dd4hep_p->SetLineStyle(7);
+   dd4hep_p->SetLineColor(kMagenta);
+   dd4hep_p->SetLineWidth(3.0);
+   dd4hep_p->Draw("same");
+
+  TLegend *l= new TLegend(0.70,0.75,0.90,0.80);
+  l->SetTextSize(0.03);
+  l->SetBorderSize(0);
+  l->AddEntry(dd4hep_p,"PWGReq","l");
+  l->Draw("same");
+ }
diff --git a/benchmarks/tracking_performances/readme.txt b/benchmarks/tracking_performances/readme.txt
@@ -0,0 +1,6 @@
+Shyam Kumar; INFN Bari, Italy; [email protected]
+Method to produce the tracking performances with ePIC tracker
+The scripts can be used to create the debug plots for the momentum resolutions.
+In the first step enter inside eic-shell and create a directory (your preferred name) there and put all the script inside that.
+Simply run the command source Script_widebin.sh
+It will produce the results with truth/realistic seeding.