From bb39a6fc562f0834abb8400986d22d30417e0736 Mon Sep 17 00:00:00 2001 From: Shyam Date: Wed, 17 Jul 2024 15:32:36 +0200 Subject: [PATCH 01/19] included software version --- .../doCompare_truth_real_widebins_mom.C | 22 ++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C index 4d417428..257960d9 100644 --- a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C @@ -10,7 +10,7 @@ #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 +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, TString epic ="24.06.0", TString eicrecon = "v1.14.0") // name = p, pt for getting p or pt dependence fitted results { //=== style of the plot========= @@ -27,6 +27,11 @@ void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1 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(); + TString symbolname = ""; + if (particle == "pi-") symbolname = "#pi^{-}"; + else symbolname = particle; + ofstream outfile; + outfile.open ("Mom_resol.txt",ios_base::app); TF1 *f1=new TF1("f1","FitMomentumResolution",0.,30.0,2); f1->SetParLimits(0,0.,0.1); @@ -44,10 +49,10 @@ void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1 TLegend *lmom; mgMom = new TMultiGraph("mgMom",";p (GeV/c); #sigmap/p %"); - lmom = new TLegend(0.70,0.80,0.90,0.93); + lmom = new TLegend(0.65,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"); + lmom->SetHeader(Form("%s ePIC(%s/%s): %1.1f < #eta < %1.1f",symbolname.Data(),epic.Data(),eicrecon.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); @@ -158,6 +163,17 @@ void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1 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)); + + // Write the numbers in output file for comparisons + outfile<<"ePIC"<GetN(); ++i){ + double x,ytrue, yreal; + gr1->GetPoint(i,x,ytrue); gr2->GetPoint(i,x,yreal); + outfile<cd(); mgMom->SetName(Form("mom_resol_%1.1f_eta_%1.1f",etamin,etamax)); mgMom->Write(); From fae69af25d4d8495646415481b190027192fe3d1 Mon Sep 17 00:00:00 2001 From: Shyam Date: Wed, 17 Jul 2024 16:35:37 +0200 Subject: [PATCH 02/19] DCAT resolution code --- .../doCompare_truth_real_widebins_dcaT.C | 224 ++++++++++++++++++ 1 file changed, 224 insertions(+) create mode 100755 benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C new file mode 100755 index 00000000..d85a73f3 --- /dev/null +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C @@ -0,0 +1,224 @@ +// 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" + +void draw_req_DCA(double etamin, double etamax, double xmin=0., double xmax=0.); +void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, TString epic ="24.06.0", TString eicrecon = "v1.14.0", 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 nptbins = 6; + double pt[nptbins] ={0.5,1.0,2.0,5.0,10.0,15.0}; + Double_t variation = 0.1; // 10 % variation + std::vector momV_truth, momV_real; + std::vector dcaTresolV_truth, err_dcaTresolV_truth, dcaTresolV_real, err_dcaTresolV_real; + momV_truth.clear(); momV_real.clear(); dcaTresolV_truth.clear(); err_dcaTresolV_truth.clear(); dcaTresolV_real.clear(); err_dcaTresolV_real.clear(); + + TString symbolname = ""; + if (particle == "pi-") symbolname = "#pi^{-}"; + else symbolname = particle; + // Write the parameters in text file (for comparison) + ofstream outfile; + outfile.open ("DCAT_resol.txt",ios_base::app); + + TF1 *f1=new TF1("f1","FitPointingAngle",0.,30.0,2); + f1->SetParLimits(0,0.,50000); + f1->SetParLimits(1,0.,50000); + + + TCanvas *c_dcaxy = new TCanvas("c_dcaxy","c_dcaxy",1400,1000); + c_dcaxy->SetMargin(0.10, 0.05 ,0.1,0.05); + c_dcaxy->SetGridy(); + + + //Reading the root file + TFile *fDCA_truth, *fDCA_real; + TGraphErrors *gr_dcaT_truth, *gr_dcaT_real, *gr_dcaZ_truth, *gr_dcaZ_real; + TMultiGraph *mgDCAT; + TLegend *lDCAT; + + mgDCAT = new TMultiGraph("mgDCAT",";p_{T} (GeV/c); #sigma_{DCA_{T}} (#mum)"); + lDCAT = new TLegend(0.65,0.80,0.90,0.93); + lDCAT->SetTextSize(0.03); + lDCAT->SetBorderSize(0); + lDCAT->SetHeader(Form("%s ePIC(%s/%s): %1.1f < #eta < %1.1f",symbolname.Data(),epic.Data(),eicrecon.Data(),etamin,etamax),"C"); + + fDCA_truth = TFile::Open(Form("./truthseed/%s/dca/final_hist_dca_truthseed.root",particle.Data())); + fDCA_real = TFile::Open(Form("./realseed/%s/dca/final_hist_dca_realseed.root",particle.Data())); + + // Truth seeding histograms + TH3D *hist_d0xy_truth = (TH3D*) fDCA_truth->Get("h_d0xy_3d"); + TH3D *hist_d0xy_real = (TH3D*) fDCA_real->Get("h_d0xy_3d"); + + // d0xy calculation for truth/real (binning are same in both cases) + Int_t etamin_bin = hist_d0xy_truth->GetYaxis()->FindBin(etamin+0.0001); + Int_t etamax_bin = hist_d0xy_truth->GetYaxis()->FindBin(etamax-0.0001); + + 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 iptbin=0; iptbinSetMargin(0.10, 0.05 ,0.1,0.07); + + double ptmin = (1.0-variation)*pt[iptbin]; // 10% range + double ptmax = (1.0+variation)*pt[iptbin]; // 10% range + + Int_t ptmin_bin = hist_d0xy_truth->GetZaxis()->FindBin(ptmin+0.0001); + Int_t ptmax_bin = hist_d0xy_truth->GetZaxis()->FindBin(ptmax-0.0001); + + TH1D *histd0xy_truth_1d = (TH1D*)hist_d0xy_truth->ProjectionX(Form("histd0xy_truth_eta%1.1f_%1.1f_pt%1.1f_%1.1f",etamin,etamax,ptmin,ptmax),etamin_bin,etamax_bin,ptmin_bin,ptmax_bin,"o"); + histd0xy_truth_1d->SetTitle(Form("d0_{xy} (truth): %1.1f <#eta< %1.1f && %1.2f SetName(Form("eta_%1.1f_%1.1f_d0xy_truth_pt_%1.1f",etamin,etamax,pt[iptbin])); + + double mu_truth = histd0xy_truth_1d->GetMean(); + double sigma_truth = histd0xy_truth_1d->GetStdDev(); + func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range + histd0xy_truth_1d->Fit(func_truth,"NR+"); + mu_truth = func_truth->GetParameter(1); + sigma_truth = func_truth->GetParameter(2); + func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range + histd0xy_truth_1d->Fit(func_truth,"R+"); + float truth_par2 = func_truth->GetParameter(2)*10000; // cm to mum 10000 factor + float truth_par2_err = func_truth->GetParError(2)*10000; + momV_truth.push_back(pt[iptbin]); + dcaTresolV_truth.push_back(truth_par2); + err_dcaTresolV_truth.push_back(truth_par2_err); + + TH1D *histd0xy_real_1d = (TH1D*)hist_d0xy_real->ProjectionX(Form("histd0xy_real_eta%1.1f_%1.1f_pt%1.1f_%1.1f",etamin,etamax,ptmin,ptmax),etamin_bin,etamax_bin,ptmin_bin,ptmax_bin,"o"); + histd0xy_real_1d->SetTitle(Form("d0_{xy} (real): %1.1f <#eta< %1.1f && %1.2f SetName(Form("eta_%1.1f_%1.1f_d0xy_real_pt_%1.1f",etamin,etamax,pt[iptbin])); + + double mu_real = histd0xy_real_1d->GetMean(); + double sigma_real = histd0xy_real_1d->GetStdDev(); + func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real); // fit with in 2 sigma range + histd0xy_real_1d->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); // fit with in 2 sigma range + histd0xy_real_1d->Fit(func_real,"R+"); + float real_par2 = func_real->GetParameter(2)*10000; // cm to mum 10000 factor + float real_par2_err = func_real->GetParError(2)*10000; + momV_real.push_back(pt[iptbin]); + dcaTresolV_real.push_back(real_par2); + err_dcaTresolV_real.push_back(real_par2_err); + + cp->cd(); + histd0xy_truth_1d->Draw(); + cp->SaveAs(Form("Debug_Plots/truth/%s/dca/truth_dcaxy_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),pt[iptbin],etamin,etamax)); + cp->Clear(); + cp->cd(); + histd0xy_real_1d->Draw(); + cp->SaveAs(Form("Debug_Plots/real/%s/dca/real_dcaxy_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),pt[iptbin],etamin,etamax)); + } // ptbin + + const int size_truth = momV_truth.size(); + double pt_truth[size_truth], err_pt_truth[size_truth], sigma_dcaxy_truth[size_truth], err_sigma_dcaxy_truth[size_truth]; + + for (int i=0; iSetName("gr_truthseed"); + gr1->SetMarkerStyle(25); + gr1->SetMarkerColor(kMagenta); + gr1->SetMarkerSize(2.0); + gr1->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{T}} (#mum)"); + gr1->GetXaxis()->CenterTitle(); + gr1->GetYaxis()->CenterTitle(); + + TGraphErrors *gr2 = new TGraphErrors(size_real,pt_real,sigma_dcaxy_real,err_pt_real,err_sigma_dcaxy_real); + gr2->SetName("gr_realseed"); + gr2->SetMarkerStyle(34); + gr2->SetMarkerColor(kRed); + gr2->SetMarkerSize(2.0); + gr2->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{T}} (#mum)"); + gr2->GetXaxis()->CenterTitle(); + gr2->GetYaxis()->CenterTitle(); + + mgDCAT->Add(gr1); + mgDCAT->Add(gr2); + c_dcaxy->cd(); +// c_dcaxy->SetLogy(); + mgDCAT->GetYaxis()->SetRangeUser(0.45,mgDCAT->GetXaxis()->GetXmax()); + mgDCAT->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr2->GetN(),gr2->GetY())); // 50% more of the maximum range + mgDCAT->Draw("AP"); + lDCAT->AddEntry(gr1,"Truth Seeding"); + lDCAT->AddEntry(gr2,"Realistic Seeding"); + lDCAT->Draw("same"); + draw_req_DCA(etamin,etamax,0.,mgDCAT->GetXaxis()->GetXmax()); + c_dcaxy->SaveAs(Form("Final_Results/%s/dca/dcaT_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax)); + // Write the numbers in output file for comparisons + outfile<<"ePIC"<GetN(); ++i){ + double x,ytrue, yreal; + gr1->GetPoint(i,x,ytrue); gr2->GetPoint(i,x,yreal); + outfile<cd(); + mgDCAT->SetName(Form("dcaT_resol_%1.1f_eta_%1.1f",etamin,etamax)); + mgDCAT->Write(); + fout->Close(); +} + +//From Yellow report from section 11.2.2 +//===Fit Pointing Resolution +float FitPointingAngle(Double_t *x, Double_t *par) +{ + float func = sqrt((par[0]*par[0])/(x[0]*x[0])+par[1]*par[1]); + return func; +} + +void draw_req_DCA(double etamin, double etamax, double xmin=0., double xmax=0.) +{ + TF1 *PWGReq_DCA2D; + if (etamin >= -3.5 && etamax <= -2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax); + else if (etamin >= -2.5 && etamax <= -1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax); + else if (etamin >= -1.0 && etamax <= 1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((20./x)^2+5.0^2)",xmin,xmax); + else if (etamin >= 1.0 && etamax <= 2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax); + else if (etamin >= 2.5 && etamax <= 3.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax); + else return; + PWGReq_DCA2D->SetLineStyle(7); + PWGReq_DCA2D->SetLineWidth(3.0); + PWGReq_DCA2D->SetLineColor(kBlue); + PWGReq_DCA2D->Draw("same"); + + TLegend *l= new TLegend(0.70,0.75,0.90,0.80); + l->SetTextSize(0.03); + l->SetBorderSize(0); + l->AddEntry(PWGReq_DCA2D,"PWGReq","l"); + l->Draw("same"); +} From 05c666c224fc6af51e4c1f484440a820106497a9 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Wed, 17 Jul 2024 16:46:10 +0200 Subject: [PATCH 03/19] Update doCompare_truth_real_widebins_dcaT.C --- .../tracking_performances/doCompare_truth_real_widebins_dcaT.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C index d85a73f3..89f58081 100755 --- a/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C @@ -9,7 +9,7 @@ #include "TMath.h" void draw_req_DCA(double etamin, double etamax, double xmin=0., double xmax=0.); -void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, TString epic ="24.06.0", TString eicrecon = "v1.14.0", Bool_t drawreq=1) // name = p, pt for getting p or pt dependence fitted results +void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, Bool_t drawreq=1, TString epic ="24.06.0", TString eicrecon = "v1.14.0") // name = p, pt for getting p or pt dependence fitted results { //=== style of the plot========= From 2c0e81add32792fcac6652e76aa5bf11c3e86e5c Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Thu, 18 Jul 2024 19:57:39 +0200 Subject: [PATCH 04/19] Update doCompare_truth_real_widebins_dcaT.C Minor fix in the code for DCAT resolutions --- .../doCompare_truth_real_widebins_dcaT.C | 115 +++++++++--------- 1 file changed, 60 insertions(+), 55 deletions(-) diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C index 89f58081..4b628944 100755 --- a/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C @@ -22,13 +22,14 @@ void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=- gStyle->SetOptFit(1); gStyle->SetOptStat(1); - const Int_t nptbins = 6; - double pt[nptbins] ={0.5,1.0,2.0,5.0,10.0,15.0}; + const Int_t nptbins = 10; + double pt[nptbins] ={0.2, 0.3, 0.5,1.0, 1.5, 2.0, 5.0, 8.0, 10., 15.0}; Double_t variation = 0.1; // 10 % variation std::vector momV_truth, momV_real; std::vector dcaTresolV_truth, err_dcaTresolV_truth, dcaTresolV_real, err_dcaTresolV_real; momV_truth.clear(); momV_real.clear(); dcaTresolV_truth.clear(); err_dcaTresolV_truth.clear(); dcaTresolV_real.clear(); err_dcaTresolV_real.clear(); + TString symbolname = ""; if (particle == "pi-") symbolname = "#pi^{-}"; else symbolname = particle; @@ -36,41 +37,41 @@ void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=- ofstream outfile; outfile.open ("DCAT_resol.txt",ios_base::app); + TF1 *f1=new TF1("f1","FitPointingAngle",0.,30.0,2); f1->SetParLimits(0,0.,50000); f1->SetParLimits(1,0.,50000); - TCanvas *c_dcaxy = new TCanvas("c_dcaxy","c_dcaxy",1400,1000); - c_dcaxy->SetMargin(0.10, 0.05 ,0.1,0.05); - c_dcaxy->SetGridy(); - + TCanvas *c_dcaxy = new TCanvas("c_dcaxy","c_dcaxy",1400,1000); + c_dcaxy->SetMargin(0.10, 0.05 ,0.1,0.05); + c_dcaxy->SetGridy(); - //Reading the root file - TFile *fDCA_truth, *fDCA_real; - TGraphErrors *gr_dcaT_truth, *gr_dcaT_real, *gr_dcaZ_truth, *gr_dcaZ_real; + //Reading the root file + TFile *fDCA_truth, *fDCA_real; + TGraphErrors *gr_dcaT_truth, *gr_dcaT_real, *gr_dcaZ_truth, *gr_dcaZ_real; + TMultiGraph *mgDCAT; - TLegend *lDCAT; - - mgDCAT = new TMultiGraph("mgDCAT",";p_{T} (GeV/c); #sigma_{DCA_{T}} (#mum)"); - lDCAT = new TLegend(0.65,0.80,0.90,0.93); - lDCAT->SetTextSize(0.03); - lDCAT->SetBorderSize(0); - lDCAT->SetHeader(Form("%s ePIC(%s/%s): %1.1f < #eta < %1.1f",symbolname.Data(),epic.Data(),eicrecon.Data(),etamin,etamax),"C"); + TLegend *lDCAT; + mgDCAT = new TMultiGraph("mgDCAT",";p_{T} (GeV/c); #sigma_{DCA_{T}} (#mum)"); + lDCAT = new TLegend(0.65,0.80,0.90,0.93); + lDCAT->SetTextSize(0.03); + lDCAT->SetBorderSize(0); + lDCAT->SetHeader(Form("%s ePIC(%s/%s): %1.1f < #eta < %1.1f",symbolname.Data(),epic.Data(),eicrecon.Data(),etamin,etamax),"C"); - fDCA_truth = TFile::Open(Form("./truthseed/%s/dca/final_hist_dca_truthseed.root",particle.Data())); - fDCA_real = TFile::Open(Form("./realseed/%s/dca/final_hist_dca_realseed.root",particle.Data())); + fDCA_truth = TFile::Open(Form("./truthseed/%s/dca/final_hist_dca_truthseed.root",particle.Data())); + fDCA_real = TFile::Open(Form("./realseed/%s/dca/final_hist_dca_realseed.root",particle.Data())); // Truth seeding histograms - TH3D *hist_d0xy_truth = (TH3D*) fDCA_truth->Get("h_d0xy_3d"); - TH3D *hist_d0xy_real = (TH3D*) fDCA_real->Get("h_d0xy_3d"); + TH3D *hist_d0xy_truth = (TH3D*) fDCA_truth->Get("h_d0xy_3d"); + TH3D *hist_d0xy_real = (TH3D*) fDCA_real->Get("h_d0xy_3d"); // d0xy calculation for truth/real (binning are same in both cases) - Int_t etamin_bin = hist_d0xy_truth->GetYaxis()->FindBin(etamin+0.0001); - Int_t etamax_bin = hist_d0xy_truth->GetYaxis()->FindBin(etamax-0.0001); + Int_t etamin_bin = hist_d0xy_truth->GetYaxis()->FindBin(etamin+0.0001); + Int_t etamax_bin = hist_d0xy_truth->GetYaxis()->FindBin(etamax-0.0001); - TF1 *func_truth = new TF1("func_truth","gaus",-0.5,0.5); - TF1 *func_real = new TF1("func_real","gaus",-0.5,0.5); + 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 iptbin=0; iptbinSetTitle(Form("d0_{xy} (truth): %1.1f <#eta< %1.1f && %1.2f SetName(Form("eta_%1.1f_%1.1f_d0xy_truth_pt_%1.1f",etamin,etamax,pt[iptbin])); + if (histd0xy_truth_1d->GetEntries()<100) continue; double mu_truth = histd0xy_truth_1d->GetMean(); double sigma_truth = histd0xy_truth_1d->GetStdDev(); func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range histd0xy_truth_1d->Fit(func_truth,"NR+"); - mu_truth = func_truth->GetParameter(1); - sigma_truth = func_truth->GetParameter(2); + mu_truth = func_truth->GetParameter(2); + sigma_truth = func_truth->GetParError(2); func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range histd0xy_truth_1d->Fit(func_truth,"R+"); float truth_par2 = func_truth->GetParameter(2)*10000; // cm to mum 10000 factor @@ -104,7 +106,8 @@ void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=- TH1D *histd0xy_real_1d = (TH1D*)hist_d0xy_real->ProjectionX(Form("histd0xy_real_eta%1.1f_%1.1f_pt%1.1f_%1.1f",etamin,etamax,ptmin,ptmax),etamin_bin,etamax_bin,ptmin_bin,ptmax_bin,"o"); histd0xy_real_1d->SetTitle(Form("d0_{xy} (real): %1.1f <#eta< %1.1f && %1.2f SetName(Form("eta_%1.1f_%1.1f_d0xy_real_pt_%1.1f",etamin,etamax,pt[iptbin])); - + + if (histd0xy_real_1d->GetEntries()<100) continue; double mu_real = histd0xy_real_1d->GetMean(); double sigma_real = histd0xy_real_1d->GetStdDev(); func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real); // fit with in 2 sigma range @@ -138,19 +141,19 @@ void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=- err_pt_truth[i] = 0.; } - const int size_real = momV_real.size(); + const int size_real = momV_real.size(); double pt_real[size_real], err_pt_real[size_real], sigma_dcaxy_real[size_real], err_sigma_dcaxy_real[size_real]; for (int i=0; iSetName("gr_truthseed"); + gr1->SetName("gr_truthseed"); gr1->SetMarkerStyle(25); gr1->SetMarkerColor(kMagenta); gr1->SetMarkerSize(2.0); @@ -158,8 +161,8 @@ void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=- gr1->GetXaxis()->CenterTitle(); gr1->GetYaxis()->CenterTitle(); - TGraphErrors *gr2 = new TGraphErrors(size_real,pt_real,sigma_dcaxy_real,err_pt_real,err_sigma_dcaxy_real); - gr2->SetName("gr_realseed"); + TGraphErrors *gr2 = new TGraphErrors(size_real,pt_real,sigma_dcaxy_real,err_pt_real,err_sigma_dcaxy_real); + gr2->SetName("gr_realseed"); gr2->SetMarkerStyle(34); gr2->SetMarkerColor(kRed); gr2->SetMarkerSize(2.0); @@ -167,31 +170,33 @@ void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=- gr2->GetXaxis()->CenterTitle(); gr2->GetYaxis()->CenterTitle(); - mgDCAT->Add(gr1); +// mgDCAT->Add(gr1); mgDCAT->Add(gr2); c_dcaxy->cd(); // c_dcaxy->SetLogy(); - mgDCAT->GetYaxis()->SetRangeUser(0.45,mgDCAT->GetXaxis()->GetXmax()); - mgDCAT->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr2->GetN(),gr2->GetY())); // 50% more of the maximum range - mgDCAT->Draw("AP"); - lDCAT->AddEntry(gr1,"Truth Seeding"); - lDCAT->AddEntry(gr2,"Realistic Seeding"); - lDCAT->Draw("same"); - draw_req_DCA(etamin,etamax,0.,mgDCAT->GetXaxis()->GetXmax()); - c_dcaxy->SaveAs(Form("Final_Results/%s/dca/dcaT_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax)); - // Write the numbers in output file for comparisons - outfile<<"ePIC"<GetN(); ++i){ - double x,ytrue, yreal; - gr1->GetPoint(i,x,ytrue); gr2->GetPoint(i,x,yreal); - outfile<cd(); - mgDCAT->SetName(Form("dcaT_resol_%1.1f_eta_%1.1f",etamin,etamax)); - mgDCAT->Write(); - fout->Close(); + mgDCAT->GetXaxis()->SetRangeUser(0.18,mgDCAT->GetXaxis()->GetXmax()); + mgDCAT->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr2->GetN(),gr2->GetY())); + mgDCAT->Draw("AP"); + // lDCAT->AddEntry(gr1,"Truth Seeding"); + lDCAT->AddEntry(gr2,"Realistic Seeding"); + lDCAT->Draw("same"); + draw_req_DCA(etamin,etamax,0.,mgDCAT->GetXaxis()->GetXmax()); + c_dcaxy->SaveAs(Form("Final_Results/%s/dca/dcaxy_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax)); + + // Write the numbers in output file for comparisons + outfile<<"ePIC"<GetN(); ++i){ + double x,ytrue, yreal; + gr2->GetPoint(i,x,yreal); + outfile<cd(); + mgDCAT->SetName(Form("dcaxy_resol_%1.1f_eta_%1.1f",etamin,etamax)); + mgDCAT->Write(); + fout->Close(); } //From Yellow report from section 11.2.2 From cf3f09dac61c075de2bf2d1d5249de66035f8f27 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Thu, 18 Jul 2024 14:46:35 -0400 Subject: [PATCH 05/19] avoid hardcoding versions in doCompare_truth_real_widebins_mom.C --- benchmarks/tracking_performances/Snakefile | 8 +++++++- .../doCompare_truth_real_widebins_mom.C | 4 ++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/benchmarks/tracking_performances/Snakefile b/benchmarks/tracking_performances/Snakefile index 5cbf9278..b3ff5373 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -93,7 +93,13 @@ rule tracking_performance_summary_at_eta: "Final_Results/pi-/mom/mom_resol_{ETA_MIN}_eta_{ETA_MAX}.root", shell: """ -root -l -b -q {input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1.)' +set +e +EPIC_VERSION="$(echo $DETECTOR_PATH | sed -n -e 's/.*epic-\([^-]\+\).*/\\1/p')" +EICRECON_VERSION="$(eicrecon -v | sed -n -e 's/.*\(v[0-9.]\+\).*/\\1/p')" +set -e +echo "ePIC version: $EPIC_VERSION" +echo "EICrecon version: $EICRECON_VERSION" +root -l -b -q {input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'$EPIC_VERSION'", "'$EICRECON_VERSION'")' """ diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C index 257960d9..41230a63 100644 --- a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C @@ -10,7 +10,7 @@ #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, TString epic ="24.06.0", TString eicrecon = "v1.14.0") // name = p, pt for getting p or pt dependence fitted results +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, TString epic ="", TString eicrecon = "") // name = p, pt for getting p or pt dependence fitted results { //=== style of the plot========= @@ -165,7 +165,7 @@ void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1 c_mom->SaveAs(Form("Final_Results/%s/mom/mom_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax)); // Write the numbers in output file for comparisons - outfile<<"ePIC"<GetN(); ++i){ double x,ytrue, yreal; From e526b7d016213a2e44072b05eac751c273dc8cb6 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Fri, 19 Jul 2024 07:23:22 +0200 Subject: [PATCH 06/19] Update doCompare_truth_real_widebins_mom.C Fixing the symbol name Pt-->p. --- .../tracking_performances/doCompare_truth_real_widebins_mom.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C index 41230a63..7b8e799b 100644 --- a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C @@ -166,7 +166,7 @@ void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1 // Write the numbers in output file for comparisons outfile<<"ePIC"<GetN(); ++i){ double x,ytrue, yreal; gr1->GetPoint(i,x,ytrue); gr2->GetPoint(i,x,yreal); From fb8ec8461b8df26470a271a74c340e5da27c821a Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 20 Jul 2024 16:27:13 -0400 Subject: [PATCH 07/19] benchmarks/tracking_performances/README.md --- .../tracking_performances/{readme.txt => README.md} | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) rename benchmarks/tracking_performances/{readme.txt => README.md} (58%) diff --git a/benchmarks/tracking_performances/readme.txt b/benchmarks/tracking_performances/README.md similarity index 58% rename from benchmarks/tracking_performances/readme.txt rename to benchmarks/tracking_performances/README.md index 4d9badac..a5610e63 100644 --- a/benchmarks/tracking_performances/readme.txt +++ b/benchmarks/tracking_performances/README.md @@ -1,6 +1,10 @@ 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 + +To run a full simulation-reconstruction-analysis chain do: +``` +snakemake --cores 1 tracking_performance +``` + It will produce the results with truth/realistic seeding. From 80a8746b2f752bcfa1489066cddbc46ca8224ea3 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 20 Jul 2024 17:50:34 -0400 Subject: [PATCH 08/19] tracking_performances: replace 15 GeV with 20 GeV (supported by campaigns) --- benchmarks/tracking_performances/Snakefile | 26 +++++++++++++------ .../Tracking_Performances.C | 6 ++--- .../doCompare_truth_real_widebins_mom.C | 4 +-- 3 files changed, 23 insertions(+), 13 deletions(-) diff --git a/benchmarks/tracking_performances/Snakefile b/benchmarks/tracking_performances/Snakefile index 5cbf9278..8bbdee90 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -48,12 +48,22 @@ rule tracking_performance_at_momentum: input: script="benchmarks/tracking_performances/Tracking_Performances.C", # TODO pass as a file list? - sim=lambda wildcards: expand( - "sim_output/tracking_performance/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", - DETECTOR_CONFIG="epic_craterlake_tracking_only", - ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", - PHASE_SPACE=["3to50deg", "45to135deg", "130to177deg"], - INDEX=range(1), + sim=lambda wildcards: branch( + wildcards.CAMPAIGN == "local", + then=expand( + "sim_output/tracking_performance/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + DETECTOR_CONFIG="epic_craterlake_tracking_only", + ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", + PHASE_SPACE=["3to50deg", "45to135deg", "130to177deg"], + INDEX=range(1), + ), + otherwise=expand( + "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + DETECTOR_CONFIG="epic_craterlake", + ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", + PHASE_SPACE=["3to50deg", "45to135deg", "130to177deg"], + INDEX=range(1), + ), ), output: "{SEEDING}/pi-/mom/Performances_mom_{MOMENTUM}_mom_resol_{SEEDING_IGNORE}_{PARTICLE}.root", @@ -80,14 +90,14 @@ rule tracking_performance_summary_at_eta: "realseed/pi-/mom/Performances_mom_{MOMENTUM:.1f}_mom_resol_realseed_pi-.root", "realseed/pi-/dca/Performances_dca_{MOMENTUM:.1f}_dca_resol_realseed_pi-.root", ], - MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 15.0], + MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], ), script="benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C", output: expand( "Debug_Plots/{SEEDING}/pi-/mom/{SEEDING}_mom_resol_mom{MOMENTUM:.1f}_{{ETA_MIN}}_eta_{{ETA_MAX}}.png", SEEDING=["real", "truth"], - MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 15.0], + MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], ), "Final_Results/pi-/mom/mom_resol_{ETA_MIN}_eta_{ETA_MAX}.png", "Final_Results/pi-/mom/mom_resol_{ETA_MIN}_eta_{ETA_MAX}.root", diff --git a/benchmarks/tracking_performances/Tracking_Performances.C b/benchmarks/tracking_performances/Tracking_Performances.C index 6d1b213c..bfb751d6 100644 --- a/benchmarks/tracking_performances/Tracking_Performances.C +++ b/benchmarks/tracking_performances/Tracking_Performances.C @@ -34,11 +34,11 @@ void Tracking_Performances(TString filename="tracking_output",TString particle=" 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}; + double pt[nbins_eta+1]={0.5,1.0,2.0,5.0,10.0,20.1}; 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.); + TH3D *h_d0xy_3d= new TH3D("h_d0xy_3d","Transverse Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,201,0.,20.1); + TH3D *h_d0z_3d= new TH3D("h_d0z_3d","Longitudinal Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,201,0.,20.1); for (int i=0; iSetOptStat(1); const Int_t nfiles = 6; - double mom[nfiles] ={0.5,1.0,2.0,5.0,10.0,15.0}; + double mom[nfiles] ={0.5,1.0,2.0,5.0,10.0,20.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(); @@ -150,7 +150,7 @@ void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1 mgMom->Add(gr1); mgMom->Add(gr2); c_mom->cd(); - mgMom->GetXaxis()->SetRangeUser(0.40,15.2); + mgMom->GetXaxis()->SetRangeUser(0.40,20.2); mgMom->GetYaxis()->SetRangeUser(0.,10.0); mgMom->Draw("AP"); lmom->AddEntry(gr1,"Truth Seeding"); From e6cff5fa63bb4746d48fd498db7714ccd5c89ab8 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Mon, 22 Jul 2024 14:23:56 +0200 Subject: [PATCH 09/19] Update doCompare_truth_real_widebins_dcaT.C --- .../tracking_performances/doCompare_truth_real_widebins_dcaT.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C index 4b628944..926a02ae 100755 --- a/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C @@ -185,7 +185,7 @@ void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=- // Write the numbers in output file for comparisons outfile<<"ePIC"<GetN(); ++i){ double x,ytrue, yreal; gr2->GetPoint(i,x,yreal); From 1a2c8fdc9933e4c897911d0b3c903325ecd140cf Mon Sep 17 00:00:00 2001 From: Shyam Date: Mon, 22 Jul 2024 14:26:56 +0200 Subject: [PATCH 10/19] DCAz resolution --- .../doCompare_truth_real_widebins_dcaz.C | 225 ++++++++++++++++++ 1 file changed, 225 insertions(+) create mode 100755 benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaz.C diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaz.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaz.C new file mode 100755 index 00000000..e8fff38c --- /dev/null +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaz.C @@ -0,0 +1,225 @@ +// 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_DCA(double etamin, double etamax, double xmin=0., double xmax=0.); +void doCompare_truth_real_widebins_dcaz(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, Bool_t drawreq=1, TString epic ="", TString eicrecon = "") // 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 nptbins = 10; + double pt[nptbins] ={0.2, 0.3, 0.5,1.0, 1.5, 2.0, 5.0, 8.0, 10., 15.0}; + Double_t variation = 0.1; // 10 % variation + std::vector momV_truth, momV_real; + std::vector dcaZresolV_truth, err_dcaZresolV_truth, dcaZresolV_real, err_dcaZresolV_real; + momV_truth.clear(); momV_real.clear(); dcaZresolV_truth.clear(); err_dcaZresolV_truth.clear(); dcaZresolV_real.clear(); err_dcaZresolV_real.clear(); + + TString symbolname = ""; + if (particle == "pi-") symbolname = "#pi^{-}"; + else symbolname = particle; + // Write the parameters in text file (for comparison) + ofstream outfile; + outfile.open ("DCAZ_resol.txt",ios_base::app); + + TF1 *f1=new TF1("f1","FitPointingAngle",0.,30.0,2); + f1->SetParLimits(0,0.,50000); + f1->SetParLimits(1,0.,50000); + + + TCanvas *c_dcaz = new TCanvas("c_dcaz","c_dcaz",1400,1000); + c_dcaz->SetMargin(0.10, 0.05 ,0.1,0.05); + c_dcaz->SetGridy(); + + //Reading the root file + TFile *fDCA_truth, *fDCA_real; + + TMultiGraph *mgDCAZ; + TLegend *lDCAZ; + mgDCAZ = new TMultiGraph("mgDCAZ",";p_{T} (GeV/c); #sigma_{DCA_{Z}} (#mum)"); + lDCAZ = new TLegend(0.70,0.80,0.90,0.93); + lDCAZ->SetTextSize(0.03); + lDCAZ->SetBorderSize(0); + lDCAZ->SetHeader(Form("%s ePIC(%s/%s): %1.1f < #eta < %1.1f",symbolname.Data(),epic.Data(),eicrecon.Data(),etamin,etamax),"C"); + + fDCA_truth = TFile::Open(Form("./truthseed/%s/dca/final_hist_dca_truthseed.root",particle.Data())); + fDCA_real = TFile::Open(Form("./realseed/%s/dca/final_hist_dca_realseed.root",particle.Data())); + + // Truth seeding histograms + TH3D *hist_d0z_truth = (TH3D*) fDCA_truth->Get("h_d0z_3d"); + TH3D *hist_d0z_real = (TH3D*) fDCA_real->Get("h_d0z_3d"); + + // d0z calculation for truth/real (binning are same in both cases) + Int_t etamin_bin = hist_d0z_truth->GetYaxis()->FindBin(etamin+0.0001); + Int_t etamax_bin = hist_d0z_truth->GetYaxis()->FindBin(etamax-0.0001); + + 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 iptbin=0; iptbinSetMargin(0.10, 0.05 ,0.1,0.07); + + double ptmin = (1.0-variation)*pt[iptbin]; // 10% range + double ptmax = (1.0+variation)*pt[iptbin]; // 10% range + + Int_t ptmin_bin = hist_d0z_truth->GetZaxis()->FindBin(ptmin+0.0001); + Int_t ptmax_bin = hist_d0z_truth->GetZaxis()->FindBin(ptmax-0.0001); + + TH1D *histd0z_truth_1d = (TH1D*)hist_d0z_truth->ProjectionX(Form("histd0z_truth_eta%1.1f_%1.1f_pt%1.1f_%1.1f",etamin,etamax,ptmin,ptmax),etamin_bin,etamax_bin,ptmin_bin,ptmax_bin,"o"); + histd0z_truth_1d->SetTitle(Form("d0_{z} (truth): %1.1f <#eta< %1.1f && %1.2f SetName(Form("eta_%1.1f_%1.1f_d0z_truth_pt_%1.1f",etamin,etamax,pt[iptbin])); + if (histd0z_truth_1d->GetEntries()<100) continue; + double mu_truth = histd0z_truth_1d->GetMean(); + double sigma_truth = histd0z_truth_1d->GetStdDev(); + func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range + histd0z_truth_1d->Fit(func_truth,"NR+"); + mu_truth = func_truth->GetParameter(2); + sigma_truth = func_truth->GetParError(2); + func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range + histd0z_truth_1d->Fit(func_truth,"R+"); + float truth_par2 = func_truth->GetParameter(2)*10000; // cm to mum 10000 factor + float truth_par2_err = func_truth->GetParError(2)*10000; + momV_truth.push_back(pt[iptbin]); + dcaZresolV_truth.push_back(truth_par2); + err_dcaZresolV_truth.push_back(truth_par2_err); + + TH1D *histd0z_real_1d = (TH1D*)hist_d0z_real->ProjectionX(Form("histd0z_real_eta%1.1f_%1.1f_pt%1.1f_%1.1f",etamin,etamax,ptmin,ptmax),etamin_bin,etamax_bin,ptmin_bin,ptmax_bin,"o"); + histd0z_real_1d->SetTitle(Form("d0_{z} (real): %1.1f <#eta< %1.1f && %1.2f SetName(Form("eta_%1.1f_%1.1f_d0z_real_pt_%1.1f",etamin,etamax,pt[iptbin])); + + if (histd0z_real_1d->GetEntries()<100) continue; + double mu_real = histd0z_real_1d->GetMean(); + double sigma_real = histd0z_real_1d->GetStdDev(); + func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real); // fit with in 2 sigma range + histd0z_real_1d->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); // fit with in 2 sigma range + histd0z_real_1d->Fit(func_real,"R+"); + float real_par2 = func_real->GetParameter(2)*10000; // cm to mum 10000 factor + float real_par2_err = func_real->GetParError(2)*10000; + momV_real.push_back(pt[iptbin]); + dcaZresolV_real.push_back(real_par2); + err_dcaZresolV_real.push_back(real_par2_err); + + cp->cd(); + histd0z_truth_1d->Draw(); + cp->SaveAs(Form("Debug_Plots/truth/%s/dca/truth_dcaz_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),pt[iptbin],etamin,etamax)); + cp->Clear(); + cp->cd(); + histd0z_real_1d->Draw(); + cp->SaveAs(Form("Debug_Plots/real/%s/dca/real_dcaz_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),pt[iptbin],etamin,etamax)); + } // ptbin + + const int size_truth = momV_truth.size(); + double pt_truth[size_truth], err_pt_truth[size_truth], sigma_dcaz_truth[size_truth], err_sigma_dcaz_truth[size_truth]; + + for (int i=0; iSetName("gr_truthseed"); + gr1->SetMarkerStyle(25); + gr1->SetMarkerColor(kMagenta); + gr1->SetMarkerSize(2.0); + gr1->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{Z}} (#mum)"); + gr1->GetXaxis()->CenterTitle(); + gr1->GetYaxis()->CenterTitle(); + + TGraphErrors *gr2 = new TGraphErrors(size_real,pt_real,sigma_dcaz_real,err_pt_real,err_sigma_dcaz_real); + gr2->SetName("gr_realseed"); + gr2->SetMarkerStyle(34); + gr2->SetMarkerColor(kRed); + gr2->SetMarkerSize(2.0); + gr2->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{Z}} (#mum)"); + gr2->GetXaxis()->CenterTitle(); + gr2->GetYaxis()->CenterTitle(); + + //mgDCAZ->Add(gr1); + mgDCAZ->Add(gr2); + c_dcaz->cd(); + //c_dcaz->SetLogy(); + mgDCAZ->GetXaxis()->SetRangeUser(0.18,mgDCAZ->GetXaxis()->GetXmax()); + mgDCAZ->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr2->GetN(),gr2->GetY())); + mgDCAZ->Draw("AP"); + // lDCAZ->AddEntry(gr1,"Truth Seeding"); + lDCAZ->AddEntry(gr2,"Realistic Seeding"); + lDCAZ->Draw("same"); + // draw_req_DCA(etamin,etamax,0.,mgDCAZ->GetXaxis()->GetXmax()); + c_dcaz->SaveAs(Form("Final_Results/%s/dca/dcaz_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax)); + + // Write the numbers in output file for comparisons + outfile<<"ePIC"<GetN(); ++i){ + double x,ytrue, yreal; + gr2->GetPoint(i,x,yreal); + outfile<cd(); + mgDCAZ->SetName(Form("dcaz_resol_%1.1f_eta_%1.1f",etamin,etamax)); + mgDCAZ->Write(); + fout->Close(); +} + +//From Yellow report from section 11.2.2 +//===Fit Pointing Resolution +float FitPointingAngle(Double_t *x, Double_t *par) +{ + float func = sqrt((par[0]*par[0])/(x[0]*x[0])+par[1]*par[1]); + return func; +} + +void draw_req_DCA(double etamin, double etamax, double xmin=0., double xmax=0.) +{ + TF1 *PWGReq_DCA2D; + if (etamin >= -3.5 && etamax <= -2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax); + else if (etamin >= -2.5 && etamax <= -1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax); + else if (etamin >= -1.0 && etamax <= 1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((20./x)^2+5.0^2)",xmin,xmax); + else if (etamin >= 1.0 && etamax <= 2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax); + else if (etamin >= 2.5 && etamax <= 3.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax); + else return; + PWGReq_DCA2D->SetLineStyle(7); + PWGReq_DCA2D->SetLineWidth(3.0); + PWGReq_DCA2D->SetLineColor(kBlue); + PWGReq_DCA2D->Draw("same"); + + TLegend *l= new TLegend(0.70,0.75,0.90,0.80); + l->SetTextSize(0.03); + l->SetBorderSize(0); + l->AddEntry(PWGReq_DCA2D,"PWGReq","l"); + l->Draw("same"); +} From c46707edc686d08a002d01b671b797292ddad5ff Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 20 Jul 2024 18:07:19 -0400 Subject: [PATCH 11/19] tracking_performance: enable running over campaigns --- .gitlab-ci.yml | 3 + benchmarks/tracking_performances/README.md | 11 +++- benchmarks/tracking_performances/Snakefile | 62 ++++++++++++++------- benchmarks/tracking_performances/config.yml | 17 +++++- 4 files changed, 70 insertions(+), 23 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 03691293..00b6aedd 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -121,6 +121,9 @@ get_data: - mkdir "${DETECTOR_CONFIG}" - ln -s "${LOCAL_DATA_PATH}/sim_output" "${DETECTOR_CONFIG}/sim_output" - ln -s "../results" "${DETECTOR_CONFIG}/results" + # cache downloaded artifacts + - if [ -d /scratch ]; then mkdir -p /scratch/EPIC; ln -sf /scratch/EPIC ./EPIC; fi + - du -hs /scratch/EPIC - ls -lrtha retry: max: 2 diff --git a/benchmarks/tracking_performances/README.md b/benchmarks/tracking_performances/README.md index a5610e63..acf8844c 100644 --- a/benchmarks/tracking_performances/README.md +++ b/benchmarks/tracking_performances/README.md @@ -4,7 +4,16 @@ The scripts can be used to create the debug plots for the momentum resolutions. To run a full simulation-reconstruction-analysis chain do: ``` -snakemake --cores 1 tracking_performance +snakemake -c2 results/tracking_performances/local +``` +or, referencing it by the rule name: +``` +snakemake --cores 1 tracking_performance_local +``` + +To process an individual campaign do: +``` +snakemake -c2 results/tracking_performances/24.04.0 ``` It will produce the results with truth/realistic seeding. diff --git a/benchmarks/tracking_performances/Snakefile b/benchmarks/tracking_performances/Snakefile index 8bbdee90..aa3d888c 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -48,27 +48,26 @@ rule tracking_performance_at_momentum: input: script="benchmarks/tracking_performances/Tracking_Performances.C", # TODO pass as a file list? - sim=lambda wildcards: branch( - wildcards.CAMPAIGN == "local", - then=expand( + sim=lambda wildcards: + expand( "sim_output/tracking_performance/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", DETECTOR_CONFIG="epic_craterlake_tracking_only", ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", PHASE_SPACE=["3to50deg", "45to135deg", "130to177deg"], INDEX=range(1), - ), - otherwise=expand( + ) + if wildcards.CAMPAIGN == "local" else + expand( "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", DETECTOR_CONFIG="epic_craterlake", ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", PHASE_SPACE=["3to50deg", "45to135deg", "130to177deg"], INDEX=range(1), ), - ), output: - "{SEEDING}/pi-/mom/Performances_mom_{MOMENTUM}_mom_resol_{SEEDING_IGNORE}_{PARTICLE}.root", - "{SEEDING}/pi-/dca/Performances_dca_{MOMENTUM}_dca_resol_{SEEDING_IGNORE}_{PARTICLE}.root", - combined_root=temp("sim_{SEEDING}_{MOMENTUM}_{SEEDING_IGNORE}_{PARTICLE}.root"), + "{CAMPAIGN}/{SEEDING}/pi-/mom/Performances_mom_{MOMENTUM}_mom_resol_{SEEDING_IGNORE}_{PARTICLE}.root", + "{CAMPAIGN}/{SEEDING}/pi-/dca/Performances_dca_{MOMENTUM}_dca_resol_{SEEDING_IGNORE}_{PARTICLE}.root", + combined_root=temp("{CAMPAIGN}/sim_{SEEDING}_{MOMENTUM}_{SEEDING_IGNORE}_{PARTICLE}.root"), shell: """ if [[ "{wildcards.SEEDING}" == "truthseed" ]]; then @@ -77,7 +76,8 @@ elif [[ "{wildcards.SEEDING}" == "realseed" ]]; then SEEDING="Seeded" fi hadd {output.combined_root} {input.sim} -root -l -b -q {input.script}'("{output.combined_root}", "{wildcards.PARTICLE}", {wildcards.MOMENTUM}, 0.15, "'$SEEDING'")' +cd {wildcards.CAMPAIGN} +root -l -b -q ../{input.script}'("../{output.combined_root}", "{wildcards.PARTICLE}", {wildcards.MOMENTUM}, 0.15, "'$SEEDING'")' """ @@ -85,25 +85,26 @@ rule tracking_performance_summary_at_eta: input: expand( [ - "truthseed/pi-/mom/Performances_mom_{MOMENTUM:.1f}_mom_resol_truth_pi-.root", - "truthseed/pi-/dca/Performances_dca_{MOMENTUM:.1f}_dca_resol_truth_pi-.root", - "realseed/pi-/mom/Performances_mom_{MOMENTUM:.1f}_mom_resol_realseed_pi-.root", - "realseed/pi-/dca/Performances_dca_{MOMENTUM:.1f}_dca_resol_realseed_pi-.root", + "{{CAMPAIGN}}/truthseed/pi-/mom/Performances_mom_{MOMENTUM:.1f}_mom_resol_truth_pi-.root", + "{{CAMPAIGN}}/truthseed/pi-/dca/Performances_dca_{MOMENTUM:.1f}_dca_resol_truth_pi-.root", + "{{CAMPAIGN}}/realseed/pi-/mom/Performances_mom_{MOMENTUM:.1f}_mom_resol_realseed_pi-.root", + "{{CAMPAIGN}}/realseed/pi-/dca/Performances_dca_{MOMENTUM:.1f}_dca_resol_realseed_pi-.root", ], MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], ), script="benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C", output: expand( - "Debug_Plots/{SEEDING}/pi-/mom/{SEEDING}_mom_resol_mom{MOMENTUM:.1f}_{{ETA_MIN}}_eta_{{ETA_MAX}}.png", + "{{CAMPAIGN}}/Debug_Plots/{SEEDING}/pi-/mom/{SEEDING}_mom_resol_mom{MOMENTUM:.1f}_{{ETA_MIN}}_eta_{{ETA_MAX}}.png", SEEDING=["real", "truth"], MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], ), - "Final_Results/pi-/mom/mom_resol_{ETA_MIN}_eta_{ETA_MAX}.png", - "Final_Results/pi-/mom/mom_resol_{ETA_MIN}_eta_{ETA_MAX}.root", + "{CAMPAIGN}/Final_Results/pi-/mom/mom_resol_{ETA_MIN}_eta_{ETA_MAX}.png", + "{CAMPAIGN}/Final_Results/pi-/mom/mom_resol_{ETA_MIN}_eta_{ETA_MAX}.root", shell: """ -root -l -b -q {input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1.)' +cd {wildcards.CAMPAIGN} +root -l -b -q ../{input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1.)' """ @@ -113,15 +114,34 @@ rule tracking_performance: input: expand( [ - "Final_Results/pi-/mom/mom_resol_{ETA_BIN}.png", - "Final_Results/pi-/mom/mom_resol_{ETA_BIN}.root", + "{{CAMPAIGN}}/Final_Results/pi-/mom/mom_resol_{ETA_BIN}.png", + "{{CAMPAIGN}}/Final_Results/pi-/mom/mom_resol_{ETA_BIN}.root", ], ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(TRACKING_PERFORMANCE_ETA_BINS[:-1], TRACKING_PERFORMANCE_ETA_BINS[1:])], ), output: - directory("results/tracking_performances") + directory("results/tracking_performances/{CAMPAIGN}/") shell: """ mkdir {output} cp {input} {output} """ + + +rule tracking_performance_local: + input: + "results/tracking_performances/local", + + +rule tracking_performance_campaigns: + input: + expand( + "results/tracking_performances/{CAMPAIGN}", + CAMPAIGN=[ + "23.10.0", + "23.11.0", + "23.12.0", + "24.03.1", + "24.04.0", + ], + ) diff --git a/benchmarks/tracking_performances/config.yml b/benchmarks/tracking_performances/config.yml index 22290396..fadeec02 100644 --- a/benchmarks/tracking_performances/config.yml +++ b/benchmarks/tracking_performances/config.yml @@ -18,7 +18,7 @@ bench:tracking_performance: needs: - ["sim:tracking_performance"] script: - - snakemake --cores 1 tracking_performance + - snakemake --cores 1 tracking_performance_local collect_results:tracking_performance: extends: .det_benchmark @@ -27,3 +27,18 @@ collect_results:tracking_performance: - "bench:tracking_performance" script: - ls -lrht + +bench:tracking_performance_campaigns: + extends: .det_benchmark + stage: benchmarks + #when: manual + script: + - snakemake --cores 1 tracking_performance_campaigns + +collect_results:tracking_performance_campaigns: + extends: .det_benchmark + stage: collect + needs: + - "bench:tracking_performance_campaigns" + script: + - ls -lrht From 03a6adcda500609eeae128a9841e529138d4d6c9 Mon Sep 17 00:00:00 2001 From: Barak Schmookler Date: Sat, 27 Jul 2024 00:33:04 -0400 Subject: [PATCH 12/19] Benchmark for low-energy photons in the ZDC LYSO (#37) --- .gitlab-ci.yml | 1 + Snakefile | 2 +- benchmarks/zdc_lyso/README.md | 11 ++ benchmarks/zdc_lyso/Snakefile | 86 ++++++++++ benchmarks/zdc_lyso/analysis/analysis.py | 195 +++++++++++++++++++++++ benchmarks/zdc_lyso/config.yml | 17 ++ benchmarks/zdc_lyso/gen_particles.cxx | 126 +++++++++++++++ 7 files changed, 437 insertions(+), 1 deletion(-) create mode 100644 benchmarks/zdc_lyso/README.md create mode 100644 benchmarks/zdc_lyso/Snakefile create mode 100644 benchmarks/zdc_lyso/analysis/analysis.py create mode 100644 benchmarks/zdc_lyso/config.yml create mode 100644 benchmarks/zdc_lyso/gen_particles.cxx diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 03691293..3b7997b0 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -135,6 +135,7 @@ include: - local: 'benchmarks/barrel_ecal/config.yml' - local: 'benchmarks/barrel_hcal/config.yml' - local: 'benchmarks/zdc/config.yml' + - local: 'benchmarks/zdc_lyso/config.yml' - local: 'benchmarks/material_maps/config.yml' - local: 'benchmarks/material_scan/config.yml' - local: 'benchmarks/pid/config.yml' diff --git a/Snakefile b/Snakefile index e7a54350..26cfc5d0 100644 --- a/Snakefile +++ b/Snakefile @@ -3,7 +3,7 @@ include: "benchmarks/barrel_ecal/Snakefile" include: "benchmarks/ecal_gaps/Snakefile" include: "benchmarks/material_scan/Snakefile" include: "benchmarks/tracking_performances/Snakefile" - +include: "benchmarks/zdc_lyso/Snakefile" rule fetch_epic: output: diff --git a/benchmarks/zdc_lyso/README.md b/benchmarks/zdc_lyso/README.md new file mode 100644 index 00000000..3d53843e --- /dev/null +++ b/benchmarks/zdc_lyso/README.md @@ -0,0 +1,11 @@ +Detector Benchmark for LYSO ZDC +=============================== + +## Overview +This benchmark generates events with single low-energy (5 MeV - 1 GeV) photons. The photons are generated with angles of 0-5.24 mRad with respect to the proton/ion beam direction. The benchmark creates performance plots of the LYSO ZDC acceptance and energy reconstruction resolution. + +## Contacts +[JiaJun Huang](jhuan328@ucr.edu) +[Barak Schmookler](baraks@ucr.edu) + + diff --git a/benchmarks/zdc_lyso/Snakefile b/benchmarks/zdc_lyso/Snakefile new file mode 100644 index 00000000..5df026df --- /dev/null +++ b/benchmarks/zdc_lyso/Snakefile @@ -0,0 +1,86 @@ +import os + + +rule ecal_lyso_sim_hepmc: + input: + script = "benchmarks/zdc_lyso/gen_particles.cxx", + output: + hepmcfile="data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + log: + "data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc.log", + params: + num_events=1000, + shell: + """ +root -l -b -q '{input.script}({params.num_events}, "{output.hepmcfile}", "{wildcards.PARTICLE}", {wildcards.THETA_MIN}, {wildcards.THETA_MAX}, 0, 360, {wildcards.BEAM_ENERGY})' +""" + + +rule ecal_lyso_sim: + input: + hepmcfile="data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", + output: + "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", + log: + "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root.log", + params: + num_events=1000, + shell: + """ +npsim \ + --runType batch \ + -v WARNING \ + --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ + --numberOfEvents {params.num_events} \ + --inputFiles {input.hepmcfile} \ + --outputFile {output} +""" + + +rule ecal_lyso_reco: + input: + "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", + output: + "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", + log: + "sim_output/zdc_lyso/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root.log", + shell: + """ +eicrecon -Ppodio:output_collections=HcalFarForwardZDCRawHits,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,EcalFarForwardZDCRawHits,EcalFarForwardZDCRecHits,EcalFarForwardZDCClusters,MCParticles {input} +mv podio_output.root {output} +""" + + +rule zdc_analysis: + input: + expand("sim_output/zdc_lyso/{{DETECTOR_CONFIG}}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", + PARTICLE=["gamma"], + BEAM_ENERGY=["0.005", "0.01", "0.05", "0.1", "0.5", "1.0"], + THETA_MIN=["0"], + THETA_MAX=["0.3"]), + script="benchmarks/zdc_lyso/analysis/analysis.py", + output: + "results/{DETECTOR_CONFIG}/zdc_lyso/plots.pdf", + shell: + """ +python {input.script} +""" + + +# Examples of invocation +rule create_all_hepmc: + input: + expand("data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + PARTICLE=["gamma"], + BEAM_ENERGY=["0.005", "0.01", "0.05", "0.1", "0.5", "1.0"], + THETA_MIN=["0"], + THETA_MAX=["0.3"]) + + +rule run_all_locally: + input: + "results/" + os.environ["DETECTOR_CONFIG"] + "/zdc_lyso/plots.pdf", + message: + "See output in {input[0]}" + diff --git a/benchmarks/zdc_lyso/analysis/analysis.py b/benchmarks/zdc_lyso/analysis/analysis.py new file mode 100644 index 00000000..04bd233f --- /dev/null +++ b/benchmarks/zdc_lyso/analysis/analysis.py @@ -0,0 +1,195 @@ +import numpy as np +import matplotlib.pyplot as plt +import mplhep as hep +import uproot +import pandas as pd +import os +from scipy.optimize import curve_fit +from matplotlib.backends.backend_pdf import PdfPages + +plt.figure() +hep.set_style(hep.style.CMS) +hep.set_style("CMS") + +def gaussian(x, amp, mean, sigma): + return amp * np.exp( -(x - mean)**2 / (2*sigma**2) ) + +def rotateY(xdata, zdata, angle): + s = np.sin(angle) + c = np.cos(angle) + rotatedz = c*zdata - s*xdata + rotatedx = s*zdata + c*xdata + return rotatedx, rotatedz + +Energy = [0.005, 0.01, 0.05, 0.1, 0.5, 1.0] +q0 = [5, 10, 40, 90, 400, 700] +q1 = [0.5, 0.5, 0.9, 5, 10, 20] + +df = pd.DataFrame({}) +for eng in Energy: + tree = uproot.open(f'sim_output/zdc_lyso/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.eicrecon.tree.edm4eic.root')['events'] + ecal_reco_energy = tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.energy'].array() + #hcal_reco_energy = tree['HcalFarForwardZDCClusters/HcalFarForwardZDCClusters.energy'].array() + + tree = uproot.open(f'sim_output/zdc_lyso/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.edm4hep.root')['events'] + ecal_sim_energy = tree['EcalFarForwardZDCHits/EcalFarForwardZDCHits.energy'].array() + hcal_sim_energy = tree['HcalFarForwardZDCHits/HcalFarForwardZDCHits.energy'].array() + + par_x = tree['MCParticles/MCParticles.momentum.x'].array()[:,2] + par_y = tree['MCParticles/MCParticles.momentum.y'].array()[:,2] + par_z = tree['MCParticles/MCParticles.momentum.z'].array()[:,2] + + eng = int(eng*1000) + + ecal_reco_energy = pd.DataFrame({f'ecal_reco_energy_{eng}': np.array(ecal_reco_energy.tolist(), dtype=object)}) + #hcal_reco_energy = pd.DataFrame({f'hcal_reco_energy_{eng}': np.array(hcal_reco_energy.tolist(), dtype=object)}) + ecal_sim_energy = pd.DataFrame({f'ecal_sim_energy_{eng}': np.array(ecal_sim_energy.tolist(), dtype=object)}) + hcal_sim_energy = pd.DataFrame({f'hcal_sim_energy_{eng}': np.array(hcal_sim_energy.tolist(), dtype=object)}) + + par_x = pd.DataFrame({f'par_x_{eng}': np.array(par_x.tolist(), dtype=object)}) + par_y = pd.DataFrame({f'par_y_{eng}': np.array(par_y.tolist(), dtype=object)}) + par_z = pd.DataFrame({f'par_z_{eng}': np.array(par_z.tolist(), dtype=object)}) + + df = pd.concat([df,ecal_reco_energy,ecal_sim_energy,hcal_sim_energy,par_x,par_y,par_z],axis=1) + +mu = [] +sigma = [] +resolution = [] + +fig1, ax = plt.subplots(3,2,figsize=(20,10)) +#fig1.suptitle('ZDC ECal Cluster Energy Reconstruction') + +plt.tight_layout() +for i in range(6): + plt.sca(ax[i%3,i//3]) + eng = int(Energy[i]*1000) + plt.title(f'Gamma Energy: {eng} MeV') + temp = np.array([sum(item) if (item != 0) else 0 for item in df[f'ecal_reco_energy_{eng}']]) + hist, x = np.histogram(np.array(temp)*1000,bins=30) + x = x[1:]/2 + x[:-1]/2 + plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o') + temp = np.array([item[0] for item in df[f'ecal_reco_energy_{eng}'] if item]) + hist, x = np.histogram(np.array(temp)*1000,bins=30) + x = x[1:]/2 + x[:-1]/2 + coeff, covar = curve_fit(gaussian,x,hist,p0=(200,q0[i],q1[i]),maxfev = 80000) + plt.plot(np.linspace(coeff[1]-5*coeff[2],coeff[1]+5*coeff[2],50),gaussian(np.linspace(coeff[1]-5*coeff[2],coeff[1]+5*coeff[2],50),*coeff) + ,label = f'$\mu$ = {coeff[1]:.3f} $\pm$ {covar[1][1]:.3f}\n$\sigma$ = {np.abs(coeff[2]):.3f} $\pm$ {covar[2][2]:.3f}\nResolution = {np.abs(coeff[2])*100/coeff[1]:.2f}%') + + plt.xlabel('Energy (MeV)') + plt.legend() + mu.append(coeff[1]) + sigma.append(coeff[2]) + resolution.append(np.abs(coeff[2])*100/coeff[1]) +#plt.savefig('results/Energy_reconstruction_cluster.pdf') +#plt.show() + +fig2, (ax1,ax2) = plt.subplots(2,1,figsize=(15,10),sharex=True) + +plt.tight_layout() +# Plot data on primary axis +ax1.scatter(np.array(Energy)*1000, mu, c='b') +ax1.plot([4.5,1000],[4.5,1000],c='black',label='x=y') +ax1.set_ylabel('Reconstructed Energy (MeV)') +ax1.set_yscale('log') +ax1.legend() +ax1.set_title('ECal Craterlake Cluster Energy Reconstruction') + +ax2.plot(np.array(Energy)*1000, resolution, c='r') +ax2.scatter(np.array(Energy)*1000, resolution, c='r') +ax2.set_ylabel('Resolution (%)') +ax2.set_xlabel('Gamma Energy (MeV)') +ax2.set_xscale('log') +#plt.savefig('results/Energy_resolution.pdf') +#plt.show() + +htower = [] +herr = [] +hmean = [] +hhits = [] +hhits_cut = [] +emean = [] +ehits = [] +etower = [] +eerr = [] +ehits_cut = [] + +fig3, ax = plt.subplots(2,3,figsize=(20,10)) +fig3.suptitle('ZDC Simulation Energy Reconstruction') +for i in range(6): + plt.sca(ax[i//3,i%3]) + eng = int(Energy[i]*1000) + + x = df[f'par_x_{eng}'].astype(float).to_numpy() + y = df[f'par_y_{eng}'].astype(float).to_numpy() + z = df[f'par_z_{eng}'].astype(float).to_numpy() + x, z = rotateY(x,z, 0.025) + theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 + condition = theta <= 3.5 + + plt.title(f'Gamma Energy: {eng} MeV') + energy1 = np.array([sum(item) if (item != 0) else 0 for item in df[f'hcal_sim_energy_{eng}']])#df.eval(f'hcal_sim_energy_{eng}').apply(lambda row: sum(row)) + hist, x = np.histogram(energy1*1000,bins=np.logspace(0,3,200)) + x = x[1:]/2 + x[:-1]/2 + plt.plot(x,hist,marker='o',label="HCal") + hhits.append(len(energy1[energy1!=0])) + condition1 = energy1!=0 + hhits_cut.append(len(energy1[condition & condition1])/len(condition[condition==True])) + energy = np.array([sum(item) if (item != 0) else 0 for item in df[f'ecal_sim_energy_{eng}']])#df.eval(f'ecal_sim_energy_{eng}').apply(lambda row: sum(row)) + hist, x = np.histogram(energy*1000,bins=np.logspace(0,3,200)) + x = x[1:]/2 + x[:-1]/2 + plt.plot(x,hist,marker='o',label="ECal") + emean.append(sum(energy[energy!=0])*1000/len(energy[energy!=0])) + hmean.append(sum(energy1[energy!=0])*1000/len(energy[energy!=0])) + condition1 = energy!=0 + ehits_cut.append(len(energy[condition & condition1])/len(condition[condition==True])) + ehits.append(len(energy[energy!=0])) + plt.legend() + plt.xscale('log') + plt.xlim(1e0,1e3) + +plt.xlabel('Energy (MeV)') +#plt.savefig('results/Energy_deposition.pdf') +#plt.show() + +fig4, ax = plt.subplots(2,1,sharex=True,gridspec_kw={'height_ratios': [2,1]}) +plt.sca(ax[0]) +plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619+np.array(emean),label='HCal+ECal',fmt='-o') +plt.errorbar(np.array(Energy)*1000,emean,label='ECal',fmt='-o') +plt.legend() +plt.yscale('log') +plt.xscale('log') +plt.ylabel('Simulation Energy (MeV)') +plt.sca(ax[1]) +plt.errorbar(np.array(Energy)*1000,(1 - np.array(emean)/(np.array(hmean)*47.619+np.array(emean)))*100,label='Total/ECal',fmt='-o') +plt.legend() +plt.ylabel('Fraction of energy\n deposited in Hcal (%)') +plt.xlabel('Truth Energy (MeV)') +#plt.savefig('results/Energy_ratio_and_Leakage.pdf') +plt.tight_layout() +#plt.show() + +fig5 = plt.figure() +plt.errorbar(np.array(Energy)*1000,np.array(hhits)/1000*100,label='HCal Hits',fmt='-o') +plt.errorbar(np.array(Energy)*1000,np.array(ehits)/1000*100,label='ECal Hits',fmt='-o') +#plt.errorbar(np.array(Energy)*1000,np.array(hhits)/np.array(ehits)*100,label='HCal / ECal',fmt='-o',c='b') + +plt.errorbar(np.array(Energy)*1000,np.array(hhits_cut)*100,label='HCal Hits with 3.5 mRad cut',fmt='-^') +plt.errorbar(np.array(Energy)*1000,np.array(ehits_cut)*100,label='ECal Hits with 3.5 mRad cut',fmt='-^') +#plt.errorbar(np.array(Energy)*1000,np.array(hhits_cut)/np.array(ehits_cut)*100,label='HCal / ECal with 3.5 mRad cut',fmt='-^',c='b') +### 3mrad cuts + +plt.legend() +plt.xlabel('Simulation Truth Gamma Energy (MeV)') +plt.ylabel('Fraction of Events with non-zero energy (%)') +#plt.savefig('results/Hits.pdf') +plt.xscale('log') +#plt.show() + +#pdfs = ['results/Energy_reconstruction_cluster.pdf','results/Energy_resolution.pdf','results/Energy_deposition.pdf','results/Energy_ratio_and_Leakage.pdf','results/Hits.pdf'] +with PdfPages(f'results/{os.environ["DETECTOR_CONFIG"]}/zdc_lyso/plots.pdf') as pdf: + pdf.savefig(fig1) + pdf.savefig(fig2) + pdf.savefig(fig3) + pdf.savefig(fig4) + pdf.savefig(fig5) + diff --git a/benchmarks/zdc_lyso/config.yml b/benchmarks/zdc_lyso/config.yml new file mode 100644 index 00000000..7810e0cf --- /dev/null +++ b/benchmarks/zdc_lyso/config.yml @@ -0,0 +1,17 @@ +sim:zdc_lyso: + extends: .det_benchmark + stage: simulate + script: + - snakemake --cores 1 run_all_locally + retry: + max: 2 + when: + - runner_system_failure + +collect_results:zdc_lyso: + extends: .det_benchmark + stage: collect + needs: + - "sim:zdc_lyso" + script: + - ls -lrht diff --git a/benchmarks/zdc_lyso/gen_particles.cxx b/benchmarks/zdc_lyso/gen_particles.cxx new file mode 100644 index 00000000..b84b1170 --- /dev/null +++ b/benchmarks/zdc_lyso/gen_particles.cxx @@ -0,0 +1,126 @@ +#include "HepMC3/GenEvent.h" +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/WriterAscii.h" +#include "HepMC3/Print.h" + +#include "TRandom3.h" +#include "TVector3.h" + +#include +#include +#include +#include +#include +#include +#include + +using namespace HepMC3; + +// Generate single electron as input to the Insert simulation. +// -- +// We generate events with a constant polar angle with respect to +// the proton direction and then rotate so that the events are given +// in normal lab coordinate system +// -- +void gen_particles( + int n_events = 10, + const char* out_fname = "gen_particles.hepmc", + TString particle_name = "e-", + double th_min = 3., // Minimum polar angle, in degrees + double th_max = 3., // Maximum polar angle, in degrees + double phi_min = 0., // Minimum azimuthal angle, in degrees + double phi_max = 360., // Maximum azimuthal angle, in degrees + double p = 10., // Momentum in GeV/c + int dist = 0 //Momentum distribution: 0=fixed, 1=uniform, 2=Gaussian + ) +{ + WriterAscii hepmc_output(out_fname); + int events_parsed = 0; + GenEvent evt(Units::GEV, Units::MM); + + // Random number generator + TRandom3 *r1 = new TRandom3(0); //Use time as random seed + + // Getting generated particle information + TDatabasePDG *pdg = new TDatabasePDG(); + TParticlePDG *particle = pdg->GetParticle(particle_name); + const double mass = particle->Mass(); + const int pdgID = particle->PdgCode(); + + for (events_parsed = 0; events_parsed < n_events; events_parsed++) { + + //Set the event number + evt.set_event_number(events_parsed); + + // FourVector(px,py,pz,e,pdgid,status) + // type 4 is beam + // pdgid 11 - electron + // pdgid 111 - pi0 + // pdgid 2212 - proton + GenParticlePtr p1 = + std::make_shared(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); + GenParticlePtr p2 = std::make_shared( + FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4); + + // Define momentum with respect to proton direction + double phi = r1->Uniform(phi_min*TMath::DegToRad(),phi_max*TMath::DegToRad()); + double th = r1->Uniform(th_min*TMath::DegToRad(),th_max*TMath::DegToRad()); + + //Total momentum distribution + double pevent = -1; + if(dist==0){ //fixed + pevent = p; + } + else if(dist==1){ //Uniform: +-50% variation + pevent = p*(1. + r1->Uniform(-0.5,0.5) ); + } + else if(dist==2){ //Gaussian: Sigma = 0.1*mean + while(pevent<0) //Avoid negative values + pevent = r1->Gaus(p,0.1*p); + } + + double px = pevent * std::cos(phi) * std::sin(th); + double py = pevent * std::sin(phi) * std::sin(th); + double pz = pevent * std::cos(th); + TVector3 pvec(px,py,pz); + + //Rotate to lab coordinate system + double cross_angle = -25./1000.; //in Rad + TVector3 pbeam_dir(sin(cross_angle),0,cos(cross_angle)); //proton beam direction + pvec.RotateY(-pbeam_dir.Theta()); // Theta is returned positive, beam in negative X + // type 1 is final state + // pdgid 11 - electron 0.510 MeV/c^2 + GenParticlePtr p3 = std::make_shared( + FourVector( + pvec.X(), pvec.Y(), pvec.Z(), + sqrt(pevent*pevent + (mass * mass))), + pdgID, 1); + + //If wanted, set non-zero vertex + double vx = 0.; + double vy = 0.; + double vz = 0.; + double vt = 0.; + + GenVertexPtr v1 = std::make_shared(); + evt.shift_position_by(FourVector(vx, vy, vz, vt)); + v1->add_particle_in(p1); + v1->add_particle_in(p2); + + v1->add_particle_out(p3); + evt.add_vertex(v1); + + if (events_parsed == 0) { + std::cout << "First event: " << std::endl; + Print::listing(evt); + } + + hepmc_output.write_event(evt); + if (events_parsed % 100 == 0) { + std::cout << "Event: " << events_parsed << std::endl; + } + evt.clear(); + } + hepmc_output.close(); + std::cout << "Events parsed and written: " << events_parsed << std::endl; +} From 0a6f4d17076802b4682616cd4a0f41265f1a0250 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Sun, 28 Jul 2024 05:58:14 +0200 Subject: [PATCH 13/19] Update doCompare_truth_real_widebins_mom.C --- .../tracking_performances/doCompare_truth_real_widebins_mom.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C index 283dffce..4d71e104 100644 --- a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C @@ -151,7 +151,7 @@ void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1 mgMom->Add(gr2); c_mom->cd(); mgMom->GetXaxis()->SetRangeUser(0.40,20.2); - mgMom->GetYaxis()->SetRangeUser(0.,10.0); + mgMom->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr2->GetN(),gr2->GetY())); // 50% more of the maximum value on yaxis mgMom->Draw("AP"); lmom->AddEntry(gr1,"Truth Seeding"); lmom->AddEntry(gr2,"Realistic Seeding"); From 45d83752f89dd59ad93c4eed264fc9b5caaf22ee Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sun, 28 Jul 2024 00:01:00 -0400 Subject: [PATCH 14/19] ecal_gaps: use DETECTOR_CONFIG from the wildcard in eicrecon (#38) --- benchmarks/ecal_gaps/Snakefile | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/benchmarks/ecal_gaps/Snakefile b/benchmarks/ecal_gaps/Snakefile index 7046eb60..91ffa7c2 100644 --- a/benchmarks/ecal_gaps/Snakefile +++ b/benchmarks/ecal_gaps/Snakefile @@ -41,7 +41,8 @@ rule ecal_gaps_recon: wildcard_constraints: INDEX="\d{4}", shell: """ -eicrecon {input} -Ppodio:output_file={output} \ +env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \ + eicrecon {input} -Ppodio:output_file={output} \ -Ppodio:output_include_collections=EcalEndcapNRecHits,EcalBarrelScFiRecHits,EcalBarrelImagingRecHits,EcalEndcapPRecHits,MCParticles """ From d73fa659375b0e16ff21cbd821ca1dea328abcc6 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 29 Jul 2024 14:01:56 -0400 Subject: [PATCH 15/19] correctly set campaign version for the plot legend --- benchmarks/tracking_performances/Snakefile | 22 +++++++++++++------ .../doCompare_truth_real_widebins_mom.C | 7 +++--- 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/benchmarks/tracking_performances/Snakefile b/benchmarks/tracking_performances/Snakefile index 860e6bc4..f9f04237 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -103,14 +103,22 @@ rule tracking_performance_summary_at_eta: "{CAMPAIGN}/Final_Results/pi-/mom/mom_resol_{ETA_MIN}_eta_{ETA_MAX}.root", shell: """ -set +e -EPIC_VERSION="$(echo $DETECTOR_PATH | sed -n -e 's/.*epic-\([^-]\+\).*/\\1/p')" -EICRECON_VERSION="$(eicrecon -v | sed -n -e 's/.*\(v[0-9.]\+\).*/\\1/p')" -set -e -echo "ePIC version: $EPIC_VERSION" -echo "EICrecon version: $EICRECON_VERSION" +if [[ "{wildcards.CAMPAIGN}" == "local" ]]; then + set +e + EPIC_VERSION="${{DETECTOR_VERSION:-}}" + EICRECON_VERSION="$(eicrecon -v | sed -n -e 's/.*\(v[0-9\.]\+\).*/\\1/p')" + # Legacy detection + : ${{EPIC_VERSION:="$(echo $DETECTOR_PATH | sed -n -e 's/.*epic-\([^-/]\+\).*/\\1/p')"}} + set -e + + echo "ePIC version: $EPIC_VERSION" + echo "EICrecon version: $EICRECON_VERSION" + EXTRA_LEGEND="ePIC $EPIC_VERSION / EICrecon $EICRECON_VERSION" +else + EXTRA_LEGEND="ePIC Simulation {wildcards.CAMPAIGN}" +fi cd {wildcards.CAMPAIGN} -root -l -b -q ../{input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'$EPIC_VERSION'", "'$EICRECON_VERSION'")' +root -l -b -q ../{input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'"$EXTRA_LEGEND"'")' """ diff --git a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C index 6416d268..153f5360 100644 --- a/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C +++ b/benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.C @@ -10,7 +10,7 @@ #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, TString epic ="", TString eicrecon = "") // name = p, pt for getting p or pt dependence fitted results +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, TString extra_legend = "") // name = p, pt for getting p or pt dependence fitted results { //=== style of the plot========= @@ -52,7 +52,8 @@ void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1 lmom = new TLegend(0.65,0.80,0.90,0.93); lmom->SetTextSize(0.03); lmom->SetBorderSize(0); - lmom->SetHeader(Form("%s ePIC(%s/%s): %1.1f < #eta < %1.1f",symbolname.Data(),epic.Data(),eicrecon.Data(),etamin,etamax),"C"); + lmom->SetHeader(extra_legend.Data(), "C"); + lmom->AddEntry((TObject*)0, Form("%s, %1.1f < #eta < %1.1f", symbolname.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); @@ -165,7 +166,7 @@ void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1 c_mom->SaveAs(Form("Final_Results/%s/mom/mom_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax)); // Write the numbers in output file for comparisons - outfile<<"ePIC"<GetN(); ++i){ double x,ytrue, yreal; From e815da28f7a47bf86215e8a3e3feeacca9d92bd8 Mon Sep 17 00:00:00 2001 From: Shyam Kumar Date: Tue, 30 Jul 2024 07:55:02 +0200 Subject: [PATCH 16/19] Update Tracking_Performances.C --- benchmarks/tracking_performances/Tracking_Performances.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/tracking_performances/Tracking_Performances.C b/benchmarks/tracking_performances/Tracking_Performances.C index bfb751d6..77fc0202 100644 --- a/benchmarks/tracking_performances/Tracking_Performances.C +++ b/benchmarks/tracking_performances/Tracking_Performances.C @@ -82,7 +82,7 @@ void Tracking_Performances(TString filename="tracking_output",TString particle=" 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 pmc = (1./charge[j])*sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec Double_t prec = 1./qoverp[j]; Double_t pzrec = prec*TMath::Cos(theta[j]); Double_t pt_rec = sqrt(prec*prec-pzrec*pzrec); From e3b0d5f2d3170b6898c363d382e859e87bb8df07 Mon Sep 17 00:00:00 2001 From: Jiajun Huang Date: Wed, 31 Jul 2024 15:55:18 -0700 Subject: [PATCH 17/19] Modified zdc lyso analysis code. --- benchmarks/zdc_lyso/analysis/analysis.py | 164 +++++++++++++++++------ 1 file changed, 125 insertions(+), 39 deletions(-) diff --git a/benchmarks/zdc_lyso/analysis/analysis.py b/benchmarks/zdc_lyso/analysis/analysis.py index 04bd233f..eace6112 100644 --- a/benchmarks/zdc_lyso/analysis/analysis.py +++ b/benchmarks/zdc_lyso/analysis/analysis.py @@ -3,9 +3,10 @@ import mplhep as hep import uproot import pandas as pd -import os from scipy.optimize import curve_fit from matplotlib.backends.backend_pdf import PdfPages +import os +import awkward as ak plt.figure() hep.set_style(hep.style.CMS) @@ -22,18 +23,21 @@ def rotateY(xdata, zdata, angle): return rotatedx, rotatedz Energy = [0.005, 0.01, 0.05, 0.1, 0.5, 1.0] -q0 = [5, 10, 40, 90, 400, 700] -q1 = [0.5, 0.5, 0.9, 5, 10, 20] + df = pd.DataFrame({}) for eng in Energy: tree = uproot.open(f'sim_output/zdc_lyso/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.eicrecon.tree.edm4eic.root')['events'] - ecal_reco_energy = tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.energy'].array() - #hcal_reco_energy = tree['HcalFarForwardZDCClusters/HcalFarForwardZDCClusters.energy'].array() + ecal_reco_energy = list(map(sum, tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.energy'].array())) + hcal_reco_energy = list(map(sum, tree['HcalFarForwardZDCClusters/HcalFarForwardZDCClusters.energy'].array())) + ecal_rec_energy = list(map(sum, tree['EcalFarForwardZDCRecHits/EcalFarForwardZDCRecHits.energy'].array())) + hcal_rec_energy = list(map(sum, tree['HcalFarForwardZDCRecHits/HcalFarForwardZDCRecHits.energy'].array())) + ecal_reco_clusters = [len(row) if len(row)>=1 else 0 for row in tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.nhits'].array()] + ecal_reco_nhits = [row[0] if len(row)>=1 else 0 for row in tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.nhits'].array()] tree = uproot.open(f'sim_output/zdc_lyso/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.edm4hep.root')['events'] - ecal_sim_energy = tree['EcalFarForwardZDCHits/EcalFarForwardZDCHits.energy'].array() - hcal_sim_energy = tree['HcalFarForwardZDCHits/HcalFarForwardZDCHits.energy'].array() + ecal_sim_energy = list(map(sum, tree['EcalFarForwardZDCHits/EcalFarForwardZDCHits.energy'].array())) + hcal_sim_energy = list(map(sum, tree['HcalFarForwardZDCHits/HcalFarForwardZDCHits.energy'].array())) par_x = tree['MCParticles/MCParticles.momentum.x'].array()[:,2] par_y = tree['MCParticles/MCParticles.momentum.y'].array()[:,2] @@ -41,66 +45,103 @@ def rotateY(xdata, zdata, angle): eng = int(eng*1000) - ecal_reco_energy = pd.DataFrame({f'ecal_reco_energy_{eng}': np.array(ecal_reco_energy.tolist(), dtype=object)}) - #hcal_reco_energy = pd.DataFrame({f'hcal_reco_energy_{eng}': np.array(hcal_reco_energy.tolist(), dtype=object)}) - ecal_sim_energy = pd.DataFrame({f'ecal_sim_energy_{eng}': np.array(ecal_sim_energy.tolist(), dtype=object)}) - hcal_sim_energy = pd.DataFrame({f'hcal_sim_energy_{eng}': np.array(hcal_sim_energy.tolist(), dtype=object)}) - + ecal_reco_energy = pd.DataFrame({f'ecal_reco_energy_{eng}': np.array(ecal_reco_energy, dtype=object)}) + hcal_reco_energy = pd.DataFrame({f'hcal_reco_energy_{eng}': np.array(hcal_reco_energy, dtype=object)}) + ecal_rec_energy = pd.DataFrame({f'ecal_rec_energy_{eng}': np.array(ecal_rec_energy, dtype=object)}) + hcal_rec_energy = pd.DataFrame({f'hcal_rec_energy_{eng}': np.array(hcal_rec_energy, dtype=object)}) + ecal_sim_energy = pd.DataFrame({f'ecal_sim_energy_{eng}': np.array(ecal_sim_energy, dtype=object)}) + hcal_sim_energy = pd.DataFrame({f'hcal_sim_energy_{eng}': np.array(hcal_sim_energy, dtype=object)}) + ecal_reco_nhits = pd.DataFrame({f'ecal_reco_nhits_{eng}': np.array(ecal_reco_nhits, dtype=object)}) + ecal_reco_clusters = pd.DataFrame({f'ecal_reco_clusters_{eng}': np.array(ecal_reco_clusters, dtype=object)}) par_x = pd.DataFrame({f'par_x_{eng}': np.array(par_x.tolist(), dtype=object)}) par_y = pd.DataFrame({f'par_y_{eng}': np.array(par_y.tolist(), dtype=object)}) par_z = pd.DataFrame({f'par_z_{eng}': np.array(par_z.tolist(), dtype=object)}) - df = pd.concat([df,ecal_reco_energy,ecal_sim_energy,hcal_sim_energy,par_x,par_y,par_z],axis=1) + + df = pd.concat([df,ecal_reco_energy,ecal_rec_energy,ecal_sim_energy,hcal_reco_energy,hcal_rec_energy,hcal_sim_energy,ecal_reco_clusters,ecal_reco_nhits,par_x,par_y,par_z],axis=1) + mu = [] sigma = [] -resolution = [] - fig1, ax = plt.subplots(3,2,figsize=(20,10)) -#fig1.suptitle('ZDC ECal Cluster Energy Reconstruction') +fig1.suptitle('ZDC ECal Cluster Energy Reconstruction') plt.tight_layout() for i in range(6): + x = df[f'par_x_{eng}'].astype(float).to_numpy() + y = df[f'par_y_{eng}'].astype(float).to_numpy() + z = df[f'par_z_{eng}'].astype(float).to_numpy() + x, z = rotateY(x,z, 0.025) + theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 + condition = theta <= 3.5 + plt.sca(ax[i%3,i//3]) eng = int(Energy[i]*1000) plt.title(f'Gamma Energy: {eng} MeV') - temp = np.array([sum(item) if (item != 0) else 0 for item in df[f'ecal_reco_energy_{eng}']]) - hist, x = np.histogram(np.array(temp)*1000,bins=30) + temp = np.array(df[f'ecal_reco_energy_{eng}'].astype(float).to_numpy()[condition])*1000 + hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) x = x[1:]/2 + x[:-1]/2 - plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o') - temp = np.array([item[0] for item in df[f'ecal_reco_energy_{eng}'] if item]) - hist, x = np.histogram(np.array(temp)*1000,bins=30) + plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Cluster') + coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0]))) + #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) + mu.append(coeff[1]) + sigma.append(coeff[2]) + + temp = np.array(df[f'ecal_rec_energy_{eng}'].astype(float).to_numpy()[condition])*1000 + hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) x = x[1:]/2 + x[:-1]/2 - coeff, covar = curve_fit(gaussian,x,hist,p0=(200,q0[i],q1[i]),maxfev = 80000) - plt.plot(np.linspace(coeff[1]-5*coeff[2],coeff[1]+5*coeff[2],50),gaussian(np.linspace(coeff[1]-5*coeff[2],coeff[1]+5*coeff[2],50),*coeff) - ,label = f'$\mu$ = {coeff[1]:.3f} $\pm$ {covar[1][1]:.3f}\n$\sigma$ = {np.abs(coeff[2]):.3f} $\pm$ {covar[2][2]:.3f}\nResolution = {np.abs(coeff[2])*100/coeff[1]:.2f}%') + plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Digitization') + coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0]))) + #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) + mu.append(coeff[1]) + sigma.append(coeff[2]) - plt.xlabel('Energy (MeV)') - plt.legend() + temp = np.array(df[f'ecal_sim_energy_{eng}'].astype(float).to_numpy()[condition])*1000 + hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) + x = x[1:]/2 + x[:-1]/2 + plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Simulation') + coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0]))) + #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) mu.append(coeff[1]) sigma.append(coeff[2]) - resolution.append(np.abs(coeff[2])*100/coeff[1]) + + plt.xlabel('Energy (MeV)') + plt.legend() + #plt.savefig('results/Energy_reconstruction_cluster.pdf') -#plt.show() + +mu = np.array(mu) +sigma = np.array(sigma) + +plt.show() fig2, (ax1,ax2) = plt.subplots(2,1,figsize=(15,10),sharex=True) plt.tight_layout() # Plot data on primary axis -ax1.scatter(np.array(Energy)*1000, mu, c='b') +ax1.scatter(np.array(Energy)*1000, mu[::3], label='cluster') +ax1.scatter(np.array(Energy)*1000, mu[1::3], label='digitization') +ax1.scatter(np.array(Energy)*1000, mu[2::3], label='simulation') + ax1.plot([4.5,1000],[4.5,1000],c='black',label='x=y') ax1.set_ylabel('Reconstructed Energy (MeV)') ax1.set_yscale('log') ax1.legend() ax1.set_title('ECal Craterlake Cluster Energy Reconstruction') -ax2.plot(np.array(Energy)*1000, resolution, c='r') -ax2.scatter(np.array(Energy)*1000, resolution, c='r') +ax2.errorbar(np.array(Energy)*1000, abs(sigma[::3]/mu[::3])*100, fmt='-o', label='cluster') +ax2.errorbar(np.array(Energy)*1000, abs(sigma[1::3]/mu[1::3])*100, fmt='-o', label='digitization') +ax2.errorbar(np.array(Energy)*1000, abs(sigma[2::3]/mu[2::3])*100, fmt='-o', label='simulation') + ax2.set_ylabel('Resolution (%)') ax2.set_xlabel('Gamma Energy (MeV)') ax2.set_xscale('log') +ax2.legend() + #plt.savefig('results/Energy_resolution.pdf') -#plt.show() + +plt.show() + htower = [] herr = [] @@ -127,14 +168,14 @@ def rotateY(xdata, zdata, angle): condition = theta <= 3.5 plt.title(f'Gamma Energy: {eng} MeV') - energy1 = np.array([sum(item) if (item != 0) else 0 for item in df[f'hcal_sim_energy_{eng}']])#df.eval(f'hcal_sim_energy_{eng}').apply(lambda row: sum(row)) + energy1 = df[f'hcal_sim_energy_{eng}'].astype(float).to_numpy()#df.eval(f'hcal_sim_energy_{eng}').apply(lambda row: sum(row)) hist, x = np.histogram(energy1*1000,bins=np.logspace(0,3,200)) x = x[1:]/2 + x[:-1]/2 plt.plot(x,hist,marker='o',label="HCal") hhits.append(len(energy1[energy1!=0])) condition1 = energy1!=0 hhits_cut.append(len(energy1[condition & condition1])/len(condition[condition==True])) - energy = np.array([sum(item) if (item != 0) else 0 for item in df[f'ecal_sim_energy_{eng}']])#df.eval(f'ecal_sim_energy_{eng}').apply(lambda row: sum(row)) + energy = df[f'ecal_sim_energy_{eng}'].astype(float).to_numpy()#df.eval(f'ecal_sim_energy_{eng}').apply(lambda row: sum(row)) hist, x = np.histogram(energy*1000,bins=np.logspace(0,3,200)) x = x[1:]/2 + x[:-1]/2 plt.plot(x,hist,marker='o',label="ECal") @@ -147,13 +188,18 @@ def rotateY(xdata, zdata, angle): plt.xscale('log') plt.xlim(1e0,1e3) + + + + plt.xlabel('Energy (MeV)') + #plt.savefig('results/Energy_deposition.pdf') -#plt.show() +plt.show() fig4, ax = plt.subplots(2,1,sharex=True,gridspec_kw={'height_ratios': [2,1]}) plt.sca(ax[0]) -plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619+np.array(emean),label='HCal+ECal',fmt='-o') +plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619+np.array(emean),label='HCal/sf+ECal',fmt='-o') plt.errorbar(np.array(Energy)*1000,emean,label='ECal',fmt='-o') plt.legend() plt.yscale('log') @@ -166,7 +212,7 @@ def rotateY(xdata, zdata, angle): plt.xlabel('Truth Energy (MeV)') #plt.savefig('results/Energy_ratio_and_Leakage.pdf') plt.tight_layout() -#plt.show() +plt.show() fig5 = plt.figure() plt.errorbar(np.array(Energy)*1000,np.array(hhits)/1000*100,label='HCal Hits',fmt='-o') @@ -183,7 +229,46 @@ def rotateY(xdata, zdata, angle): plt.ylabel('Fraction of Events with non-zero energy (%)') #plt.savefig('results/Hits.pdf') plt.xscale('log') -#plt.show() +plt.show() + +fig6, ax = plt.subplots(2,3,figsize=(20,10)) +fig6.suptitle('ZDC Clustering') +fig6.tight_layout(pad=1.8) +for i in range(6): + plt.sca(ax[i//3,i%3]) + eng = int(Energy[i]*1000) + + x = df[f'par_x_{eng}'].astype(float).to_numpy() + y = df[f'par_y_{eng}'].astype(float).to_numpy() + z = df[f'par_z_{eng}'].astype(float).to_numpy() + x, z = rotateY(x,z, 0.025) + theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 + condition = theta <= 3.5 + + plt.hist(df[f'ecal_reco_clusters_{eng}'][condition],bins=np.linspace(0,5,6)) + plt.xlabel('Number of Clusters') + plt.title(f'Gamma Energy: {eng} MeV') +plt.show() + +fig7, ax = plt.subplots(2,3,figsize=(20,10)) +fig7.suptitle('ZDC Towering in Clusters') +fig7.tight_layout(pad=1.8) +for i in range(6): + plt.sca(ax[i//3,i%3]) + eng = int(Energy[i]*1000) + + x = df[f'par_x_{eng}'].astype(float).to_numpy() + y = df[f'par_y_{eng}'].astype(float).to_numpy() + z = df[f'par_z_{eng}'].astype(float).to_numpy() + x, z = rotateY(x,z, 0.025) + theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 + condition = theta <= 3.5 + + plt.hist(df[f'ecal_reco_nhits_{eng}'][condition],bins=np.linspace(0,max(df[f'ecal_reco_nhits_{eng}'][condition]),max(df[f'ecal_reco_nhits_{eng}'][condition])+1)) + plt.xlabel('Number of tower in Clusters') + plt.title(f'Gamma Energy: {eng} MeV') +plt.show() + #pdfs = ['results/Energy_reconstruction_cluster.pdf','results/Energy_resolution.pdf','results/Energy_deposition.pdf','results/Energy_ratio_and_Leakage.pdf','results/Hits.pdf'] with PdfPages(f'results/{os.environ["DETECTOR_CONFIG"]}/zdc_lyso/plots.pdf') as pdf: @@ -192,4 +277,5 @@ def rotateY(xdata, zdata, angle): pdf.savefig(fig3) pdf.savefig(fig4) pdf.savefig(fig5) - + pdf.savefig(fig6) + pdf.savefig(fig7) From 9b53b9db03eaae5eff9898926799c7f096220ab6 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Mon, 5 Aug 2024 16:13:49 -0500 Subject: [PATCH 18/19] fix: bench:tracking_performance: snakemake --cores 3 --- benchmarks/tracking_performances/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/tracking_performances/config.yml b/benchmarks/tracking_performances/config.yml index fadeec02..44df84d3 100644 --- a/benchmarks/tracking_performances/config.yml +++ b/benchmarks/tracking_performances/config.yml @@ -18,7 +18,7 @@ bench:tracking_performance: needs: - ["sim:tracking_performance"] script: - - snakemake --cores 1 tracking_performance_local + - snakemake --cores 3 tracking_performance_local collect_results:tracking_performance: extends: .det_benchmark From a07c5f7cc6471cd65912f3eb57cda6241e029df0 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 5 Aug 2024 18:29:11 -0400 Subject: [PATCH 19/19] tracking_performances: upload debug plots as artifacts --- benchmarks/tracking_performances/Snakefile | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/benchmarks/tracking_performances/Snakefile b/benchmarks/tracking_performances/Snakefile index f9f04237..14ae80b3 100644 --- a/benchmarks/tracking_performances/Snakefile +++ b/benchmarks/tracking_performances/Snakefile @@ -122,6 +122,24 @@ root -l -b -q ../{input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX} """ +rule tracking_performance_debug_montage: + input: + expand( + [ + "{{CAMPAIGN}}/Debug_Plots/{SEEDING}/pi-/mom/{SEEDING}_mom_resol_mom{MOMENTUM:.1f}_{{ETA_BIN}}.png", + ], + MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], + SEEDING=["truth", "real"], + ), + output: + "{CAMPAIGN}/Debug_Plots/pi-/mom/mom_resol_debug_{ETA_BIN}.png", + shell: + """ +montage -mode concatenate {input} {output} || true +ls {output} +""" + + TRACKING_PERFORMANCE_ETA_BINS = [-3.5, -2.5, -1.0, 1.0, 2.5, 3.5] rule tracking_performance: @@ -130,6 +148,7 @@ rule tracking_performance: [ "{{CAMPAIGN}}/Final_Results/pi-/mom/mom_resol_{ETA_BIN}.png", "{{CAMPAIGN}}/Final_Results/pi-/mom/mom_resol_{ETA_BIN}.root", + "{{CAMPAIGN}}/Debug_Plots/pi-/mom/mom_resol_debug_{ETA_BIN}.png", ], ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(TRACKING_PERFORMANCE_ETA_BINS[:-1], TRACKING_PERFORMANCE_ETA_BINS[1:])], ),