Skip to content
Snippets Groups Projects
Unverified Commit 0a85ff9d authored by Shyam Kumar's avatar Shyam Kumar Committed by GitHub
Browse files

Merge pull request #27 from eic/tracking_benchmark 

Tracking benchmark
parents c93e12d1 0d0c9c65
No related branches found
No related tags found
No related merge requests found
Pipeline #96779 failed
Stage: status-pending
Stage: config
Stage: initialize
Stage: data_init
Stage: simulate
Stage: calibrate
Stage: benchmarks
Stage: collect
Stage: deploy
Stage: status-report
Hide whitespace changes
Inline Side-by-side
benchmarks/tracking_performances/Script_widebin.sh 0 → 100644
+ 64
− 0
View file @ 0a85ff9d
#!/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
benchmarks/tracking_performances/Tracking_Performances.C 0 → 100644
+ 118
− 0
View file @ 0a85ff9d
// 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; 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();
}
benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C 0 → 100644
+ 193
− 0
View file @ 0a85ff9d
// 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<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");
}
benchmarks/tracking_performances/readme.txt 0 → 100644
+ 6
− 0
View file @ 0a85ff9d
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.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment