Merge branch 'master' into pr/tracking_performance_campaigns

.gitlab-ci.yml

@@ -138,6 +138,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'

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:

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
 """
 

benchmarks/tracking_performances/Snakefile

@@ -103,8 +103,14 @@ 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"
 cd {wildcards.CAMPAIGN}
-root -l -b -q ../{input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1.)'
+root -l -b -q ../{input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'$EPIC_VERSION'", "'$EICRECON_VERSION'")'
 """
 
 

benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C

+// 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, 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=========
+   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<double> momV_truth, momV_real;
+   std::vector<double> 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; iptbin<nptbins; ++iptbin){
+    	
+   TCanvas *cp = new TCanvas("cp","cp",1400,1000);
+   cp->SetMargin(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 <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));   
+    histd0xy_truth_1d->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(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
+   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 <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));   
+    histd0xy_real_1d->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
+   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; i<size_truth; i++){
+	pt_truth[i] = momV_truth.at(i);
+	sigma_dcaxy_truth[i] = dcaTresolV_truth.at(i);
+	err_sigma_dcaxy_truth[i] = err_dcaTresolV_truth.at(i);
+	err_pt_truth[i] = 0.;
+	}
+	
+     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; i<size_real; i++){
+	pt_real[i] = momV_real.at(i);
+     sigma_dcaxy_real[i] = dcaTresolV_real.at(i);
+	err_sigma_dcaxy_real[i] = err_dcaTresolV_real.at(i);
+	err_pt_real[i] = 0.;
+	}
+
+     TFile *fout = new TFile(Form("Final_Results/%s/dca/dcaxy_resol_%1.1f_eta_%1.1f.root",particle.Data(),etamin,etamax),"recreate");
+	TGraphErrors *gr1 = new TGraphErrors(size_truth,pt_truth,sigma_dcaxy_truth,err_pt_truth,err_sigma_dcaxy_truth);
+     gr1->SetName("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->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"<<setw(20)<<epic.Data()<<setw(20)<<"EICRecon"<<setw(20)<<eicrecon.Data()<<endl;
+   outfile<<"Etamin"<<setw(20)<<"Etamax"<<setw(20)<<"Pt (GeV/c) \t"<<setw(20)<<"DCAT Resol #mum (Real)"<<endl;
+   for (Int_t i = 0; i<gr1->GetN(); ++i){
+   double x,ytrue, yreal;
+   gr2->GetPoint(i,x,yreal);  
+   outfile<<etamin<<setw(20)<<etamax<<setw(20)<<x<<setw(20)<<yreal<<endl;
+   } 
+    outfile.close();
+     
+     fout->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
+//===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");
+}

benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaz.C

+// 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<double> momV_truth, momV_real;
+   std::vector<double> 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; iptbin<nptbins; ++iptbin){
+    	
+   TCanvas *cp = new TCanvas("cp","cp",1400,1000);
+   cp->SetMargin(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 <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));   
+    histd0z_truth_1d->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 <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));   
+    histd0z_real_1d->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; i<size_truth; i++){
+	pt_truth[i] = momV_truth.at(i);
+	sigma_dcaz_truth[i] = dcaZresolV_truth.at(i);
+	err_sigma_dcaz_truth[i] = err_dcaZresolV_truth.at(i);
+	err_pt_truth[i] = 0.;
+	}
+	
+     const int size_real = momV_real.size();
+	double pt_real[size_real], err_pt_real[size_real], sigma_dcaz_real[size_real], err_sigma_dcaz_real[size_real]; 
+	
+	for (int i=0; i<size_real; i++){
+	pt_real[i] = momV_real.at(i);
+     sigma_dcaz_real[i] = dcaZresolV_real.at(i);
+	err_sigma_dcaz_real[i] = err_dcaZresolV_real.at(i);
+	err_pt_real[i] = 0.;
+	}
+
+     TFile *fout = new TFile(Form("Final_Results/%s/dca/dcaz_resol_%1.1f_eta_%1.1f.root",particle.Data(),etamin,etamax),"recreate");
+	TGraphErrors *gr1 = new TGraphErrors(size_truth,pt_truth,sigma_dcaz_truth,err_pt_truth,err_sigma_dcaz_truth);
+     gr1->SetName("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"<<setw(20)<<epic.Data()<<setw(20)<<"EICRecon"<<setw(20)<<eicrecon.Data()<<endl;
+   outfile<<"Etamin"<<setw(20)<<"Etamax"<<setw(20)<<"Pt (GeV/c) \t"<<setw(20)<<"DCAz Resol #mum (Real)"<<endl;
+   for (Int_t i = 0; i<gr1->GetN(); ++i){
+   double x,ytrue, yreal;
+   gr2->GetPoint(i,x,yreal);  
+   outfile<<etamin<<setw(20)<<etamax<<setw(20)<<x<<setw(20)<<yreal<<endl;
+   } 
+    outfile.close();
+    fout->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");
+}

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 ="", TString eicrecon = "") // 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,20.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();
+   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"<<setw(20)<<epic.Data()<<setw(20)<<"EICrecon"<<setw(20)<<eicrecon.Data()<<endl;
+  outfile<<"Etamin"<<setw(20)<<"Etamax"<<setw(20)<<"p (GeV/c) \t"<<setw(20)<<"Resol  #mum (Truth)"<<setw(20)<<"Resol #mum (Real)"<<endl;
+  for (Int_t i = 0; i<gr1->GetN(); ++i){
+  double x,ytrue, yreal;
+  gr1->GetPoint(i,x,ytrue);    gr2->GetPoint(i,x,yreal);  
+  outfile<<etamin<<setw(20)<<etamax<<setw(20)<<x<<setw(20)<<ytrue<<setw(20)<<yreal<<endl;
+  }
+  outfile.close();
+	
 	fout->cd();
 	mgMom->SetName(Form("mom_resol_%1.1f_eta_%1.1f",etamin,etamax));
 	mgMom->Write();

benchmarks/zdc_lyso/README.md

+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)
+
+

benchmarks/zdc_lyso/Snakefile

+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]}"
+

benchmarks/zdc_lyso/analysis/analysis.py

+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)
+

benchmarks/zdc_lyso/config.yml

+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

benchmarks/zdc_lyso/gen_particles.cxx

+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "TRandom3.h"
+#include "TVector3.h"
+
+#include <iostream>
+#include <random>
+#include <cmath>
+#include <math.h>
+#include <TMath.h>
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+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<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+    GenParticlePtr p2 = std::make_shared<GenParticle>(
+        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<GenParticle>(
+        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<GenVertex>();
+    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;
+}