diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 002bb608794dcb43a4a08fd6e8a9d1263aeff90a..2e836ec66971160841267b79c98dc1d7143119b9 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -145,6 +145,10 @@ include:
- local: 'benchmarks/pid/config.yml'
- local: 'benchmarks/timing/config.yml'
- local: 'benchmarks/b0_tracker/config.yml'
+ - local: 'benchmarks/zdc_sigma/config.yml'
+ - local: 'benchmarks/zdc_lambda/config.yml'
+ - local: 'benchmarks/insert_neutron/config.yml'
+
deploy_results:
allow_failure: true
@@ -159,6 +163,9 @@ deploy_results:
- "collect_results:pid"
- "collect_results:tracking_performance"
- "collect_results:tracking_performance_campaigns"
+ - "collect_results:zdc_sigma"
+ - "collect_results:zdc_lambda"
+ - "collect_results:insert_neutron"
- "collect_results:tracking_performances_dis"
- "collect_results:zdc"
- "collect_results:zdc_lyso"
diff --git a/Snakefile b/Snakefile
index 738bb5e5289c458a0f957fdf02f3a2047800caed..4a12021b8a726b4dc9b2d57084b6088b7b29e21c 100644
--- a/Snakefile
+++ b/Snakefile
@@ -8,6 +8,10 @@ include: "benchmarks/material_scan/Snakefile"
include: "benchmarks/tracking_performances/Snakefile"
include: "benchmarks/tracking_performances_dis/Snakefile"
include: "benchmarks/zdc_lyso/Snakefile"
+include: "benchmarks/zdc_lambda/Snakefile"
+include: "benchmarks/zdc_sigma/Snakefile"
+include: "benchmarks/insert_neutron/Snakefile"
+
use_s3 = config["remote_provider"].lower() == "s3"
use_xrootd = config["remote_provider"].lower() == "xrootd"
diff --git a/benchmarks/insert_neutron/Snakefile b/benchmarks/insert_neutron/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..55e3ff9e09e06786c99d89ff414d5fa5c3e0bf37
--- /dev/null
+++ b/benchmarks/insert_neutron/Snakefile
@@ -0,0 +1,58 @@
+rule insert_neutron_generate:
+ input:
+ script="benchmarks/insert_neutron/analysis/gen_particles.cxx",
+ params:
+ NEVENTS_GEN=1000,
+ th_max=5.7,
+ th_min=2.0
+ output:
+ GEN_FILE="sim_output/insert_neutron/neutron_{P}GeV.hepmc"
+ shell:
+ """
+mkdir -p sim_output/insert_neutron
+root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "neutron", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+"""
+
+rule insert_neutron_simulate:
+ input:
+ GEN_FILE="sim_output/insert_neutron/neutron_{P}GeV.hepmc"
+ params:
+ PHYSICS_LIST="FTFP_BERT"
+ output:
+ SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV.edm4hep.root"
+ shell:
+ """
+NEVENTS_SIM=1000
+# Running simulation
+npsim \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --numberOfEvents $NEVENTS_SIM \
+ --physicsList {params.PHYSICS_LIST} \
+ --inputFiles {input.GEN_FILE} \
+ --outputFile {output.SIM_FILE}
+"""
+
+rule insert_neutron_recon:
+ input:
+ SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV.edm4hep.root"
+ output:
+ REC_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV.edm4eic.root"
+ shell:
+ """
+NEVENTS_REC=1000
+eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalEndcapPInsertRecHits,HcalEndcapPInsertClusters,HcalEndcapPInsertSubcellHits,EcalEndcapPInsertRecHits,EcalEndcapPInsertClusters -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule insert_neutron_analysis:
+ input:
+ expand("sim_output/insert_neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV.edm4eic.root",
+ P=[20, 30, 40, 50, 60, 70, 80],
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ script="benchmarks/insert_neutron/analysis/neutron_plots.py",
+ output:
+ results_dir=directory("results/{DETECTOR_CONFIG}/insert_neutron"),
+ shell:
+ """
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""
diff --git a/benchmarks/insert_neutron/analysis/gen_particles.cxx b/benchmarks/insert_neutron/analysis/gen_particles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0877a85223bbc1af27af91e7dedfb212de7cb20b
--- /dev/null
+++ b/benchmarks/insert_neutron/analysis/gen_particles.cxx
@@ -0,0 +1,127 @@
+#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 = 1000,
+ 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
+ int useCrossingAngle=1 // 0= no rotation, 1 = -25 mrad
+ )
+{
+ 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.*useCrossingAngle; //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;
+}
diff --git a/benchmarks/insert_neutron/analysis/neutron_plots.py b/benchmarks/insert_neutron/analysis/neutron_plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..9314a9745027c3fd77782c09a8636ab8797f6fa0
--- /dev/null
+++ b/benchmarks/insert_neutron/analysis/neutron_plots.py
@@ -0,0 +1,281 @@
+import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys, uproot as ur
+import mplhep as hep
+hep.style.use("CMS")
+
+plt.rcParams['figure.facecolor']='white'
+plt.rcParams['savefig.facecolor']='white'
+plt.rcParams['savefig.bbox']='tight'
+
+plt.rcParams["figure.figsize"] = (7, 7)
+config=sys.argv[1].split("/")[1] #results/{config}/neutron
+outdir=sys.argv[1]+"/"
+try:
+ import os
+ os.mkdir(outdir[:-1])
+except:
+ pass
+
+#read files
+arrays_sim={}
+for p in 20,30,40,50,60,70,80:
+ arrays_sim[p] = ur.open(f'sim_output/insert_neutron/{config}_rec_neutron_{p}GeV.edm4eic.root:events')\
+ .arrays()
+
+def gauss(x, A,mu, sigma):
+ return A * np.exp(-(x-mu)**2/(2*sigma**2))
+
+#get the truth pseudorapidity and energy
+for array in arrays_sim.values():
+ tilt=-0.025
+ px=array['MCParticles.momentum.x'][:,2]
+ py=array['MCParticles.momentum.y'][:,2]
+ pz=array['MCParticles.momentum.z'][:,2]
+ p=np.sqrt(px**2+py**2+pz**2)
+
+ pxp=px*np.cos(tilt)-pz*np.sin(tilt)
+ pyp=py
+ pzp=pz*np.cos(tilt)+px*np.sin(tilt)
+
+ array['E_truth']=np.hypot(p, 0.9406)
+ array['eta_truth']=1/2*np.log((p+pzp)/(p-pzp))
+ array['theta_truth']=np.arccos(pzp/p)
+
+#
+# get reconstructed theta as avg of theta of cluster centers, weighted by energy
+for array in arrays_sim.values():
+ tilt=-0.025
+ x=array['HcalEndcapPInsertClusters.position.x']
+ y=array['HcalEndcapPInsertClusters.position.y']
+ z=array['HcalEndcapPInsertClusters.position.z']
+ E=array['HcalEndcapPInsertClusters.energy']
+ r=np.sqrt(x**2+y**2+z**2)
+
+ xp=x*np.cos(tilt)-z*np.sin(tilt)
+ yp=y
+ zp=z*np.cos(tilt)+x*np.sin(tilt)
+
+ w=E
+
+ array['theta_recon']=np.sum(np.arccos(zp/r)*w, axis=-1)/np.sum(w, axis=-1)
+ array['eta_recon']=-np.log(np.tan(array['theta_recon']/2))
+
+
+#plot theta residuals:
+print("making theta recon plot")
+from scipy.optimize import curve_fit
+
+fig, axs=plt.subplots(1,2, figsize=(16,8))
+plt.sca(axs[0])
+p=40
+eta_min=3.4; eta_max=3.6
+y,x,_=plt.hist(1000*(arrays_sim[p]['theta_recon']-arrays_sim[p]['theta_truth'])\
+ [(arrays_sim[p]['eta_truth']>eta_min)&(arrays_sim[p]['eta_truth']<eta_max)], bins=50,
+ range=(-10,10), histtype='step')
+bc=(x[1:]+x[:-1])/2
+slc=abs(bc)<3
+# try:
+fnc=gauss
+sigma=np.sqrt(y[slc])+(y[slc]==0)
+p0=(100, 0, 5)
+coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+xx=np.linspace(-5,5,100)
+plt.plot(xx,fnc(xx,*coeff))
+# except:
+# pass
+plt.xlabel("$\\theta_{rec}-\\theta_{truth}$ [mrad]")
+plt.ylabel("events")
+plt.title(f"$p={p}$ GeV, ${eta_min}<\\eta<{eta_max}$")
+
+r=[3.2,3.4,3.6,3.8,4.0]
+for eta_min, eta_max in zip(r[:-1],r[1:]):
+ xvals=[]
+ sigmas=[]
+ dsigmas=[]
+ for p in 20,30,40, 50, 60, 70, 80:
+ y,x=np.histogram(1000*(arrays_sim[p]['theta_recon']-arrays_sim[p]['theta_truth'])\
+ [(arrays_sim[p]['eta_truth']>eta_min)&(arrays_sim[p]['eta_truth']<eta_max)],
+ bins=50, range=(-10,10))
+ bc=(x[1:]+x[:-1])/2
+ slc=abs(bc)<3
+ fnc=gauss
+ p0=(100, 0, 5)
+ #print(bc[slc],y[slc])
+ sigma=np.sqrt(y[slc])+(y[slc]==0)
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+ sigmas.append(np.abs(coeff[2]))
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ pass
+ plt.sca(axs[1])
+ plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', label=f"${eta_min}<\\eta<{eta_max}$")
+plt.xlabel("$p_{n}$ [GeV]")
+plt.ylabel("$\\sigma[\\theta]$ [mrad]")
+plt.ylim(0)
+plt.legend()
+plt.tight_layout()
+plt.savefig(outdir+"neutron_theta_recon.pdf")
+
+#now determine the energy recon parameters
+pvals=[]
+resvals=[]
+reserrs=[]
+wvals=[]
+svals=[]
+Euncorr=[]
+
+print("determining the energy recon correction parameters")
+fig,axs=plt.subplots(1,2, figsize=(20,10))
+eta_min=3.4;eta_max=3.6
+for p in 20, 30,40,50,60,70, 80:
+ best_res=1000
+ res_err=1000
+ best_s=1000
+ wrange=np.linspace(30, 70, 41)*0.0257
+ coeff_best=None
+
+ wbest=0
+ a=arrays_sim[p]
+ h=np.sum(a[f'HcalEndcapPInsertClusters.energy'], axis=-1)
+ e=np.sum(a[f'EcalEndcapPInsertClusters.energy'], axis=-1)
+ for w in wrange:
+
+ r=(e/w+h)[(h>0)&(a['eta_truth']>eta_min)&(a['eta_truth']<eta_max)]
+ y,x=np.histogram(r,bins=50)
+ bcs=(x[1:]+x[:-1])/2
+ fnc=gauss
+ slc=abs(bcs-np.mean(r)*1.25)<2*np.std(r)
+ sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
+ p0=(100, np.mean(r), np.std(r))
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+ res=np.abs(coeff[2]/coeff[1])
+
+ if res<best_res:
+ best_res=res
+ coeff_best=coeff
+ res_err=res*np.hypot(np.sqrt(var_matrix[2][2])/coeff[2], np.sqrt(var_matrix[1][1])/coeff[1])
+ wbest=w
+ best_s=np.hypot(p,0.9406)/coeff[1]
+ Euncorr_best=coeff[1]
+ except :
+ print("fit failed")
+
+ if p==50:
+ r=(e/wbest+h)[(h>0)&(a['eta_truth']>3.4)&(a['eta_truth']<3.6)]
+ plt.sca(axs[0])
+ y, x, _= plt.hist(r, histtype='step', bins=50)
+ xx=np.linspace(20, 55, 100)
+ plt.plot(xx,fnc(xx, *coeff_best), ls='-')
+ plt.xlabel("$E_{uncorr}=E_{Hcal}+E_{Ecal}/w$ [GeV]")
+ plt.title(f"p=50 GeV, ${eta_min}<\\eta<{eta_max}$, w={wbest:.2f}")
+ plt.axvline(np.sqrt(50**2+.9406**2), color='g', ls=':')
+ plt.text(40, max(y)*0.9, "generated\nenergy", color='g', fontsize=20)
+
+ Euncorr.append(Euncorr_best)
+ resvals.append(best_res)
+ reserrs.append(res_err)
+ pvals.append(p)
+ svals.append(best_s)
+ wvals.append(wbest)
+
+pvals=np.array(pvals)
+svals=np.array(svals)
+Euncorr=np.array(Euncorr)
+plt.sca(axs[1])
+plt.plot(Euncorr,wvals, label="w")
+w_avg=np.mean(wvals)
+plt.axhline(w_avg, label=f'w avg: {w_avg:.2f}', ls=':')
+plt.plot(Euncorr,svals, label="s")
+m=(np.sum(svals*Euncorr)*len(Euncorr)-np.sum(Euncorr)*np.sum(svals))/(np.sum(Euncorr**2)*len(Euncorr)-np.sum(Euncorr)**2)
+b=np.mean(svals)-np.mean(Euncorr)*m
+plt.plot(Euncorr,Euncorr*m+b, label=f"s fit: ${m:.4f}E_{{uncorr}}+{b:.2f}$", ls=':')
+
+plt.xlabel("$E_{uncorr}=E_{Hcal}+E_{Ecal}/w$ [GeV]")
+plt.title("$E_{n,recon}=s\\times(E_{Hcal}+E_{Ecal}/w)$")
+plt.ylabel('parameter values')
+plt.legend()
+plt.ylim(0)
+plt.savefig(outdir+"neutron_energy_params.pdf")
+
+#now make the energy plot
+print("making energy recon plot")
+fig, axs=plt.subplots(1,3, figsize=(24,8))
+partitions=[3.2,3.4, 3.6, 3.8, 4.0]
+for eta_min, eta_max in zip(partitions[:-1],partitions[1:]):
+ pvals=[]
+ resvals=[]
+ reserrs=[]
+ scalevals=[]
+ scaleerrs=[]
+ for p in 20, 30,40,50,60,70, 80:
+ best_res=1000
+ res_err=1000
+
+
+ wrange=np.linspace(30, 70, 30)*0.0257
+
+ w=w_avg
+ a=arrays_sim[p]
+ h=np.sum(a[f'HcalEndcapPInsertClusters.energy'], axis=-1)
+ e=np.sum(a[f'EcalEndcapPInsertClusters.energy'], axis=-1)
+ #phi=a['phi_truth']
+ uncorr=(e/w+h)
+ s=-0.0064*uncorr+1.80
+ r=uncorr*s #reconstructed energy with correction
+ r=r[(h>0)&(a['eta_truth']>eta_min)&(a['eta_truth']<eta_max)]#&(abs(phi)>np.pi/2)]
+ y,x=np.histogram(r,bins=50)
+ bcs=(x[1:]+x[:-1])/2
+ fnc=gauss
+ slc=abs(bcs-np.mean(r)*1.25)<2*np.std(r)
+ sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
+ p0=(100, np.mean(r), np.std(r))
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+ res=np.abs(coeff[2]/coeff[1])
+
+ if res<best_res:
+ best_res=res
+ res_err=res*np.hypot(np.sqrt(var_matrix[2][2])/coeff[2], np.sqrt(var_matrix[1][1])/coeff[1])
+ wbest=w
+ scale=coeff[1]/np.sqrt(p**2+0.9406**2)
+ dscale=np.sqrt(var_matrix[1][1]/np.sqrt(p**2+0.9406**2))
+ except :
+ print("fit failed")
+ if p==50 and eta_min==3.4:
+ plt.sca(axs[0])
+ plt.errorbar(bcs, y, np.sqrt(y)+(y==0),marker='o', ls='', )
+ plt.title(f'p={p} GeV, ${eta_min}<\\eta<{eta_max}$')
+ plt.xlabel("$E_{recon}$ [GeV]")
+ plt.ylabel("events")
+ #plt.ylim(0)
+ xx=np.linspace(40, 70, 50)
+ plt.plot(xx, fnc(xx, *coeff))
+ resvals.append(best_res)
+ reserrs.append(res_err)
+ scalevals.append(scale)
+ scaleerrs.append(dscale)
+ pvals.append(p)
+ plt.sca(axs[1])
+ plt.errorbar(pvals, resvals, reserrs, marker='o', ls='', label=f"${eta_min}<\\eta<{eta_max}$")
+ #plt.ylim(0)
+ plt.ylabel("$\\sigma[E]/\\mu[E]$")
+ plt.xlabel("$p_{n}$ [GeV]")
+
+ plt.sca(axs[2])
+ plt.errorbar(pvals, scalevals, scaleerrs, marker='o', ls='', label=f"${eta_min}<\\eta<{eta_max}$")
+
+
+ plt.ylabel("$\\mu[E]/E$")
+
+
+axs[2].set_xlabel("$p_n$ [GeV]")
+axs[2].axhline(1, ls='--', color='0.5', alpha=0.7)
+axs[0].set_ylim(0)
+axs[1].set_ylim(0, 0.35)
+axs[2].set_ylim(0)
+axs[1].legend()
+axs[2].legend()
+plt.tight_layout()
+plt.savefig(outdir+"neutron_energy_recon.pdf")
diff --git a/benchmarks/insert_neutron/config.yml b/benchmarks/insert_neutron/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..aac0a7327aac706ea388f4b90d8b1b432ab840e1
--- /dev/null
+++ b/benchmarks/insert_neutron/config.yml
@@ -0,0 +1,35 @@
+sim:insert_neutron:
+ stage: simulate
+ extends: .det_benchmark
+ parallel:
+ matrix:
+ - P: 20
+ - P: 30
+ - P: 40
+ - P: 50
+ - P: 60
+ - P: 70
+ - P: 80
+ timeout: 1 hours
+ script:
+ - snakemake --cores 1 sim_output/insert_neutron/epic_craterlake_rec_neutron_${P}GeV.edm4eic.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:insert_neutron:
+ stage: benchmarks
+ extends: .det_benchmark
+ needs: ["sim:insert_neutron"]
+ script:
+ - snakemake --cores 1 results/epic_craterlake/insert_neutron
+
+collect_results:insert_neutron:
+ stage: collect
+ needs: ["bench:insert_neutron"]
+ script:
+ - ls -al
+ - mv results{,_save}/ # move results directory out of the way to preserve it
+ - snakemake --cores 1 --delete-all-output results/epic_craterlake/insert_neutron
+ - mv results{_save,}/
diff --git a/benchmarks/zdc_lambda/Snakefile b/benchmarks/zdc_lambda/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..21563330bb0683079b22b7099870ff269976c261
--- /dev/null
+++ b/benchmarks/zdc_lambda/Snakefile
@@ -0,0 +1,63 @@
+rule zdc_lambda_generate:
+ input:
+ script="benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx",
+ params:
+ NEVENTS_GEN=100000,
+ output:
+ GEN_FILE="sim_output/zdc_lambda/lambda_decay_{P}GeV.hepmc"
+ shell:
+ """
+root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildcards.P},{wildcards.P})'
+"""
+
+rule zdc_lambda_simulate:
+ input:
+ GEN_FILE="sim_output/zdc_lambda/lambda_decay_{P}GeV.hepmc"
+ params:
+ PHYSICS_LIST="FTFP_BERT"
+ output:
+ SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV.edm4hep.root"
+ shell:
+ """
+if [[ {wildcards.P} -gt 220 ]]; then
+ NEVENTS_SIM=500
+else
+ NEVENTS_SIM=1000
+fi
+# Running simulation
+npsim \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --numberOfEvents $NEVENTS_SIM \
+ --physicsList {params.PHYSICS_LIST} \
+ --inputFiles {input.GEN_FILE} \
+ --outputFile {output.SIM_FILE}
+"""
+
+rule zdc_lambda_recon:
+ input:
+ SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV.edm4hep.root"
+ output:
+ REC_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV.edm4eic.root"
+ shell:
+ """
+if [[ {wildcards.P} -gt 220 ]]; then
+ NEVENTS_REC=500
+else
+ NEVENTS_REC=1000
+fi
+eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalFarForwardZDCClusters,HcalFarForwardZDCRecHits,HcalFarForwardZDCSubcellHits -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule zdc_lambda_analysis:
+ input:
+ expand("sim_output/zdc_lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV.edm4eic.root",
+ P=[100, 125, 150,175, 200, 225, 250, 275],
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ script="benchmarks/zdc_lambda/analysis/lambda_plots.py",
+ output:
+ results_dir=directory("results/{DETECTOR_CONFIG}/zdc_lambda"),
+ shell:
+ """
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""
diff --git a/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx b/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..567eda5bcf5497f1b64596745ef5820e7a54eff9
--- /dev/null
+++ b/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx
@@ -0,0 +1,239 @@
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "TRandom3.h"
+#include "TVector3.h"
+
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+#include <iostream>
+#include <random>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+std::tuple<double, int, double> GetParticleInfo(TDatabasePDG* pdg, TString particle_name)
+{
+ TParticlePDG *particle = pdg->GetParticle(particle_name);
+ const double mass = particle->Mass();
+ const int pdgID = particle->PdgCode();
+ const double lifetime = particle->Lifetime();
+ return std::make_tuple(mass, pdgID, lifetime);
+}
+// Calculates the decay length of a particle. Samples from an exponential decay.
+double GetDecayLength(TRandom3* r1, double lifetime, double mass, double momentum_magnitude)
+{
+ double c_speed = TMath::C() * 1000.; // speed of light im mm/sec
+ double average_decay_length = (momentum_magnitude/mass) * lifetime * c_speed;
+ return r1->Exp(average_decay_length);
+}
+
+// Generate single lambda mesons and decay them to a neutron + 2 photons
+void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "lambda_decay.hepmc",
+ double p_min = 100., // in GeV/c
+ double p_max = 275.) // in GeV/c
+{
+
+ const double theta_min = 0.0; // in mRad
+ const double theta_max = 3.0; // in mRad
+ //const double p_min = 100.; // in GeV/c
+ //const double p_max = 275.; // in GeV/c
+
+ WriterAscii hepmc_output(out_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // Random number generator
+ TRandom3 *r1 = new TRandom3(seed); //Default = 0, which uses clock to set seed
+ cout<<"Random number seed is "<<r1->GetSeed()<<"!"<<endl;
+
+ // Getting generated particle information
+ TDatabasePDG *pdg = new TDatabasePDG();
+
+ auto lambda_info = GetParticleInfo(pdg, "Lambda0");
+ double lambda_mass = std::get<0>(lambda_info);
+ int lambda_pdgID = std::get<1>(lambda_info);
+ double lambda_lifetime = std::get<2>(lambda_info);
+
+ auto neutron_info = GetParticleInfo(pdg, "neutron");
+ double neutron_mass = std::get<0>(neutron_info);
+ int neutron_pdgID = std::get<1>(neutron_info);
+
+ auto pi0_info = GetParticleInfo(pdg, "pi0");
+ double pi0_mass = std::get<0>(pi0_info);
+ int pi0_pdgID = std::get<1>(pi0_info);
+ double pi0_lifetime = std::get<2>(pi0_info);
+
+ auto photon_info = GetParticleInfo(pdg, "gamma");
+ double photon_mass = std::get<0>(photon_info);
+ int photon_pdgID = std::get<1>(photon_info);
+
+ 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 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 EIC proton beam direction
+ Double_t lambda_p = r1->Uniform(p_min, p_max);
+ Double_t lambda_phi = r1->Uniform(0.0, 2.0 * M_PI);
+ Double_t lambda_th = r1->Uniform(theta_min/1000., theta_max/1000.); // Divide by 1000 for radians
+ Double_t lambda_px = lambda_p * TMath::Cos(lambda_phi) * TMath::Sin(lambda_th);
+ Double_t lambda_py = lambda_p * TMath::Sin(lambda_phi) * TMath::Sin(lambda_th);
+ Double_t lambda_pz = lambda_p * TMath::Cos(lambda_th);
+ Double_t lambda_E = TMath::Sqrt(lambda_p*lambda_p + lambda_mass*lambda_mass);
+
+ // Rotate to lab coordinate system
+ TVector3 lambda_pvec(lambda_px, lambda_py, lambda_pz);
+ double cross_angle = -25./1000.; // in Rad
+ TVector3 pbeam_dir(TMath::Sin(cross_angle), 0, TMath::Cos(cross_angle)); //proton beam direction
+ lambda_pvec.RotateY(cross_angle); // Theta is returned positive, beam in negative X
+
+ // type 2 is state that will decay
+ GenParticlePtr p_lambda = std::make_shared<GenParticle>(
+ FourVector(lambda_pvec.X(), lambda_pvec.Y(), lambda_pvec.Z(), lambda_E), lambda_pdgID, 2 );
+
+ // Generating lambda particle, will be generated at origin
+ // Must have input electron + proton for vertex
+ GenVertexPtr lambda_initial_vertex = std::make_shared<GenVertex>();
+ lambda_initial_vertex->add_particle_in(p1);
+ lambda_initial_vertex->add_particle_in(p2);
+ lambda_initial_vertex->add_particle_out(p_lambda);
+ evt.add_vertex(lambda_initial_vertex);
+
+ // Generate neutron + pi0 in lambda rest frame
+ TLorentzVector neutron_rest, pi0_rest;
+
+ // Generating uniformly along a sphere
+ double cost_neutron_rest = r1->Uniform(-1,1);
+ double th_neutron_rest = TMath::ACos(cost_neutron_rest);
+ double sint_neutron_rest = TMath::Sin(th_neutron_rest);
+
+ double phi_neutron_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
+ double cosp_neutron_rest = TMath::Cos(phi_neutron_rest);
+ double sinp_neutron_rest = TMath::Sin(phi_neutron_rest);
+
+ // Calculate energy of each particle in the lambda rest frame
+ // See problem 3.19 in Introduction to Elementary Particles, 2nd edition by D. Griffiths
+ double E_neutron_rest = (-TMath::Power(pi0_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(neutron_mass, 2.) ) / (2. * lambda_mass) ;
+ double E_pi0_rest = (-TMath::Power(neutron_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(pi0_mass, 2.) ) / (2. * lambda_mass) ;
+
+ // Both particles will have the same momentum, so just use neutron variables
+ double momentum_rest = TMath::Sqrt( E_neutron_rest*E_neutron_rest - neutron_mass*neutron_mass );
+
+ neutron_rest.SetE(E_neutron_rest);
+ neutron_rest.SetPx( momentum_rest * sint_neutron_rest * cosp_neutron_rest );
+ neutron_rest.SetPy( momentum_rest * sint_neutron_rest * sinp_neutron_rest );
+ neutron_rest.SetPz( momentum_rest * cost_neutron_rest );
+
+ pi0_rest.SetE(E_pi0_rest);
+ pi0_rest.SetPx( -neutron_rest.Px() );
+ pi0_rest.SetPy( -neutron_rest.Py() );
+ pi0_rest.SetPz( -neutron_rest.Pz() );
+
+ // Boost neutron & pion to lab frame
+ TLorentzVector lambda_lab(lambda_pvec.X(), lambda_pvec.Y(), lambda_pvec.Z(), lambda_E);
+ TVector3 lambda_boost = lambda_lab.BoostVector();
+ TLorentzVector neutron_lab, pi0_lab;
+ neutron_lab = neutron_rest;
+ neutron_lab.Boost(lambda_boost);
+ pi0_lab = pi0_rest;
+ pi0_lab.Boost(lambda_boost);
+
+ // Calculating position for lambda decay
+ TVector3 lambda_unit = lambda_lab.Vect().Unit();
+ double lambda_decay_length = GetDecayLength(r1, lambda_lifetime, lambda_mass, lambda_lab.P());
+ TVector3 lambda_decay_position = lambda_unit * lambda_decay_length;
+ double lambda_decay_time = lambda_decay_length / lambda_lab.Beta() ; // Decay time in lab frame in length units (mm)
+
+ // Generating vertex for lambda decay
+ GenParticlePtr p_neutron = std::make_shared<GenParticle>(
+ FourVector(neutron_lab.Px(), neutron_lab.Py(), neutron_lab.Pz(), neutron_lab.E()), neutron_pdgID, 1 );
+
+ GenParticlePtr p_pi0 = std::make_shared<GenParticle>(
+ FourVector(pi0_lab.Px(), pi0_lab.Py(), pi0_lab.Pz(), pi0_lab.E()), pi0_pdgID, 2 );
+
+ GenVertexPtr v_lambda_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
+ v_lambda_decay->add_particle_in(p_lambda);
+ v_lambda_decay->add_particle_out(p_neutron);
+ v_lambda_decay->add_particle_out(p_pi0);
+
+ evt.add_vertex(v_lambda_decay);
+
+ // Generate two photons from pi0 decay
+ TLorentzVector gamma1_rest, gamma2_rest;
+
+ // Generating uniformly along a sphere
+ double cost_gamma1_rest = r1->Uniform(-1,1);
+ double th_gamma1_rest = TMath::ACos(cost_gamma1_rest);
+ double sint_gamma1_rest = TMath::Sin(th_gamma1_rest);
+
+ double phi_gamma1_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
+ double cosp_gamma1_rest = TMath::Cos(phi_gamma1_rest);
+ double sinp_gamma1_rest = TMath::Sin(phi_gamma1_rest);
+
+ // Photons are massless so they each get equal energies
+ gamma1_rest.SetE(pi0_mass/2.);
+ gamma1_rest.SetPx( (pi0_mass/2.)*sint_gamma1_rest*cosp_gamma1_rest );
+ gamma1_rest.SetPy( (pi0_mass/2.)*sint_gamma1_rest*sinp_gamma1_rest );
+ gamma1_rest.SetPz( (pi0_mass/2.)*cost_gamma1_rest );
+
+ gamma2_rest.SetE(pi0_mass/2.);
+ gamma2_rest.SetPx( -gamma1_rest.Px() );
+ gamma2_rest.SetPy( -gamma1_rest.Py() );
+ gamma2_rest.SetPz( -gamma1_rest.Pz() );
+
+ // Boost neutron & pion to lab frame
+ TVector3 pi0_boost = pi0_lab.BoostVector();
+ TLorentzVector gamma1_lab, gamma2_lab;
+ gamma1_lab = gamma1_rest;
+ gamma1_lab.Boost(pi0_boost);
+ gamma2_lab = gamma2_rest;
+ gamma2_lab.Boost(pi0_boost);
+
+ GenParticlePtr p_gamma1 = std::make_shared<GenParticle>(
+ FourVector(gamma1_lab.Px(), gamma1_lab.Py(), gamma1_lab.Pz(), gamma1_lab.E()), photon_pdgID, 1 );
+
+ GenParticlePtr p_gamma2 = std::make_shared<GenParticle>(
+ FourVector(gamma2_lab.Px(), gamma2_lab.Py(), gamma2_lab.Pz(), gamma2_lab.E()), photon_pdgID, 1 );
+
+ // Generate pi0 at same position as the lambda. Approximating pi0 decay as instantaneous
+ GenVertexPtr v_pi0_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
+ v_pi0_decay->add_particle_in(p_pi0);
+ v_pi0_decay->add_particle_out(p_gamma1);
+ v_pi0_decay->add_particle_out(p_gamma2);
+
+ //std::cout<< lambda_pvec.Angle(pbeam_dir) << " " << neutron_lab.Angle(pbeam_dir) << " " << gamma1_lab.Angle(pbeam_dir) << " " << gamma2_lab.Angle(pbeam_dir) << std::endl;
+
+ evt.add_vertex(v_pi0_decay);
+
+ if (events_parsed == 0) {
+ std::cout << "First event: " << std::endl;
+ Print::listing(evt);
+ }
+ double zdc_z=35800;
+ TVector3 extrap_neutron=lambda_decay_position+neutron_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(neutron_lab.Vect())));
+ TVector3 extrap_gamma1=lambda_decay_position+gamma1_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma1_lab.Vect())));
+ TVector3 extrap_gamma2=lambda_decay_position+gamma2_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma2_lab.Vect())));
+ if (extrap_neutron.Angle(pbeam_dir)<0.004 && extrap_gamma1.Angle(pbeam_dir)<0.004 && extrap_gamma2.Angle(pbeam_dir)<0.004 && lambda_decay_position.Dot(pbeam_dir)<zdc_z)
+ hepmc_output.write_event(evt);
+ if (events_parsed % 1000 == 0) {
+ std::cout << "Event: " << events_parsed << std::endl;
+ }
+ evt.clear();
+ }
+ hepmc_output.close();
+
+ std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}
diff --git a/benchmarks/zdc_lambda/analysis/lambda_plots.py b/benchmarks/zdc_lambda/analysis/lambda_plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..88aff7f38b130967a0062ec5d193562e50ced06a
--- /dev/null
+++ b/benchmarks/zdc_lambda/analysis/lambda_plots.py
@@ -0,0 +1,372 @@
+import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys
+import mplhep as hep
+hep.style.use("CMS")
+
+plt.rcParams['figure.facecolor']='white'
+plt.rcParams['savefig.facecolor']='white'
+plt.rcParams['savefig.bbox']='tight'
+
+plt.rcParams["figure.figsize"] = (7, 7)
+
+outdir=sys.argv[1]+"/"
+config=outdir.split("/")[1]
+try:
+ import os
+ os.mkdir(outdir[:-1])
+except:
+ pass
+import uproot as ur
+arrays_sim={}
+momenta=100, 125, 150, 175,200,225,250,275
+for p in momenta:
+ filename=f'sim_output/zdc_lambda/{config}_rec_lambda_dec_{p}GeV.edm4eic.root'
+ print("opening file", filename)
+ events = ur.open(filename+':events')
+ arrays_sim[p] = events.arrays()[:-1] #remove last event, which for some reason is blank
+ import gc
+ gc.collect()
+ print("read", filename)
+
+def gauss(x, A,mu, sigma):
+ return A * np.exp(-(x-mu)**2/(2*sigma**2))
+
+#keep track of the number of clusters per event
+nclusters={}
+
+for p in momenta:
+ plt.figure()
+ nclusters[p]=[]
+ for i in range(len(arrays_sim[p])):
+ nclusters[p].append(len(arrays_sim[p]["HcalFarForwardZDCClusters.position.x"][i]))
+ nclusters[p]=np.array(nclusters[p])
+ plt.hist(nclusters[p],bins=20, range=(0,20))
+ plt.xlabel("number of clusters")
+ plt.yscale('log')
+ plt.title(f"$p_\Lambda={p}$ GeV")
+ plt.ylim(1)
+ plt.savefig(outdir+f"nclust_{p}GeV_recon.pdf")
+ print("saved file ", outdir+f"nclust_{p}GeV_recon.pdf")
+
+
+
+pt_truth={}
+theta_truth={}
+
+for p in momenta:
+ #get the truth value of theta* and pt*
+ px=arrays_sim[p]["MCParticles.momentum.x"][:,2]
+ py=arrays_sim[p]["MCParticles.momentum.y"][:,2]
+ pz=arrays_sim[p]["MCParticles.momentum.z"][:,2]
+ tilt=-0.025
+ pt_truth[p]=np.hypot(px*np.cos(tilt)-pz*np.sin(tilt), py)
+ theta_truth[p]=np.arctan2(pt_truth[p],pz*np.cos(tilt)+px*np.sin(tilt))
+
+
+#create an array with the same shape as the cluster-level arrays
+is_neutron_cand={}
+for p in momenta:
+ is_neutron_cand[p]=(0*arrays_sim[p][f"HcalFarForwardZDCClusters.energy"]).to_list()
+
+ #largest_eigenvalues
+ for i in range(len(arrays_sim[p])):
+ pars=arrays_sim[p]["_HcalFarForwardZDCClusters_shapeParameters"][i]
+ index_of_max=-1
+ max_val=0
+ eigs=[]
+ shape_params_per_cluster=7
+ for j in range(len(pars)//shape_params_per_cluster):
+ largest_eigenvalue=max(pars[shape_params_per_cluster*j+4:shape_params_per_cluster*j+7])
+ eigs.append(largest_eigenvalue)
+ if(largest_eigenvalue>max_val):
+ max_val=largest_eigenvalue
+ index_of_max=j
+ if index_of_max >=0:
+ is_neutron_cand[p][i][index_of_max]=1
+ eigs.sort()
+
+ is_neutron_cand[p]=ak.Array(is_neutron_cand[p])
+
+
+#with the position of the vertex determined by assuming the mass of the pi0
+#corrected pt* and theta* recon
+pt_recon_corr={}
+theta_recon_corr={}
+mass_recon_corr={}
+mass_pi0_recon_corr={}
+pi0_converged={}
+zvtx_recon={}
+
+maxZ=35800
+for p in momenta:
+ xvtx=0
+ yvtx=0
+ zvtx=0
+
+ for iteration in range(20):
+
+ #compute the value of theta* using the clusters in the ZDC
+ xc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.x"]
+ yc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.y"]
+ zc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.z"]
+ E=arrays_sim[p][f"HcalFarForwardZDCClusters.energy"]
+ #apply correction to the neutron candidates only
+ A,B,C=-0.0756, -1.91, 2.30
+ neutron_corr=(1-is_neutron_cand[p])+is_neutron_cand[p]/(1+A+B/np.sqrt(E)+C/E)
+ E=E*neutron_corr
+
+ E_recon=np.sum(E, axis=-1)
+ pabs=np.sqrt(E**2-is_neutron_cand[p]*.9406**2)
+ tilt=-0.025
+ xcp=xc*np.cos(tilt)-zc*np.sin(tilt)
+ ycp=yc
+ zcp=zc*np.cos(tilt)+xc*np.sin(tilt)
+ rcp=np.sqrt(xcp**2+ycp**2+zcp**2)
+
+ ux=(xcp-xvtx)
+ uy=(ycp-yvtx)
+ uz=(zcp-zvtx)
+
+ norm=np.sqrt(ux**2+uy**2+uz**2)
+ ux=ux/norm
+ uy=uy/norm
+ uz=uz/norm
+
+ px_recon,py_recon,pz_recon=np.sum(pabs*ux, axis=-1),np.sum(pabs*uy, axis=-1),np.sum(pabs*uz, axis=-1)
+
+ pt_recon_corr[p]=np.hypot(px_recon,py_recon)
+ theta_recon_corr[p]=np.arctan2(pt_recon_corr[p], pz_recon)
+
+ mass_recon_corr[p]=np.sqrt((E_recon)**2\
+ -(px_recon)**2\
+ -(py_recon)**2\
+ -(pz_recon)**2)
+ mass_pi0_recon_corr[p]=np.sqrt(np.sum(pabs*(1-is_neutron_cand[p]), axis=-1)**2\
+ -np.sum(pabs*ux*(1-is_neutron_cand[p]), axis=-1)**2\
+ -np.sum(pabs*uy*(1-is_neutron_cand[p]), axis=-1)**2\
+ -np.sum(pabs*uz*(1-is_neutron_cand[p]), axis=-1)**2)
+ alpha=1
+ if iteration==0:
+ u=np.sqrt(px_recon**2+py_recon**2+pz_recon**2)
+ ux=px_recon/u
+ uy=py_recon/u
+ uz=pz_recon/u
+ zeta=1/2
+ zvtx=maxZ*np.power(zeta,alpha)
+ xvtx=ux/uz*zvtx
+ yvtx=uy/uz*zvtx
+ else :
+ u=np.sqrt(px_recon**2+py_recon**2+pz_recon**2)
+ ux=px_recon/u
+ uy=py_recon/u
+ uz=pz_recon/u
+ s=2*(mass_pi0_recon_corr[p]<0.135)-1
+ zeta=np.power(zvtx/maxZ, 1/alpha)
+ zeta=zeta+s*1/2**(1+iteration)
+ zvtx=maxZ*np.power(zeta,alpha)
+ xvtx=ux/uz*zvtx
+ yvtx=uy/uz*zvtx
+ #print(zvtx)
+ pi0_converged[p]=np.abs(mass_pi0_recon_corr[p]-0.135)<0.01
+ zvtx_recon[p]=zvtx
+
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+plt.sca(axs[0])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+x=[]
+y=[]
+for p in momenta:
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ x+=list(theta_truth[p][accept]*1000)
+ y+=list(theta_recon_corr[p][accept]*1000)
+plt.scatter(x,y)
+plt.xlabel("$\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]")
+plt.ylabel("$\\theta^{*\\rm recon}_{\\Lambda}$ [mrad]")
+plt.xlim(0,3.2)
+plt.ylim(0,3.2)
+
+plt.sca(axs[1])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+y,x,_=plt.hist(y-np.array(x), bins=50, range=(-1,1))
+bc=(x[1:]+x[:-1])/2
+
+from scipy.optimize import curve_fit
+slc=abs(bc)<0.3
+fnc=gauss
+p0=[100, 0, 0.05]
+coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
+ sigma=np.sqrt(y[slc])+(y[slc]==0))
+x=np.linspace(-1, 1)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+plt.xlabel("$\\theta^{*\\rm recon}_{\\Lambda}-\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]")
+plt.ylabel("events")
+
+plt.sca(axs[2])
+sigmas=[]
+dsigmas=[]
+xvals=[]
+for p in momenta:
+
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ y,x=np.histogram((theta_recon_corr[p]-theta_truth[p])[accept]*1000, bins=100, range=(-1,1))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc)<0.3
+ fnc=gauss
+ p0=(100, 0, 0.06)
+ #print(bc[slc],y[slc])
+ sigma=np.sqrt(y[slc])+(y[slc]==0)
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+ sigmas.append(coeff[2])
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print("fit failed")
+plt.ylim(0, 0.3)
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+x=np.linspace(100, 275, 100)
+plt.plot(x, 3/np.sqrt(x), color='tab:orange')
+plt.text(170, .23, "YR requirement:\n 3 mrad/$\\sqrt{E}$")
+plt.xlabel("$E_{\\Lambda}$ [GeV]")
+plt.ylabel("$\\sigma[\\theta^*_{\\Lambda}]$ [mrad]")
+plt.tight_layout()
+plt.savefig(outdir+"thetastar_recon.pdf")
+#plt.show()
+
+
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+plt.sca(axs[0])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+x=[]
+y=[]
+for p in momenta:
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ x+=list(arrays_sim[p]['MCParticles.vertex.z'][:,3][accept]/1000)
+ y+=list(zvtx_recon[p][accept]/1000)
+plt.scatter(x,y)
+#print(x,y)
+plt.xlabel("$z^{\\rm truth}_{\\rm vtx}$ [m]")
+plt.ylabel("$z^{\\rm recon}_{\\rm vtx}$ [m]")
+plt.xlim(0,40)
+plt.ylim(0,40)
+
+plt.sca(axs[1])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+y,x,_=plt.hist(y-np.array(x), bins=50, range=(-10,10))
+bc=(x[1:]+x[:-1])/2
+
+from scipy.optimize import curve_fit
+slc=abs(bc)<5
+fnc=gauss
+p0=[100, 0, 1]
+coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
+ sigma=np.sqrt(y[slc])+(y[slc]==0))
+x=np.linspace(-5, 5)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+print(coeff[2], np.sqrt(var_matrix[2][2]))
+plt.xlabel("$z^{*\\rm recon}_{\\rm vtx}-z^{*\\rm truth}_{\\rm vtx}$ [m]")
+plt.ylabel("events")
+
+plt.sca(axs[2])
+sigmas=[]
+dsigmas=[]
+xvals=[]
+for p in momenta:
+
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ a=list((zvtx_recon[p]-arrays_sim[p]['MCParticles.vertex.z'][:,3])[accept]/1000)
+ y,x=np.histogram(a, bins=100, range=(-10,10))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc)<5
+ fnc=gauss
+ p0=(100, 0, 1)
+ #print(bc[slc],y[slc])
+ sigma=np.sqrt(y[slc])+(y[slc]==0)
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+ sigmas.append(coeff[2])
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print("fit failed")
+plt.ylim(0, 2)
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+x=np.linspace(100, 275, 100)
+
+avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
+plt.axhline(avg, color='tab:orange')
+plt.text(150, 1.25,f"$\\sigma\\approx${avg:.1f} m")
+
+plt.xlabel("$E_{\\Lambda}$ [GeV]")
+plt.ylabel("$\\sigma[z_{\\rm vtx}]$ [m]")
+plt.tight_layout()
+plt.savefig(outdir+"zvtx_recon.pdf")
+#plt.show()
+
+p=100
+fig,axs=plt.subplots(1,2, figsize=(16, 8))
+plt.sca(axs[0])
+lambda_mass=1.115683
+vals=[]
+for p in momenta:
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ vals+=list(mass_recon_corr[p][accept])
+
+y,x,_= plt.hist(vals, bins=100, range=(1.0, 1.25))
+bc=(x[1:]+x[:-1])/2
+plt.axvline(lambda_mass, ls='--', color='tab:green', lw=3)
+plt.text(lambda_mass+.01, np.max(y)*1.05, "PDG mass", color='tab:green')
+plt.xlabel("$m_{\\Lambda}^{\\rm recon}$ [GeV]")
+plt.ylim(0, np.max(y)*1.2)
+plt.xlim(1.0, 1.25)
+
+from scipy.optimize import curve_fit
+slc=abs(bc-lambda_mass)<0.07
+fnc=gauss
+p0=[100, lambda_mass, 0.04]
+coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
+ sigma=np.sqrt(y[slc])+(y[slc]==0))
+x=np.linspace(0.8, 1.3, 200)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+print(coeff[2], np.sqrt(var_matrix[2][2]))
+plt.xlabel("$m^{\\rm recon}_{\\Lambda}$ [GeV]")
+plt.ylabel("events")
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+
+plt.sca(axs[1])
+xvals=[]
+sigmas=[]
+dsigmas=[]
+for p in momenta:
+ accept=(nclusters[p]==3) &(pi0_converged[p])
+ y,x= np.histogram(mass_recon_corr[p][accept], bins=100, range=(0.6,1.4))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc-lambda_mass)<0.07
+ fnc=gauss
+ p0=[100, lambda_mass, 0.05]
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+ x=np.linspace(0.8, 1.3, 200)
+ sigmas.append(coeff[2])
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print("fit failed")
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
+plt.axhline(avg, color='tab:orange')
+plt.text(150, 0.01,f"$\\sigma\\approx${avg*1000:.0f} MeV")
+plt.xlabel("$E_{\\Lambda}$ [GeV]")
+plt.ylabel("$\\sigma[m_{\\Lambda}]$ [GeV]")
+plt.ylim(0, 0.02)
+plt.tight_layout()
+plt.savefig(outdir+"lambda_mass_rec.pdf")
diff --git a/benchmarks/zdc_lambda/config.yml b/benchmarks/zdc_lambda/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..396b192491e929267a8d56cc4809bab37cdb1987
--- /dev/null
+++ b/benchmarks/zdc_lambda/config.yml
@@ -0,0 +1,36 @@
+sim:zdc_lambda:
+ stage: simulate
+ extends: .det_benchmark
+ parallel:
+ matrix:
+ - P: 100
+ - P: 125
+ - P: 150
+ - P: 175
+ - P: 200
+ - P: 225
+ - P: 250
+ - P: 275
+ timeout: 1 hours
+ script:
+ - snakemake --cores 1 sim_output/zdc_lambda/epic_zdc_sipm_on_tile_only_rec_lambda_dec_${P}GeV.edm4eic.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:zdc_lambda:
+ stage: benchmarks
+ extends: .det_benchmark
+ needs: ["sim:zdc_lambda"]
+ script:
+ - snakemake --cores 1 results/epic_zdc_sipm_on_tile_only/zdc_lambda
+
+collect_results:zdc_lambda:
+ stage: collect
+ needs: ["bench:zdc_lambda"]
+ script:
+ - ls -al
+ - mv results{,_save}/ # move results directory out of the way to preserve it
+ - snakemake --cores 1 --delete-all-output results/epic_zdc_sipm_on_tile_only/zdc_lambda
+ - mv results{_save,}/
diff --git a/benchmarks/zdc_sigma/Snakefile b/benchmarks/zdc_sigma/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..a92a37819bae72e8842ada3bb1a2995fe37694af
--- /dev/null
+++ b/benchmarks/zdc_sigma/Snakefile
@@ -0,0 +1,63 @@
+rule zdc_sigma_generate:
+ input:
+ script="benchmarks/zdc_sigma/analysis/gen_sigma_decay.cxx",
+ params:
+ NEVENTS_GEN=100000,
+ output:
+ GEN_FILE="sim_output/zdc_sigma/sigma_decay_{P}GeV.hepmc"
+ shell:
+ """
+root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildcards.P},{wildcards.P})'
+"""
+
+rule zdc_sigma_simulate:
+ input:
+ GEN_FILE="sim_output/zdc_sigma/sigma_decay_{P}GeV.hepmc"
+ params:
+ PHYSICS_LIST="FTFP_BERT"
+ output:
+ SIM_FILE="sim_output/zdc_sigma/{DETECTOR_CONFIG}_sim_sigma_dec_{P}GeV.edm4hep.root"
+ shell:
+ """
+if [[ {wildcards.P} -gt 220 ]]; then
+ NEVENTS_SIM=500
+else
+ NEVENTS_SIM=1000
+fi
+# Running simulation
+npsim \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --numberOfEvents $NEVENTS_SIM \
+ --physicsList {params.PHYSICS_LIST} \
+ --inputFiles {input.GEN_FILE} \
+ --outputFile {output.SIM_FILE}
+"""
+
+rule zdc_sigma_recon:
+ input:
+ SIM_FILE="sim_output/zdc_sigma/{DETECTOR_CONFIG}_sim_sigma_dec_{P}GeV.edm4hep.root"
+ output:
+ REC_FILE="sim_output/zdc_sigma/{DETECTOR_CONFIG}_rec_sigma_dec_{P}GeV.edm4eic.root"
+ shell:
+ """
+if [[ {wildcards.P} -gt 220 ]]; then
+ NEVENTS_REC=500
+else
+ NEVENTS_REC=1000
+fi
+eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalFarForwardZDCClusters,HcalFarForwardZDCRecHits,HcalFarForwardZDCSubcellHits -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule zdc_sigma_analysis:
+ input:
+ expand("sim_output/zdc_sigma/{DETECTOR_CONFIG}_rec_sigma_dec_{P}GeV.edm4eic.root",
+ P=[100, 125, 150,175, 200, 225, 250, 275],
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ script="benchmarks/zdc_sigma/analysis/sigma_plots.py",
+ output:
+ results_dir=directory("results/{DETECTOR_CONFIG}/zdc_sigma"),
+ shell:
+ """
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""
diff --git a/benchmarks/zdc_sigma/analysis/gen_sigma_decay.cxx b/benchmarks/zdc_sigma/analysis/gen_sigma_decay.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4da92801257e861bca80cc6638722fff360019c7
--- /dev/null
+++ b/benchmarks/zdc_sigma/analysis/gen_sigma_decay.cxx
@@ -0,0 +1,305 @@
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "TRandom3.h"
+#include "TVector3.h"
+
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+#include <iostream>
+#include <random>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+std::tuple<double, int, double> GetParticleInfo(TDatabasePDG* pdg, TString particle_name)
+{
+ TParticlePDG *particle = pdg->GetParticle(particle_name);
+ const double mass = particle->Mass();
+ const int pdgID = particle->PdgCode();
+ const double lifetime = particle->Lifetime();
+ return std::make_tuple(mass, pdgID, lifetime);
+}
+// Calculates the decay length of a particle. Samples from an exponential decay.
+double GetDecayLength(TRandom3* r1, double lifetime, double mass, double momentum_magnitude)
+{
+ double c_speed = TMath::C() * 1000.; // speed of light im mm/sec
+ double average_decay_length = (momentum_magnitude/mass) * lifetime * c_speed;
+ return r1->Exp(average_decay_length);
+}
+
+// Generate single sigma baryons and decay them to a neutron + 2 photons
+void gen_sigma_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "sigma_decay.hepmc",
+ double p_min = 100., // in GeV/c
+ double p_max = 275.) // in GeV/c
+{
+ int accepted_events=0;
+ const double theta_min = 0.0; // in mRad
+ const double theta_max = 3.0; // in mRad
+ //const double p_min = 100.; // in GeV/c
+ //const double p_max = 275.; // in GeV/c
+
+ WriterAscii hepmc_output(out_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // Random number generator
+ TRandom3 *r1 = new TRandom3(seed); //Default = 0, which uses clock to set seed
+ cout<<"Random number seed is "<<r1->GetSeed()<<"!"<<endl;
+
+ // Getting generated particle information
+ TDatabasePDG *pdg = new TDatabasePDG();
+
+ auto sigma_info = GetParticleInfo(pdg, "Sigma0");
+ double sigma_mass = std::get<0>(sigma_info);
+ int sigma_pdgID = std::get<1>(sigma_info);
+ double sigma_lifetime = std::get<2>(sigma_info);
+
+ auto lambda_info = GetParticleInfo(pdg, "Lambda0");
+ double lambda_mass = std::get<0>(lambda_info);
+ int lambda_pdgID = std::get<1>(lambda_info);
+ double lambda_lifetime = std::get<2>(lambda_info);
+
+ auto neutron_info = GetParticleInfo(pdg, "neutron");
+ double neutron_mass = std::get<0>(neutron_info);
+ int neutron_pdgID = std::get<1>(neutron_info);
+
+ auto pi0_info = GetParticleInfo(pdg, "pi0");
+ double pi0_mass = std::get<0>(pi0_info);
+ int pi0_pdgID = std::get<1>(pi0_info);
+ double pi0_lifetime = std::get<2>(pi0_info);
+
+ auto photon_info = GetParticleInfo(pdg, "gamma");
+ double photon_mass = std::get<0>(photon_info);
+ int photon_pdgID = std::get<1>(photon_info);
+
+ 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 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 EIC proton beam direction
+ Double_t sigma_p = r1->Uniform(p_min, p_max);
+ Double_t sigma_phi = r1->Uniform(0.0, 2.0 * M_PI);
+ Double_t sigma_th = r1->Uniform(theta_min/1000., theta_max/1000.); // Divide by 1000 for radians
+ Double_t sigma_px = sigma_p * TMath::Cos(sigma_phi) * TMath::Sin(sigma_th);
+ Double_t sigma_py = sigma_p * TMath::Sin(sigma_phi) * TMath::Sin(sigma_th);
+ Double_t sigma_pz = sigma_p * TMath::Cos(sigma_th);
+ Double_t sigma_E = TMath::Sqrt(sigma_p*sigma_p + sigma_mass*sigma_mass);
+
+
+ TVector3 sigma_pvec(sigma_px, sigma_py, sigma_pz);
+
+ double cross_angle = -25./1000.; // in Rad
+ TVector3 pbeam_dir(TMath::Sin(cross_angle), 0, TMath::Cos(cross_angle)); //proton beam direction
+ sigma_pvec.RotateY(cross_angle); // Theta is returned positive, beam in negative X
+
+ // type 2 is state that will decay
+ GenParticlePtr p_sigma = std::make_shared<GenParticle>(
+ FourVector(sigma_pvec.X(), sigma_pvec.Y(), sigma_pvec.Z(), sigma_E), sigma_pdgID, 2 );
+ // Generating sigma particle, will be generated at origin
+ // Must have input electron + proton for vertex
+ GenVertexPtr sigma_initial_vertex = std::make_shared<GenVertex>();
+ sigma_initial_vertex->add_particle_in(p1);
+ sigma_initial_vertex->add_particle_in(p2);
+ sigma_initial_vertex->add_particle_out(p_sigma);
+ evt.add_vertex(sigma_initial_vertex);
+
+ // Generate lambda + gamma in sigma rest frame
+ TLorentzVector lambda_rest, gamma_rest;
+
+ // Generating uniformly along a sphere
+ double cost_lambda_rest = r1->Uniform(-1,1);
+ double th_lambda_rest = TMath::ACos(cost_lambda_rest);
+ double sint_lambda_rest = TMath::Sin(th_lambda_rest);
+
+ double phi_lambda_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
+ double cosp_lambda_rest = TMath::Cos(phi_lambda_rest);
+ double sinp_lambda_rest = TMath::Sin(phi_lambda_rest);
+
+ // Calculate energy of each particle in the sigma rest frame
+ // See problem 3.19 in Introduction to Elementary Particles, 2nd edition by D. Griffiths
+ double E_lambda_rest = (-TMath::Power(photon_mass, 2.) + TMath::Power(sigma_mass, 2.) + TMath::Power(lambda_mass, 2.) ) / (2. * sigma_mass) ;
+ double E_gamma_rest = (-TMath::Power(lambda_mass, 2.) + TMath::Power(sigma_mass, 2.) + TMath::Power(photon_mass, 2.) ) / (2. * sigma_mass) ;
+
+ // Both particles will have the same momentum, so just use lambda variables
+ double momentum_rest = TMath::Sqrt( E_lambda_rest*E_lambda_rest - lambda_mass*lambda_mass );
+
+ lambda_rest.SetE(E_lambda_rest);
+ lambda_rest.SetPx( momentum_rest * sint_lambda_rest * cosp_lambda_rest );
+ lambda_rest.SetPy( momentum_rest * sint_lambda_rest * sinp_lambda_rest );
+ lambda_rest.SetPz( momentum_rest * cost_lambda_rest );
+
+ gamma_rest.SetE(E_gamma_rest);
+ gamma_rest.SetPx( -lambda_rest.Px() );
+ gamma_rest.SetPy( -lambda_rest.Py() );
+ gamma_rest.SetPz( -lambda_rest.Pz() );
+
+ // Boost lambda & pion to lab frame
+ TLorentzVector sigma_lab(sigma_pvec.X(), sigma_pvec.Y(), sigma_pvec.Z(), sigma_E);
+ TVector3 sigma_boost = sigma_lab.BoostVector();
+ TLorentzVector lambda_lab, gamma_lab;
+ lambda_lab = lambda_rest;
+ lambda_lab.Boost(sigma_boost);
+ gamma_lab = gamma_rest;
+ gamma_lab.Boost(sigma_boost);
+
+ // Calculating position for sigma decay
+ TVector3 sigma_unit = sigma_lab.Vect().Unit();
+ double sigma_decay_length = GetDecayLength(r1, sigma_lifetime, sigma_mass, sigma_lab.P());
+ TVector3 sigma_decay_position = sigma_unit * sigma_decay_length;
+ double sigma_decay_time = sigma_decay_length / sigma_lab.Beta() ; // Decay time in lab frame in length units (mm)
+
+ // Generating vertex for sigma decay
+ GenParticlePtr p_lambda = std::make_shared<GenParticle>(
+ FourVector(lambda_lab.Px(), lambda_lab.Py(), lambda_lab.Pz(), lambda_lab.E()), lambda_pdgID, 2 );
+
+ GenParticlePtr p_gamma = std::make_shared<GenParticle>(
+ FourVector(gamma_lab.Px(), gamma_lab.Py(), gamma_lab.Pz(), gamma_lab.E()), photon_pdgID, 1 );
+
+ GenVertexPtr v_sigma_decay = std::make_shared<GenVertex>(FourVector(sigma_decay_position.X(), sigma_decay_position.Y(), sigma_decay_position.Z(), sigma_decay_time));
+ v_sigma_decay->add_particle_in(p_sigma);
+ v_sigma_decay->add_particle_out(p_lambda);
+ v_sigma_decay->add_particle_out(p_gamma);
+
+ evt.add_vertex(v_sigma_decay);
+
+ // Generate neutron + pi0 in lambda rest frame
+ TLorentzVector neutron_rest, pi0_rest;
+
+ // Generating uniformly along a sphere
+ double cost_neutron_rest = r1->Uniform(-1,1);
+ double th_neutron_rest = TMath::ACos(cost_neutron_rest);
+ double sint_neutron_rest = TMath::Sin(th_neutron_rest);
+
+ double phi_neutron_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
+ double cosp_neutron_rest = TMath::Cos(phi_neutron_rest);
+ double sinp_neutron_rest = TMath::Sin(phi_neutron_rest);
+
+ // Calculate energy of each particle in the lambda rest frame
+ // See problem 3.19 in Introduction to Elementary Particles, 2nd edition by D. Griffiths
+ double E_neutron_rest = (-TMath::Power(pi0_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(neutron_mass, 2.) ) / (2. * lambda_mass) ;
+ double E_pi0_rest = (-TMath::Power(neutron_mass, 2.) + TMath::Power(lambda_mass, 2.) + TMath::Power(pi0_mass, 2.) ) / (2. * lambda_mass) ;
+
+ // Both particles will have the same momentum, so just use neutron variables
+ momentum_rest = TMath::Sqrt( E_neutron_rest*E_neutron_rest - neutron_mass*neutron_mass );
+
+ neutron_rest.SetE(E_neutron_rest);
+ neutron_rest.SetPx( momentum_rest * sint_neutron_rest * cosp_neutron_rest );
+ neutron_rest.SetPy( momentum_rest * sint_neutron_rest * sinp_neutron_rest );
+ neutron_rest.SetPz( momentum_rest * cost_neutron_rest );
+
+ pi0_rest.SetE(E_pi0_rest);
+ pi0_rest.SetPx( -neutron_rest.Px() );
+ pi0_rest.SetPy( -neutron_rest.Py() );
+ pi0_rest.SetPz( -neutron_rest.Pz() );
+
+ // Boost neutron & pion to lab frame
+ TVector3 lambda_boost = lambda_lab.BoostVector();
+ TLorentzVector neutron_lab, pi0_lab;
+ neutron_lab = neutron_rest;
+ neutron_lab.Boost(lambda_boost);
+ pi0_lab = pi0_rest;
+ pi0_lab.Boost(lambda_boost);
+
+ // Calculating position for lambda decay
+ TVector3 lambda_unit = lambda_lab.Vect().Unit();
+ double lambda_decay_length = GetDecayLength(r1, lambda_lifetime, lambda_mass, lambda_lab.P());
+ TVector3 lambda_decay_position = lambda_unit * lambda_decay_length;
+ double lambda_decay_time = lambda_decay_length / lambda_lab.Beta() ; // Decay time in lab frame in length units (mm)
+
+ // Generating vertex for lambda decay
+ GenParticlePtr p_neutron = std::make_shared<GenParticle>(
+ FourVector(neutron_lab.Px(), neutron_lab.Py(), neutron_lab.Pz(), neutron_lab.E()), neutron_pdgID, 1 );
+
+ GenParticlePtr p_pi0 = std::make_shared<GenParticle>(
+ FourVector(pi0_lab.Px(), pi0_lab.Py(), pi0_lab.Pz(), pi0_lab.E()), pi0_pdgID, 2 );
+
+ GenVertexPtr v_lambda_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
+ v_lambda_decay->add_particle_in(p_lambda);
+ v_lambda_decay->add_particle_out(p_neutron);
+ v_lambda_decay->add_particle_out(p_pi0);
+
+ evt.add_vertex(v_lambda_decay);
+
+ // Generate two photons from pi0 decay
+ TLorentzVector gamma1_rest, gamma2_rest;
+
+ // Generating uniformly along a sphere
+ double cost_gamma1_rest = r1->Uniform(-1,1);
+ double th_gamma1_rest = TMath::ACos(cost_gamma1_rest);
+ double sint_gamma1_rest = TMath::Sin(th_gamma1_rest);
+
+ double phi_gamma1_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
+ double cosp_gamma1_rest = TMath::Cos(phi_gamma1_rest);
+ double sinp_gamma1_rest = TMath::Sin(phi_gamma1_rest);
+
+ // Photons are massless so they each get equal energies
+ gamma1_rest.SetE(pi0_mass/2.);
+ gamma1_rest.SetPx( (pi0_mass/2.)*sint_gamma1_rest*cosp_gamma1_rest );
+ gamma1_rest.SetPy( (pi0_mass/2.)*sint_gamma1_rest*sinp_gamma1_rest );
+ gamma1_rest.SetPz( (pi0_mass/2.)*cost_gamma1_rest );
+
+ gamma2_rest.SetE(pi0_mass/2.);
+ gamma2_rest.SetPx( -gamma1_rest.Px() );
+ gamma2_rest.SetPy( -gamma1_rest.Py() );
+ gamma2_rest.SetPz( -gamma1_rest.Pz() );
+
+ // Boost neutron & pion to lab frame
+ TVector3 pi0_boost = pi0_lab.BoostVector();
+ TLorentzVector gamma1_lab, gamma2_lab;
+ gamma1_lab = gamma1_rest;
+ gamma1_lab.Boost(pi0_boost);
+ gamma2_lab = gamma2_rest;
+ gamma2_lab.Boost(pi0_boost);
+
+ GenParticlePtr p_gamma1 = std::make_shared<GenParticle>(
+ FourVector(gamma1_lab.Px(), gamma1_lab.Py(), gamma1_lab.Pz(), gamma1_lab.E()), photon_pdgID, 1 );
+
+ GenParticlePtr p_gamma2 = std::make_shared<GenParticle>(
+ FourVector(gamma2_lab.Px(), gamma2_lab.Py(), gamma2_lab.Pz(), gamma2_lab.E()), photon_pdgID, 1 );
+
+ // Generate pi0 at same position as the lambda. Approximating pi0 decay as instantaneous
+ GenVertexPtr v_pi0_decay = std::make_shared<GenVertex>(FourVector(lambda_decay_position.X(), lambda_decay_position.Y(), lambda_decay_position.Z(), lambda_decay_time));
+ v_pi0_decay->add_particle_in(p_pi0);
+ v_pi0_decay->add_particle_out(p_gamma1);
+ v_pi0_decay->add_particle_out(p_gamma2);
+
+ //std::cout<< lambda_pvec.Angle(pbeam_dir) << " " << neutron_lab.Angle(pbeam_dir) << " " << gamma1_lab.Angle(pbeam_dir) << " " << gamma2_lab.Angle(pbeam_dir) << std::endl;
+
+ evt.add_vertex(v_pi0_decay);
+
+ if (events_parsed == 0) {
+ std::cout << "First event: " << std::endl;
+ Print::listing(evt);
+ }
+ double zdc_z=35800;
+ TVector3 extrap_gamma=sigma_decay_position+gamma_lab.Vect()*((zdc_z-pbeam_dir.Dot(sigma_decay_position))/(pbeam_dir.Dot(gamma_lab.Vect())));
+ TVector3 extrap_neutron=lambda_decay_position+neutron_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(neutron_lab.Vect())));
+ TVector3 extrap_gamma1=lambda_decay_position+gamma1_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma1_lab.Vect())));
+ TVector3 extrap_gamma2=lambda_decay_position+gamma2_lab.Vect()*((zdc_z-pbeam_dir.Dot(lambda_decay_position))/(pbeam_dir.Dot(gamma2_lab.Vect())));
+ if (extrap_neutron.Angle(pbeam_dir)<0.004 && extrap_gamma1.Angle(pbeam_dir)<0.004 && extrap_gamma2.Angle(pbeam_dir)<0.004 && extrap_gamma.Angle(pbeam_dir)<0.004 && lambda_decay_position.Dot(pbeam_dir)<zdc_z){
+ hepmc_output.write_event(evt);
+ accepted_events ++;
+ }
+ if (events_parsed % 1000 == 0) {
+ std::cout << "Event: " << events_parsed << " ("<<accepted_events<<" accepted)"<< std::endl;
+ }
+ evt.clear();
+ }
+ hepmc_output.close();
+
+ std::cout << "Events generated: " << events_parsed << " ("<<accepted_events<<" accepted)"<< std::endl;
+}
diff --git a/benchmarks/zdc_sigma/analysis/sigma_plots.py b/benchmarks/zdc_sigma/analysis/sigma_plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..6488b128ee74387df22079729a2d2c522505e1fd
--- /dev/null
+++ b/benchmarks/zdc_sigma/analysis/sigma_plots.py
@@ -0,0 +1,481 @@
+import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys
+import mplhep as hep
+hep.style.use("CMS")
+
+plt.rcParams['figure.facecolor']='white'
+plt.rcParams['savefig.facecolor']='white'
+plt.rcParams['savefig.bbox']='tight'
+
+plt.rcParams["figure.figsize"] = (7, 7)
+
+outdir=sys.argv[1]+"/"
+config=outdir.split("/")[1]
+try:
+ import os
+ os.mkdir(outdir[:-1])
+except:
+ pass
+import uproot as ur
+arrays_sim={}
+momenta=100, 125, 150, 175,200,225,250,275
+for p in momenta:
+ filename=f'sim_output/zdc_sigma/{config}_rec_sigma_dec_{p}GeV.edm4eic.root'
+ print("opening file", filename)
+ events = ur.open(filename+':events')
+ arrays_sim[p] = events.arrays()[:-1] #remove last event, which for some reason is blank
+ import gc
+ gc.collect()
+ print("read", filename)
+
+def gauss(x, A,mu, sigma):
+ return A * np.exp(-(x-mu)**2/(2*sigma**2))
+
+#keep track of the number of clusters per event
+nclusters={}
+
+for p in momenta:
+ plt.figure()
+ nclusters[p]=[]
+ for i in range(len(arrays_sim[p])):
+ nclusters[p].append(len(arrays_sim[p]["HcalFarForwardZDCClusters.position.x"][i]))
+ nclusters[p]=np.array(nclusters[p])
+ plt.hist(nclusters[p],bins=20, range=(0,20))
+ plt.xlabel("number of clusters")
+ plt.yscale('log')
+ plt.title(f"$p_\Sigma={p}$ GeV")
+ plt.ylim(1)
+ plt.savefig(outdir+f"nclust_{p}GeV_recon.pdf")
+ print("saved file ", outdir+f"nclust_{p}GeV_recon.pdf")
+
+
+
+pt_truth={}
+theta_truth={}
+
+for p in momenta:
+ #get the truth value of theta* and pt*
+ px=arrays_sim[p]["MCParticles.momentum.x"][:,2]
+ py=arrays_sim[p]["MCParticles.momentum.y"][:,2]
+ pz=arrays_sim[p]["MCParticles.momentum.z"][:,2]
+ tilt=-0.025
+ pt_truth[p]=np.hypot(px*np.cos(tilt)-pz*np.sin(tilt), py)
+ theta_truth[p]=np.arctan2(pt_truth[p],pz*np.cos(tilt)+px*np.sin(tilt))
+
+#create an array with the same shape as the cluster-level arrays
+is_neutron_cand={}
+for p in momenta:
+ is_neutron_cand[p]=(0*arrays_sim[p][f"HcalFarForwardZDCClusters.energy"]).to_list()
+
+ #largest_eigenvalues
+ for i in range(len(arrays_sim[p])):
+ pars=arrays_sim[p]["_HcalFarForwardZDCClusters_shapeParameters"][i]
+ index_of_max=-1
+ max_val=0
+ eigs=[]
+ #Must make sure this doesn't get messed up if someone changes the number of shape parameters in EICrecon.
+ nClust=nclusters[p][i]
+ nShapePars=len(pars)//nClust
+ for j in range(nClust):
+ largest_eigenvalue=max(pars[nShapePars*j+4:nShapePars*j+7])
+ eigs.append(largest_eigenvalue)
+ if(largest_eigenvalue>max_val):
+ max_val=largest_eigenvalue
+ index_of_max=j
+ if index_of_max >=0:
+ is_neutron_cand[p][i][index_of_max]=1
+ eigs.sort()
+
+ is_neutron_cand[p]=ak.Array(is_neutron_cand[p])
+
+import ROOT
+
+lambda_mass=1.115683
+pi0_mass=0.1349768
+pt_recon_corr={}
+theta_recon_corr={}
+mass_recon_corr={}
+mass_lambda_recon_corr={}
+mass_pi0_recon_corr={}
+pi0_converged={}
+zvtx_recon={}
+
+#run this event-by-event:
+maxZ=35800
+for p in momenta:
+ pt_recon_corr[p]=[]
+ theta_recon_corr[p]=[]
+ mass_recon_corr[p]=[]
+ mass_lambda_recon_corr[p]=[]
+ mass_pi0_recon_corr[p]=[]
+ zvtx_recon[p]=[]
+ for evt in range(len(arrays_sim[p])):
+ if nclusters[p][evt]!=4:
+ nan=-1
+ pt_recon_corr[p].append(nan)
+ theta_recon_corr[p].append(nan)
+ mass_recon_corr[p].append(nan)
+ mass_lambda_recon_corr[p].append(nan)
+ mass_pi0_recon_corr[p].append(nan)
+ zvtx_recon[p].append(nan)
+ continue
+ xc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.x"][evt]
+ yc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.y"][evt]
+ zc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.z"][evt]
+ E=arrays_sim[p][f"HcalFarForwardZDCClusters.energy"][evt]
+
+ #apply correction to the neutron candidates only
+ A,B,C=-0.0756, -1.91, 2.30
+ neutron_corr=(1-is_neutron_cand[p][evt])+is_neutron_cand[p][evt]/(1+A+B/np.sqrt(E)+C/E)
+ E=E*neutron_corr
+
+ pabs=np.sqrt(E**2-is_neutron_cand[p][evt]*.9406**2)
+ tilt=-0.025
+ xcp=xc*np.cos(tilt)-zc*np.sin(tilt)
+ ycp=yc
+ zcp=zc*np.cos(tilt)+xc*np.sin(tilt)
+
+ #search for the combination of photons that would give the best lambda mass
+ pt_best=-999
+ theta_best=-999
+ mass_lambda_best=-999
+ mass_sigma_best=-999
+ mass_pi0_best=-999
+ zvtx_best=-999
+ for hypothesis in range(4):
+ if is_neutron_cand[p][evt][hypothesis]:
+ continue
+
+ xvtx=0
+ yvtx=0
+ zvtx=0
+ #find the vertex position that reconstructs the pi0 mass
+ for iteration in range(20):
+ tot=ROOT.TLorentzVector(0,0,0,0)
+ Lambda=ROOT.TLorentzVector(0,0,0,0)
+ pi0=ROOT.TLorentzVector(0,0,0,0)
+
+ for i in range(4):
+
+ if i!=hypothesis:
+ ux=xcp[i]-xvtx
+ uy=ycp[i]-yvtx
+ uz=zcp[i]-zvtx
+ else:
+ ux=xcp[i]
+ uy=ycp[i]
+ uz=zcp[i]
+ u=np.sqrt(ux**2+uy**2+uz**2)
+ ux/=u
+ uy/=u
+ uz/=u
+
+ P=ROOT.TLorentzVector(pabs[i]*ux, pabs[i]*uy, pabs[i]*uz, E[i])
+ tot+=P
+ if not is_neutron_cand[p][evt][i] and i!=hypothesis:
+ pi0+=P
+ if i!=hypothesis:
+ Lambda+=P
+ alpha=1
+ if iteration==0:
+ zeta=1/2
+ zvtx=maxZ*np.power(zeta,alpha)
+ xvtx=Lambda.X()/Lambda.Z()*zvtx
+ yvtx=Lambda.Y()/Lambda.Z()*zvtx
+ else :
+ s=2*(pi0.M()<pi0_mass)-1
+ zeta=np.power(zvtx/maxZ, 1/alpha)
+ zeta=zeta+s*1/2**(1+iteration)
+ zvtx=maxZ*np.power(zeta,alpha)
+ xvtx=Lambda.X()/Lambda.Z()*zvtx
+ yvtx=Lambda.Y()/Lambda.Z()*zvtx
+
+ if abs(Lambda.M()-lambda_mass)< abs(mass_lambda_best-lambda_mass):
+ pt_best=tot.Pt()
+ theta_best=tot.Theta()
+ mass_lambda_best=Lambda.M()
+ mass_sigma_best=tot.M()
+ mass_pi0_best=pi0.M()
+ zvtx_best=zvtx
+
+ pt_recon_corr[p].append(pt_best)
+ theta_recon_corr[p].append(theta_best)
+ mass_recon_corr[p].append(mass_sigma_best)
+ mass_lambda_recon_corr[p].append(mass_lambda_best)
+ mass_pi0_recon_corr[p].append(mass_pi0_best)
+ zvtx_recon[p].append(zvtx_best)
+ pt_recon_corr[p]=ak.Array(pt_recon_corr[p])
+ theta_recon_corr[p]=ak.Array(theta_recon_corr[p])
+ mass_recon_corr[p]=ak.Array(mass_recon_corr[p])
+ mass_lambda_recon_corr[p]=ak.Array(mass_lambda_recon_corr[p])
+ mass_pi0_recon_corr[p]=ak.Array(mass_pi0_recon_corr[p])
+ zvtx_recon[p]=ak.Array(zvtx_recon[p])
+
+#now make plots
+
+#reconstructed vertex position plot
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+plt.sca(axs[0])
+plt.title(f"$E_{{\\Sigma}}=100-275$ GeV")
+x=[]
+y=[]
+for p in momenta:
+ accept=(nclusters[p]==4)# &(pi0_converged[p])
+ x+=list(theta_truth[p][accept]*1000)
+ y+=list(theta_recon_corr[p][accept]*1000)
+#print(x)
+plt.scatter(x,y)
+plt.xlabel("$\\theta^{*\\rm truth}_{\\Sigma}$ [mrad]")
+plt.ylabel("$\\theta^{*\\rm recon}_{\\Sigma}$ [mrad]")
+plt.xlim(0,3.2)
+plt.ylim(0,3.2)
+
+plt.sca(axs[1])
+plt.title(f"$E_{{\\Sigma}}=100-275$ GeV")
+y,x,_=plt.hist(y-np.array(x), bins=50, range=(-1,1))
+bc=(x[1:]+x[:-1])/2
+
+from scipy.optimize import curve_fit
+slc=abs(bc)<0.6
+fnc=gauss
+p0=[100, 0, 0.5]
+coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
+ sigma=np.sqrt(y[slc])+(y[slc]==0))
+x=np.linspace(-1, 1)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+plt.xlabel("$\\theta^{*\\rm recon}_{\\Sigma}-\\theta^{*\\rm truth}_{\\Sigma}$ [mrad]")
+plt.ylabel("events")
+
+plt.sca(axs[2])
+sigmas=[]
+dsigmas=[]
+xvals=[]
+for p in momenta:
+
+ accept=(nclusters[p]==4)
+ y,x=np.histogram((theta_recon_corr[p]-theta_truth[p])[accept]*1000, bins=100, range=(-0.5,0.5))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc)<0.3
+ fnc=gauss
+ p0=(100, 0, 0.15)
+ #print(bc[slc],y[slc])
+ sigma=np.sqrt(y[slc])+(y[slc]==0)
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+ sigmas.append(np.abs(coeff[2]))
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print(f"fit failed for p={p}")
+print(xvals)
+print(sigmas)
+plt.ylim(0, 0.3)
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+x=np.linspace(100, 275, 100)
+plt.plot(x, 3/np.sqrt(x), color='tab:orange')
+plt.text(170, .23, "YR requirement:\n 3 mrad/$\\sqrt{E}$")
+plt.xlabel("$E_{\\Sigma}$ [GeV]")
+plt.ylabel("$\\sigma[\\theta^*_{\\Sigma}]$ [mrad]")
+plt.tight_layout()
+plt.savefig(outdir+"thetastar_recon.pdf")
+
+#reconstructed vertex position plot
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+plt.sca(axs[0])
+plt.title(f"$E_{{\\Sigma}}=100-275$ GeV")
+x=[]
+y=[]
+for p in momenta:
+ accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+ tilt=-0.025
+ x+=list(arrays_sim[p]['MCParticles.vertex.z'][:,5][accept]*np.cos(tilt)/1000
+ +np.sin(tilt)*arrays_sim[p]['MCParticles.vertex.z'][:,5][accept]/1000)
+ y+=list(zvtx_recon[p][accept]/1000)
+plt.scatter(x,y)
+#print(x,y)
+plt.xlabel("$z^{\\rm truth}_{\\rm vtx}$ [m]")
+plt.ylabel("$z^{\\rm recon}_{\\rm vtx}$ [m]")
+plt.xlim(0,40)
+plt.ylim(0,40)
+
+plt.sca(axs[1])
+plt.title(f"$E_{{\\Sigma}}=100-275$ GeV")
+y,x,_=plt.hist(y-np.array(x), bins=50, range=(-10,10))
+bc=(x[1:]+x[:-1])/2
+
+from scipy.optimize import curve_fit
+slc=abs(bc)<5
+fnc=gauss
+p0=[100, 0, 1]
+coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
+ sigma=np.sqrt(y[slc])+(y[slc]==0))
+x=np.linspace(-5, 5)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+plt.xlabel("$z^{*\\rm recon}_{\\rm vtx}-z^{*\\rm truth}_{\\rm vtx}$ [m]")
+plt.ylabel("events")
+
+plt.sca(axs[2])
+sigmas=[]
+dsigmas=[]
+xvals=[]
+for p in momenta:
+
+ accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+ a=list((zvtx_recon[p]-arrays_sim[p]['MCParticles.vertex.z'][:,5])[accept]/1000)
+ y,x=np.histogram(a, bins=100, range=(-10,10))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc)<5
+ fnc=gauss
+ p0=(100, 0, 1)
+ #print(bc[slc],y[slc])
+ sigma=np.sqrt(y[slc])+(y[slc]==0)
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+ sigmas.append(abs(coeff[2]))
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print(f"fit failed for p={p}")
+plt.ylim(0, 3)
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+x=np.linspace(100, 275, 100)
+
+avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
+plt.axhline(avg, color='tab:orange')
+plt.text(150, 1.25,f"$\\sigma\\approx${avg:.1f} m")
+
+plt.xlabel("$E_{\\Sigma}$ [GeV]")
+plt.ylabel("$\\sigma[z_{\\rm vtx}]$ [m]")
+plt.tight_layout()
+plt.savefig(outdir+"zvtx_recon.pdf")
+
+#lambda mass reconstruction
+fig,axs=plt.subplots(1,2, figsize=(16, 8))
+plt.sca(axs[0])
+lambda_mass=1.115683
+vals=[]
+for p in momenta:
+ accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+ vals+=list(mass_lambda_recon_corr[p][accept])
+
+y,x,_= plt.hist(vals, bins=100, range=(0.9, 1.3))
+bc=(x[1:]+x[:-1])/2
+plt.axvline(lambda_mass, ls='--', color='tab:green', lw=3)
+plt.text(lambda_mass+.01, np.max(y)*1.05, "PDG mass", color='tab:green')
+plt.xlabel("$m_{\\Lambda}^{\\rm recon}$ [GeV]")
+plt.ylim(0, np.max(y)*1.2)
+plt.xlim(0.9, 1.3)
+
+from scipy.optimize import curve_fit
+slc=abs(bc-lambda_mass)<0.05
+fnc=gauss
+p0=[100, lambda_mass, 0.03]
+coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
+ sigma=np.sqrt(y[slc])+(y[slc]==0))
+x=np.linspace(0.8, 1.3, 200)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+print(coeff[2], np.sqrt(var_matrix[2][2]))
+plt.xlabel("$m^{\\rm recon}_{\\Lambda}$ [GeV]")
+plt.ylabel("events")
+plt.title(f"$E_{{\\Sigma}}=100-275$ GeV")
+
+plt.sca(axs[1])
+xvals=[]
+sigmas=[]
+dsigmas=[]
+for p in momenta:
+ accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+ y,x= np.histogram(mass_lambda_recon_corr[p][accept], bins=100, range=(0.6,1.4))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc-lambda_mass)<0.05
+ fnc=gauss
+ p0=[100, lambda_mass, 0.03]
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+ x=np.linspace(0.8, 1.3, 200)
+ sigmas.append(coeff[2])
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print("fit failed")
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
+plt.axhline(avg, color='tab:orange')
+plt.text(150, 0.01,f"$\\sigma\\approx${avg*1000:.0f} MeV")
+plt.xlabel("$E_{\\Sigma}$ [GeV]")
+plt.ylabel("$\\sigma[m_{\\Lambda}]$ [GeV]")
+plt.ylim(0, 0.02)
+plt.tight_layout()
+plt.savefig(outdir+"lambda_mass_rec_from_sigma_decay.pdf")
+
+#sigma mass reconstruction
+p=100
+fig,axs=plt.subplots(1,2, figsize=(16, 8))
+plt.sca(axs[0])
+sigma_mass=1.192
+vals=[]
+for p in momenta:
+ accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+ vals+=list(mass_recon_corr[p][accept])
+
+y,x,_= plt.hist(vals, bins=100, range=(1.0, 1.4))
+bc=(x[1:]+x[:-1])/2
+plt.axvline(sigma_mass, ls='--', color='tab:green', lw=3)
+plt.text(sigma_mass+.01, np.max(y)*1.05, "PDG mass", color='tab:green')
+plt.xlabel("$m_{\\Sigma}^{\\rm recon}$ [GeV]")
+plt.ylim(0, np.max(y)*1.2)
+plt.xlim(1.0, 1.45)
+
+from scipy.optimize import curve_fit
+slc=abs(bc-sigma_mass)<0.02
+fnc=gauss
+p0=[100, sigma_mass, 0.03]
+coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
+ sigma=np.sqrt(y[slc])+(y[slc]==0))
+x=np.linspace(0.8, 1.3, 200)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+print(coeff[2], np.sqrt(var_matrix[2][2]))
+plt.xlabel("$m^{\\rm recon}_{\\Sigma}$ [GeV]")
+plt.ylabel("events")
+plt.title(f"$E_{{\\Sigma}}=100-275$ GeV")
+
+plt.sca(axs[1])
+xvals=[]
+sigmas=[]
+dsigmas=[]
+for p in momenta:
+ accept=(nclusters[p]==4)&(abs(mass_pi0_recon_corr[p]-pi0_mass)<.01)
+ y,x= np.histogram(mass_recon_corr[p][accept], bins=100, range=(1.0,1.4))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc-sigma_mass)<0.02
+ fnc=gauss
+ p0=[100, sigma_mass, 0.03]
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+ sigmas.append(abs(coeff[2]))
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print("fit failed")
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
+plt.axhline(avg, color='tab:orange')
+plt.text(150, 0.01,f"$\\sigma\\approx${avg*1000:.0f} MeV")
+plt.xlabel("$E_{\\Sigma}$ [GeV]")
+plt.ylabel("$\\sigma[m_{\\Sigma}]$ [GeV]")
+plt.ylim(0, 0.1)
+plt.tight_layout()
+plt.savefig(outdir+"sigma_mass_rec.pdf")
diff --git a/benchmarks/zdc_sigma/config.yml b/benchmarks/zdc_sigma/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..fc81428a64275692721d236bec268f5cb087e534
--- /dev/null
+++ b/benchmarks/zdc_sigma/config.yml
@@ -0,0 +1,36 @@
+sim:zdc_sigma:
+ stage: simulate
+ extends: .det_benchmark
+ parallel:
+ matrix:
+ - P: 100
+ - P: 125
+ - P: 150
+ - P: 175
+ - P: 200
+ - P: 225
+ - P: 250
+ - P: 275
+ timeout: 1 hours
+ script:
+ - snakemake --cores 1 sim_output/zdc_sigma/epic_zdc_sipm_on_tile_only_rec_sigma_dec_${P}GeV.edm4eic.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:zdc_sigma:
+ stage: benchmarks
+ extends: .det_benchmark
+ needs: ["sim:zdc_sigma"]
+ script:
+ - snakemake --cores 1 results/epic_zdc_sipm_on_tile_only/zdc_sigma
+
+collect_results:zdc_sigma:
+ stage: collect
+ needs: ["bench:zdc_sigma"]
+ script:
+ - ls -al
+ - mv results{,_save}/ # move results directory out of the way to preserve it
+ - snakemake --cores 1 --delete-all-output results/epic_zdc_sipm_on_tile_only/zdc_sigma
+ - mv results{_save,}/