diff --git a/benchmarks/tracking_performances/Script_widebin.sh b/benchmarks/tracking_performances/Script_widebin.sh new file mode 100644 index 00000000..e6335920 --- /dev/null +++ b/benchmarks/tracking_performances/Script_widebin.sh @@ -0,0 +1,64 @@ +#!/bin/bash +# Script; shyam kumar; INFN Bari, Italy +# shyam.kumar@ba.infn.it; shyam055119@gmail.com +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 new file mode 100644 index 00000000..0122406d --- /dev/null +++ b/benchmarks/tracking_performances/Tracking_Performances.C @@ -0,0 +1,118 @@ +// Code to extract the Tracking Performances +// Shyam Kumar; INFN Bari, Italy +// shyam.kumar@ba.infn.it; shyam.kumar@cern.ch + +#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; iSetTitle(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: "<GetName()<<"\t NEvents: "< charge(myReader, "MCParticles.charge"); + TTreeReaderArray vx_mc(myReader, "MCParticles.vertex.x"); + TTreeReaderArray vy_mc(myReader, "MCParticles.vertex.y"); + TTreeReaderArray vz_mc(myReader, "MCParticles.vertex.z"); + TTreeReaderArray px_mc(myReader, "MCParticles.momentum.x"); + TTreeReaderArray py_mc(myReader, "MCParticles.momentum.y"); + TTreeReaderArray pz_mc(myReader, "MCParticles.momentum.z"); + TTreeReaderArray status(myReader, "MCParticles.generatorStatus"); + TTreeReaderArray pdg(myReader, "MCParticles.PDG"); + TTreeReaderArray match_flag(myReader, Form("CentralCKF%sTrackParameters.type",name.Data())); + TTreeReaderArray d0xy(myReader, Form("CentralCKF%sTrackParameters.loc.a",name.Data())); + TTreeReaderArray d0z(myReader, Form("CentralCKF%sTrackParameters.loc.b",name.Data())); + TTreeReaderArray theta(myReader, Form("CentralCKF%sTrackParameters.theta",name.Data())); + TTreeReaderArray phi(myReader, Form("CentralCKF%sTrackParameters.phi",name.Data())); + TTreeReaderArray 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; ibineta[ibin] && etamcFill(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; ibinWrite(); + 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 new file mode 100644 index 00000000..a3871f28 --- /dev/null +++ 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; shyam.kumar@ba.infn.it; shyam055119@gmail.com + +#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 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; iSetMargin(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; iSetName("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 new file mode 100644 index 00000000..4d9badac --- /dev/null +++ b/benchmarks/tracking_performances/readme.txt @@ -0,0 +1,6 @@ +Shyam Kumar; INFN Bari, Italy; shyam055119@gmail.com +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.