From c72c399b3bba78d55114f6cafd18a7ebe79cd47b Mon Sep 17 00:00:00 2001
From: Dmitry Kalinkin <dmitry.kalinkin@gmail.com>
Date: Wed, 2 Oct 2024 08:39:47 -0400
Subject: [PATCH] parallelize zdc and insert simulations (#74)
---
benchmarks/insert_muon/Snakefile | 20 ++-
benchmarks/insert_muon/analysis/muon_plots.py | 169 +++++++++---------
benchmarks/insert_muon/config.yml | 3 +-
benchmarks/insert_neutron/Snakefile | 20 ++-
.../insert_neutron/analysis/neutron_plots.py | 6 +-
benchmarks/insert_neutron/config.yml | 3 +-
benchmarks/zdc_lambda/Snakefile | 30 ++--
.../zdc_lambda/analysis/gen_lambda_decay.cxx | 15 +-
.../zdc_lambda/analysis/lambda_plots.py | 11 +-
benchmarks/zdc_lambda/config.yml | 3 +-
benchmarks/zdc_photon/Snakefile | 34 ++--
.../zdc_photon/analysis/zdc_photon_plots.py | 13 +-
benchmarks/zdc_photon/config.yml | 3 +-
benchmarks/zdc_pi0/Snakefile | 20 ++-
benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py | 11 +-
benchmarks/zdc_pi0/config.yml | 3 +-
benchmarks/zdc_sigma/Snakefile | 28 ++-
.../zdc_sigma/analysis/gen_sigma_decay.cxx | 15 +-
benchmarks/zdc_sigma/analysis/sigma_plots.py | 11 +-
benchmarks/zdc_sigma/config.yml | 3 +-
20 files changed, 199 insertions(+), 222 deletions(-)
diff --git a/benchmarks/insert_muon/Snakefile b/benchmarks/insert_muon/Snakefile
index b36026b0..a9a32168 100644
--- a/benchmarks/insert_muon/Snakefile
+++ b/benchmarks/insert_muon/Snakefile
@@ -14,17 +14,19 @@ root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "mu-", {
rule insert_muon_simulate:
input:
- GEN_FILE="sim_output/insert_muon/mu-_{P}GeV.hepmc"
+ GEN_FILE="sim_output/insert_muon/mu-_{P}GeV.hepmc",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
params:
PHYSICS_LIST="FTFP_BERT"
output:
- SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV_{INDEX}.edm4hep.root"
shell:
"""
-NEVENTS_SIM=5000
+NEVENTS_SIM=1000
# Running simulation
npsim \
--compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \
--numberOfEvents $NEVENTS_SIM \
--physicsList {params.PHYSICS_LIST} \
--inputFiles {input.GEN_FILE} \
@@ -33,20 +35,22 @@ npsim \
rule insert_muon_recon:
input:
- SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV_{INDEX}.edm4hep.root"
output:
- REC_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_rec_mu-_{P}GeV.edm4hep.root"
+ REC_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_rec_mu-_{P}GeV_{INDEX}.edm4hep.root"
shell:
"""
-NEVENTS_REC=5000
+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_muon_analysis:
input:
- expand("sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV.edm4hep.root",
+ expand("sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV_{INDEX}.edm4hep.root",
P=[50],
- DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
+ INDEX=range(5),
+ ),
script="benchmarks/insert_muon/analysis/muon_plots.py",
output:
results_dir=directory("results/{DETECTOR_CONFIG}/insert_muon"),
diff --git a/benchmarks/insert_muon/analysis/muon_plots.py b/benchmarks/insert_muon/analysis/muon_plots.py
index 4c45aac3..5c81b52d 100644
--- a/benchmarks/insert_muon/analysis/muon_plots.py
+++ b/benchmarks/insert_muon/analysis/muon_plots.py
@@ -1,86 +1,83 @@
-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)
-
-config=sys.argv[1].split("/")[1] #results/{config}/{benchmark_name}
-outdir=sys.argv[1]+"/"
-try:
- import os
- os.mkdir(outdir[:-1])
-except:
- pass
-
-def Landau(x, normalization,location,stdev):
- #print(type(x))
- u=(x-location)*3.591/stdev/2.355
- renormalization = 1.64872*normalization
- return renormalization * np.exp(-u/2 - np.exp(-u)/2)
-
-import uproot as ur
-arrays_sim={}
-momenta=50,
-for p in momenta:
- filename=f'sim_output/insert_muon/{config}_sim_mu-_{p}GeV.edm4hep.root'
- print("opening file", filename)
- events = ur.open(filename+':events')
- arrays_sim[p] = events.arrays()
- import gc
- gc.collect()
- print("read", filename)
-
-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['eta_truth']=1/2*np.log((p+pzp)/(p-pzp))
- array['phi_truth']=np.arctan2(pyp,pxp)
-
-for p in 50,:
- E=arrays_sim[p]["HcalEndcapPInsertHits.energy"]
- y, x,_=plt.hist(1e3*ak.flatten(E),bins=100, range=(0, 1.2), histtype='step')
- bc=(x[1:]+x[:-1])/2
- from scipy.optimize import curve_fit
- slc=abs(bc-.48)<.15
- fnc=Landau
- p0=[100, .5, .05]
- #print(list(y), list(x))
- coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
- sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
- print(coeff)
- xx=np.linspace(0,.7, 100)
- MIP=coeff[1]/1000
- plt.plot(xx, fnc(xx,*coeff), label=f'Landau fit:\nMIP={coeff[1]*1e3:.0f}$\\pm${1e3*np.sqrt(var_matrix[1][1]):.0f} keV')
- plt.xlabel("hit energy [MeV]")
- plt.ylabel("hits")
- plt.title(f"$E_{{\\mu^-}}=${p} GeV")
- plt.legend(fontsize=20)
- plt.savefig(outdir+"/MIP.pdf")
-
- plt.figure(figsize=(10,7))
- array=arrays_sim[p]
- bins=30; r=((-np.pi, np.pi),(2.8, 4.2))
- selection=np.sum(array["HcalEndcapPInsertHits.energy"]>0.5*MIP,axis=-1)>0
- h1, xedges, yedges = np.histogram2d(list(array[selection]['phi_truth']),list(array[selection]['eta_truth']), bins=bins, range=r)
- h2, xedges, yedges = np.histogram2d(list(array['phi_truth']),list(array['eta_truth']), bins=bins, range=r)
-
- h = h1 / h2
- pc=plt.pcolor(xedges, yedges, h.T,linewidth=0)
- plt.xlabel("$\\phi^*$ [rad]")
- plt.ylabel("$\\eta^*$")
- cb = plt.colorbar(pc)
- cb.set_label("acceptance")
- plt.title(f"$E_{{\\mu^-}}=${p} GeV")
- plt.savefig(outdir+"/acceptance.pdf")
+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)
+
+config=sys.argv[1].split("/")[1] #results/{config}/{benchmark_name}
+outdir=sys.argv[1]+"/"
+try:
+ import os
+ os.mkdir(outdir[:-1])
+except:
+ pass
+
+def Landau(x, normalization,location,stdev):
+ #print(type(x))
+ u=(x-location)*3.591/stdev/2.355
+ renormalization = 1.64872*normalization
+ return renormalization * np.exp(-u/2 - np.exp(-u)/2)
+
+import uproot as ur
+arrays_sim={}
+momenta=50,
+for p in momenta:
+ arrays_sim[p] = ur.concatenate({
+ f'sim_output/insert_muon/{config}_sim_mu-_{p}GeV_{index}.edm4hep.root': 'events'
+ for index in range(5)
+ })
+
+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['eta_truth']=1/2*np.log((p+pzp)/(p-pzp))
+ array['phi_truth']=np.arctan2(pyp,pxp)
+
+for p in 50,:
+ E=arrays_sim[p]["HcalEndcapPInsertHits.energy"]
+ y, x,_=plt.hist(1e3*ak.flatten(E),bins=100, range=(0, 1.2), histtype='step')
+ bc=(x[1:]+x[:-1])/2
+ from scipy.optimize import curve_fit
+ slc=abs(bc-.48)<.15
+ fnc=Landau
+ p0=[100, .5, .05]
+ #print(list(y), list(x))
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+ print(coeff)
+ xx=np.linspace(0,.7, 100)
+ MIP=coeff[1]/1000
+ plt.plot(xx, fnc(xx,*coeff), label=f'Landau fit:\nMIP={coeff[1]*1e3:.0f}$\\pm${1e3*np.sqrt(var_matrix[1][1]):.0f} keV')
+ plt.xlabel("hit energy [MeV]")
+ plt.ylabel("hits")
+ plt.title(f"$E_{{\\mu^-}}=${p} GeV")
+ plt.legend(fontsize=20)
+ plt.savefig(outdir+"/MIP.pdf")
+
+ plt.figure(figsize=(10,7))
+ array=arrays_sim[p]
+ bins=30; r=((-np.pi, np.pi),(2.8, 4.2))
+ selection=np.sum(array["HcalEndcapPInsertHits.energy"]>0.5*MIP,axis=-1)>0
+ h1, xedges, yedges = np.histogram2d(list(array[selection]['phi_truth']),list(array[selection]['eta_truth']), bins=bins, range=r)
+ h2, xedges, yedges = np.histogram2d(list(array['phi_truth']),list(array['eta_truth']), bins=bins, range=r)
+
+ h = h1 / h2
+ pc=plt.pcolor(xedges, yedges, h.T,linewidth=0)
+ plt.xlabel("$\\phi^*$ [rad]")
+ plt.ylabel("$\\eta^*$")
+ cb = plt.colorbar(pc)
+ cb.set_label("acceptance")
+ plt.title(f"$E_{{\\mu^-}}=${p} GeV")
+ plt.savefig(outdir+"/acceptance.pdf")
diff --git a/benchmarks/insert_muon/config.yml b/benchmarks/insert_muon/config.yml
index 6f7681f1..252cd6e8 100644
--- a/benchmarks/insert_muon/config.yml
+++ b/benchmarks/insert_muon/config.yml
@@ -4,9 +4,8 @@ sim:insert_muon:
parallel:
matrix:
- P: 50
- timeout: 1 hours
script:
- - snakemake --cores 1 sim_output/insert_muon/epic_craterlake_sim_mu-_${P}GeV.edm4hep.root
+ - snakemake --cores 5 sim_output/insert_muon/epic_craterlake_sim_mu-_${P}GeV_{0,1,2,3,4}.edm4hep.root
retry:
max: 2
when:
diff --git a/benchmarks/insert_neutron/Snakefile b/benchmarks/insert_neutron/Snakefile
index 55e3ff9e..136c56e2 100644
--- a/benchmarks/insert_neutron/Snakefile
+++ b/benchmarks/insert_neutron/Snakefile
@@ -15,17 +15,19 @@ root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "neutron
rule insert_neutron_simulate:
input:
- GEN_FILE="sim_output/insert_neutron/neutron_{P}GeV.hepmc"
+ GEN_FILE="sim_output/insert_neutron/neutron_{P}GeV.hepmc",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
params:
PHYSICS_LIST="FTFP_BERT"
output:
- SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV_{INDEX}.edm4hep.root"
shell:
"""
-NEVENTS_SIM=1000
+NEVENTS_SIM=200
# Running simulation
npsim \
--compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \
--numberOfEvents $NEVENTS_SIM \
--physicsList {params.PHYSICS_LIST} \
--inputFiles {input.GEN_FILE} \
@@ -34,20 +36,22 @@ npsim \
rule insert_neutron_recon:
input:
- SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV_{INDEX}.edm4hep.root"
output:
- REC_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV.edm4eic.root"
+ REC_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV_{INDEX}.edm4eic.root"
shell:
"""
-NEVENTS_REC=1000
+NEVENTS_REC=200
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",
+ expand("sim_output/insert_neutron/{DETECTOR_CONFIG}_rec_neutron_{P}GeV_{INDEX}.edm4eic.root",
P=[20, 30, 40, 50, 60, 70, 80],
- DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
+ INDEX=range(5),
+ ),
script="benchmarks/insert_neutron/analysis/neutron_plots.py",
output:
results_dir=directory("results/{DETECTOR_CONFIG}/insert_neutron"),
diff --git a/benchmarks/insert_neutron/analysis/neutron_plots.py b/benchmarks/insert_neutron/analysis/neutron_plots.py
index 9314a974..b9d2002a 100644
--- a/benchmarks/insert_neutron/analysis/neutron_plots.py
+++ b/benchmarks/insert_neutron/analysis/neutron_plots.py
@@ -18,8 +18,10 @@ except:
#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()
+ arrays_sim[p] = ur.concatenate({
+ f'sim_output/insert_neutron/{config}_rec_neutron_{p}GeV_{index}.edm4eic.root': 'events'
+ for index in range(5)
+ })
def gauss(x, A,mu, sigma):
return A * np.exp(-(x-mu)**2/(2*sigma**2))
diff --git a/benchmarks/insert_neutron/config.yml b/benchmarks/insert_neutron/config.yml
index aac0a732..4c723e68 100644
--- a/benchmarks/insert_neutron/config.yml
+++ b/benchmarks/insert_neutron/config.yml
@@ -10,9 +10,8 @@ sim:insert_neutron:
- 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
+ - snakemake --cores 5 sim_output/insert_neutron/epic_craterlake_rec_neutron_${P}GeV_{0,1,2,3,4}.edm4eic.root
retry:
max: 2
when:
diff --git a/benchmarks/zdc_lambda/Snakefile b/benchmarks/zdc_lambda/Snakefile
index 21563330..9ffd1701 100644
--- a/benchmarks/zdc_lambda/Snakefile
+++ b/benchmarks/zdc_lambda/Snakefile
@@ -2,7 +2,7 @@ rule zdc_lambda_generate:
input:
script="benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx",
params:
- NEVENTS_GEN=100000,
+ NEVENTS_GEN=1000,
output:
GEN_FILE="sim_output/zdc_lambda/lambda_decay_{P}GeV.hepmc"
shell:
@@ -12,21 +12,19 @@ root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildca
rule zdc_lambda_simulate:
input:
- GEN_FILE="sim_output/zdc_lambda/lambda_decay_{P}GeV.hepmc"
+ GEN_FILE="sim_output/zdc_lambda/lambda_decay_{P}GeV.hepmc",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
params:
PHYSICS_LIST="FTFP_BERT"
output:
- SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV_{INDEX}.edm4hep.root"
shell:
"""
-if [[ {wildcards.P} -gt 220 ]]; then
- NEVENTS_SIM=500
-else
- NEVENTS_SIM=1000
-fi
+NEVENTS_SIM=200
# Running simulation
npsim \
--compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \
--numberOfEvents $NEVENTS_SIM \
--physicsList {params.PHYSICS_LIST} \
--inputFiles {input.GEN_FILE} \
@@ -35,24 +33,22 @@ npsim \
rule zdc_lambda_recon:
input:
- SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_sim_lambda_dec_{P}GeV_{INDEX}.edm4hep.root"
output:
- REC_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV.edm4eic.root"
+ REC_FILE="sim_output/zdc_lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV_{INDEX}.edm4eic.root"
shell:
"""
-if [[ {wildcards.P} -gt 220 ]]; then
- NEVENTS_REC=500
-else
- NEVENTS_REC=1000
-fi
+NEVENTS_REC=200
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",
+ expand("sim_output/zdc_lambda/{DETECTOR_CONFIG}_rec_lambda_dec_{P}GeV_{INDEX}.edm4eic.root",
P=[100, 125, 150,175, 200, 225, 250, 275],
- DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
+ INDEX=range(5),
+ ),
script="benchmarks/zdc_lambda/analysis/lambda_plots.py",
output:
results_dir=directory("results/{DETECTOR_CONFIG}/zdc_lambda"),
diff --git a/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx b/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx
index 567eda5b..1b399dda 100644
--- a/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx
+++ b/benchmarks/zdc_lambda/analysis/gen_lambda_decay.cxx
@@ -43,7 +43,6 @@ void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname =
//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
@@ -71,10 +70,11 @@ void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname =
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++) {
+ int accepted_events = 0;
+ while (accepted_events < n_events) {
//Set the event number
- evt.set_event_number(events_parsed);
+ evt.set_event_number(accepted_events);
// FourVector(px,py,pz,e,pdgid,status)
// type 4 is beam
@@ -218,7 +218,7 @@ void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname =
evt.add_vertex(v_pi0_decay);
- if (events_parsed == 0) {
+ if (accepted_events == 0) {
std::cout << "First event: " << std::endl;
Print::listing(evt);
}
@@ -227,13 +227,14 @@ void gen_lambda_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname =
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)
+ accepted_events++;
hepmc_output.write_event(evt);
- if (events_parsed % 1000 == 0) {
- std::cout << "Event: " << events_parsed << std::endl;
+ if (accepted_events % 1000 == 0) {
+ std::cout << "Event: " << accepted_events << std::endl;
}
evt.clear();
}
hepmc_output.close();
- std::cout << "Events parsed and written: " << events_parsed << std::endl;
+ std::cout << "Events parsed and written: " << accepted_events << std::endl;
}
diff --git a/benchmarks/zdc_lambda/analysis/lambda_plots.py b/benchmarks/zdc_lambda/analysis/lambda_plots.py
index 88aff7f3..bf855d9d 100644
--- a/benchmarks/zdc_lambda/analysis/lambda_plots.py
+++ b/benchmarks/zdc_lambda/analysis/lambda_plots.py
@@ -19,13 +19,10 @@ 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)
+ arrays_sim[p] = ur.concatenate({
+ f'sim_output/zdc_lambda/{config}_rec_lambda_dec_{p}GeV_{index}.edm4eic.root': 'events'
+ for index in range(5)
+ })
def gauss(x, A,mu, sigma):
return A * np.exp(-(x-mu)**2/(2*sigma**2))
diff --git a/benchmarks/zdc_lambda/config.yml b/benchmarks/zdc_lambda/config.yml
index 396b1924..e62b0d9c 100644
--- a/benchmarks/zdc_lambda/config.yml
+++ b/benchmarks/zdc_lambda/config.yml
@@ -11,9 +11,8 @@ sim:zdc_lambda:
- 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
+ - snakemake --cores 5 sim_output/zdc_lambda/epic_zdc_sipm_on_tile_only_rec_lambda_dec_${P}GeV_{0,1,2,3,4}.edm4eic.root
retry:
max: 2
when:
diff --git a/benchmarks/zdc_photon/Snakefile b/benchmarks/zdc_photon/Snakefile
index cfb51242..2ca2ad3d 100644
--- a/benchmarks/zdc_photon/Snakefile
+++ b/benchmarks/zdc_photon/Snakefile
@@ -8,33 +8,27 @@ rule zdc_photon_generate:
GEN_FILE="sim_output/zdc_photon/zdc_photon_{P}GeV.hepmc"
shell:
"""
-if [[ {wildcards.P} -gt 140 ]]; then
- NEVENTS_GEN=300
-else
- NEVENTS_GEN=1000
-fi
+NEVENTS_GEN=200
mkdir -p sim_output/zdc_photon
root -l -b -q '{input.script}('$NEVENTS_GEN',"{output.GEN_FILE}", "gamma", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
"""
rule zdc_photon_simulate:
input:
- GEN_FILE="sim_output/zdc_photon/zdc_photon_{P}GeV.hepmc"
+ GEN_FILE="sim_output/zdc_photon/zdc_photon_{P}GeV.hepmc",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
params:
PHYSICS_LIST="FTFP_BERT"
output:
- SIM_FILE="sim_output/zdc_photon/{DETECTOR_CONFIG}_sim_zdc_photon_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/zdc_photon/{DETECTOR_CONFIG}_sim_zdc_photon_{P}GeV_{INDEX}.edm4hep.root"
shell:
"""
# Running simulation
-if [[ {wildcards.P} -gt 140 ]]; then
- NEVENTS_SIM=300
-else
- NEVENTS_SIM=1000
-fi
+NEVENTS_SIM=200
npsim \
--compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
--physicsList {params.PHYSICS_LIST} \
+ --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \
--numberOfEvents $NEVENTS_SIM \
--inputFiles {input.GEN_FILE} \
--outputFile {output.SIM_FILE}
@@ -42,24 +36,22 @@ npsim \
rule zdc_photon_recon:
input:
- SIM_FILE="sim_output/zdc_photon/{DETECTOR_CONFIG}_sim_zdc_photon_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/zdc_photon/{DETECTOR_CONFIG}_sim_zdc_photon_{P}GeV_{INDEX}.edm4hep.root"
output:
- REC_FILE="sim_output/zdc_photon/{DETECTOR_CONFIG}_rec_zdc_photon_{P}GeV.edm4hep.root"
+ REC_FILE="sim_output/zdc_photon/{DETECTOR_CONFIG}_rec_zdc_photon_{P}GeV_{INDEX}.edm4eic.root"
shell:
"""
-if [[ {wildcards.P} -gt 140 ]]; then
- NEVENTS_REC=300
-else
- NEVENTS_REC=1000
-fi
+NEVENTS_REC=200
eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,HcalFarForwardZDCSubcellHits -Pjana:nevents=$NEVENTS_REC
"""
rule zdc_photon_analysis:
input:
- expand("sim_output/zdc_photon/{DETECTOR_CONFIG}_rec_zdc_photon_{P}GeV.edm4hep.root",
+ expand("sim_output/zdc_photon/{DETECTOR_CONFIG}_rec_zdc_photon_{P}GeV_{INDEX}.edm4eic.root",
P=[20, 30, 50, 70, 100, 150, 200, 275],
- DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
+ INDEX=range(5),
+ ),
script="benchmarks/zdc_photon/analysis/zdc_photon_plots.py",
output:
results_dir=directory("results/{DETECTOR_CONFIG}/zdc_photon"),
diff --git a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
index b06cf434..94a93101 100644
--- a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
+++ b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
@@ -23,14 +23,11 @@ import uproot as ur
arrays_sim={}
momenta=20, 30, 50, 70, 100, 150, 200, 275
for p in momenta:
- filename=f'sim_output/zdc_photon/{config}_rec_zdc_photon_{p}GeV.edm4hep.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)
-
+ arrays_sim[p] = ur.concatenate({
+ f'sim_output/zdc_photon/{config}_rec_zdc_photon_{p}GeV_{index}.edm4eic.root': 'events'
+ for index in range(5)
+ })
+
fig,axs=plt.subplots(1,3, figsize=(24, 8))
pvals=[]
resvals=[]
diff --git a/benchmarks/zdc_photon/config.yml b/benchmarks/zdc_photon/config.yml
index 946650a3..4d205eb2 100644
--- a/benchmarks/zdc_photon/config.yml
+++ b/benchmarks/zdc_photon/config.yml
@@ -11,9 +11,8 @@ sim:zdc_photon:
- P: 150
- P: 200
- P: 275
- timeout: 1 hours
script:
- - snakemake --cores 1 sim_output/zdc_photon/epic_zdc_sipm_on_tile_only_rec_zdc_photon_${P}GeV.edm4hep.root
+ - snakemake --cores 5 sim_output/zdc_photon/epic_zdc_sipm_on_tile_only_rec_zdc_photon_${P}GeV_{0,1,2,3,4}.edm4eic.root
retry:
max: 2
when:
diff --git a/benchmarks/zdc_pi0/Snakefile b/benchmarks/zdc_pi0/Snakefile
index 57fd2c90..c2a67993 100644
--- a/benchmarks/zdc_pi0/Snakefile
+++ b/benchmarks/zdc_pi0/Snakefile
@@ -13,17 +13,19 @@ root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildca
rule zdc_pi0_simulate:
input:
- GEN_FILE="sim_output/zdc_pi0/zdc_pi0_{P}GeV.hepmc"
+ GEN_FILE="sim_output/zdc_pi0/zdc_pi0_{P}GeV.hepmc",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
params:
PHYSICS_LIST="FTFP_BERT"
output:
- SIM_FILE="sim_output/zdc_pi0/{DETECTOR_CONFIG}_sim_zdc_pi0_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/zdc_pi0/{DETECTOR_CONFIG}_sim_zdc_pi0_{P}GeV_{INDEX}.edm4hep.root"
shell:
"""
-NEVENTS_SIM=1000
+NEVENTS_SIM=200
# Running simulation
npsim \
--compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \
--numberOfEvents $NEVENTS_SIM \
--physicsList {params.PHYSICS_LIST} \
--inputFiles {input.GEN_FILE} \
@@ -32,20 +34,22 @@ npsim \
rule zdc_pi0_recon:
input:
- SIM_FILE="sim_output/zdc_pi0/{DETECTOR_CONFIG}_sim_zdc_pi0_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/zdc_pi0/{DETECTOR_CONFIG}_sim_zdc_pi0_{P}GeV_{INDEX}.edm4hep.root"
output:
- REC_FILE="sim_output/zdc_pi0/{DETECTOR_CONFIG}_rec_zdc_pi0_{P}GeV.edm4hep.root"
+ REC_FILE="sim_output/zdc_pi0/{DETECTOR_CONFIG}_rec_zdc_pi0_{P}GeV_{INDEX}.edm4eic.root"
shell:
"""
-NEVENTS_REC=1000
+NEVENTS_REC=200
eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,HcalFarForwardZDCSubcellHits -Pjana:nevents=$NEVENTS_REC
"""
rule zdc_pi0_analysis:
input:
- expand("sim_output/zdc_pi0/{DETECTOR_CONFIG}_rec_zdc_pi0_{P}GeV.edm4hep.root",
+ expand("sim_output/zdc_pi0/{DETECTOR_CONFIG}_rec_zdc_pi0_{P}GeV_{INDEX}.edm4eic.root",
P=[60, 80, 100, 130, 160],
- DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
+ INDEX=range(5),
+ ),
script="benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py",
output:
results_dir=directory("results/{DETECTOR_CONFIG}/zdc_pi0"),
diff --git a/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py b/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py
index 92adac55..1b54430e 100644
--- a/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py
+++ b/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py
@@ -21,13 +21,10 @@ import uproot as ur
arrays_sim={}
momenta=60, 80, 100, 130, 160,
for p in momenta:
- filename=f'sim_output/zdc_pi0/{config}_rec_zdc_pi0_{p}GeV.edm4hep.root'
- print("opening file", filename)
- events = ur.open(filename+':events')
- arrays_sim[p] = events.arrays()
- import gc
- gc.collect()
- print("read", filename)
+ arrays_sim[p] = ur.concatenate({
+ f'sim_output/zdc_pi0/{config}_rec_zdc_pi0_{p}GeV_{index}.edm4eic.root': 'events'
+ for index in range(5)
+ })
#energy res plot
fig,axs=plt.subplots(1,3, figsize=(24, 8))
diff --git a/benchmarks/zdc_pi0/config.yml b/benchmarks/zdc_pi0/config.yml
index a2fcb5de..ae2c1fd9 100644
--- a/benchmarks/zdc_pi0/config.yml
+++ b/benchmarks/zdc_pi0/config.yml
@@ -8,9 +8,8 @@ sim:zdc_pi0:
- P: 100
- P: 130
- P: 160
- timeout: 1 hours
script:
- - snakemake --cores 1 sim_output/zdc_pi0/epic_zdc_sipm_on_tile_only_rec_zdc_pi0_${P}GeV.edm4hep.root
+ - snakemake --cores 5 sim_output/zdc_pi0/epic_zdc_sipm_on_tile_only_rec_zdc_pi0_${P}GeV_{0,1,2,3,4}.edm4eic.root
retry:
max: 2
when:
diff --git a/benchmarks/zdc_sigma/Snakefile b/benchmarks/zdc_sigma/Snakefile
index a92a3781..80e8b25d 100644
--- a/benchmarks/zdc_sigma/Snakefile
+++ b/benchmarks/zdc_sigma/Snakefile
@@ -12,21 +12,19 @@ root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildca
rule zdc_sigma_simulate:
input:
- GEN_FILE="sim_output/zdc_sigma/sigma_decay_{P}GeV.hepmc"
+ GEN_FILE="sim_output/zdc_sigma/sigma_decay_{P}GeV.hepmc",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
params:
PHYSICS_LIST="FTFP_BERT"
output:
- SIM_FILE="sim_output/zdc_sigma/{DETECTOR_CONFIG}_sim_sigma_dec_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/zdc_sigma/{DETECTOR_CONFIG}_sim_sigma_dec_{P}GeV_{INDEX}.edm4hep.root"
shell:
"""
-if [[ {wildcards.P} -gt 220 ]]; then
- NEVENTS_SIM=500
-else
- NEVENTS_SIM=1000
-fi
+NEVENTS_SIM=200
# Running simulation
npsim \
--compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \
--numberOfEvents $NEVENTS_SIM \
--physicsList {params.PHYSICS_LIST} \
--inputFiles {input.GEN_FILE} \
@@ -35,24 +33,22 @@ npsim \
rule zdc_sigma_recon:
input:
- SIM_FILE="sim_output/zdc_sigma/{DETECTOR_CONFIG}_sim_sigma_dec_{P}GeV.edm4hep.root"
+ SIM_FILE="sim_output/zdc_sigma/{DETECTOR_CONFIG}_sim_sigma_dec_{P}GeV_{INDEX}.edm4hep.root"
output:
- REC_FILE="sim_output/zdc_sigma/{DETECTOR_CONFIG}_rec_sigma_dec_{P}GeV.edm4eic.root"
+ REC_FILE="sim_output/zdc_sigma/{DETECTOR_CONFIG}_rec_sigma_dec_{P}GeV_{INDEX}.edm4eic.root"
shell:
"""
-if [[ {wildcards.P} -gt 220 ]]; then
- NEVENTS_REC=500
-else
- NEVENTS_REC=1000
-fi
+NEVENTS_REC=200
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",
+ expand("sim_output/zdc_sigma/{DETECTOR_CONFIG}_rec_sigma_dec_{P}GeV_{INDEX}.edm4eic.root",
P=[100, 125, 150,175, 200, 225, 250, 275],
- DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
+ INDEX=range(5),
+ ),
script="benchmarks/zdc_sigma/analysis/sigma_plots.py",
output:
results_dir=directory("results/{DETECTOR_CONFIG}/zdc_sigma"),
diff --git a/benchmarks/zdc_sigma/analysis/gen_sigma_decay.cxx b/benchmarks/zdc_sigma/analysis/gen_sigma_decay.cxx
index 4da92801..4e808434 100644
--- a/benchmarks/zdc_sigma/analysis/gen_sigma_decay.cxx
+++ b/benchmarks/zdc_sigma/analysis/gen_sigma_decay.cxx
@@ -36,14 +36,12 @@ void gen_sigma_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "
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
@@ -76,10 +74,11 @@ void gen_sigma_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "
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++) {
+ int accepted_events = 0;
+ while (accepted_events < n_events) {
//Set the event number
- evt.set_event_number(events_parsed);
+ evt.set_event_number(accepted_events);
// FourVector(px,py,pz,e,pdgid,status)
// type 4 is beam
@@ -281,7 +280,7 @@ void gen_sigma_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "
evt.add_vertex(v_pi0_decay);
- if (events_parsed == 0) {
+ if (accepted_events == 0) {
std::cout << "First event: " << std::endl;
Print::listing(evt);
}
@@ -294,12 +293,12 @@ void gen_sigma_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "
hepmc_output.write_event(evt);
accepted_events ++;
}
- if (events_parsed % 1000 == 0) {
- std::cout << "Event: " << events_parsed << " ("<<accepted_events<<" accepted)"<< std::endl;
+ if (accepted_events % 1000 == 0) {
+ std::cout << "Event: " << accepted_events << std::endl;
}
evt.clear();
}
hepmc_output.close();
- std::cout << "Events generated: " << events_parsed << " ("<<accepted_events<<" accepted)"<< std::endl;
+ std::cout << "Events parsed and written: " << accepted_events << std::endl;
}
diff --git a/benchmarks/zdc_sigma/analysis/sigma_plots.py b/benchmarks/zdc_sigma/analysis/sigma_plots.py
index 6488b128..fc422ed6 100644
--- a/benchmarks/zdc_sigma/analysis/sigma_plots.py
+++ b/benchmarks/zdc_sigma/analysis/sigma_plots.py
@@ -19,13 +19,10 @@ 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)
+ arrays_sim[p] = ur.concatenate({
+ f'sim_output/zdc_sigma/{config}_rec_sigma_dec_{p}GeV_{index}.edm4eic.root': 'events'
+ for index in range(5)
+ })
def gauss(x, A,mu, sigma):
return A * np.exp(-(x-mu)**2/(2*sigma**2))
diff --git a/benchmarks/zdc_sigma/config.yml b/benchmarks/zdc_sigma/config.yml
index fc81428a..5a490510 100644
--- a/benchmarks/zdc_sigma/config.yml
+++ b/benchmarks/zdc_sigma/config.yml
@@ -11,9 +11,8 @@ sim:zdc_sigma:
- 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
+ - snakemake --cores 5 sim_output/zdc_sigma/epic_zdc_sipm_on_tile_only_rec_sigma_dec_${P}GeV_{0,1,2,3,4}.edm4eic.root
retry:
max: 2
when:
--
GitLab