From 97ec00d29168962d2868f1d70e00cfa01a28e1a4 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 30 Sep 2024 06:26:24 -0400 Subject: [PATCH] parallelize zdc and insert simulations --- benchmarks/insert_muon/Snakefile | 17 +- benchmarks/insert_muon/analysis/muon_plots.py | 169 +++++++++--------- benchmarks/insert_muon/config.yml | 3 +- benchmarks/insert_neutron/Snakefile | 17 +- .../insert_neutron/analysis/neutron_plots.py | 6 +- benchmarks/insert_neutron/config.yml | 3 +- benchmarks/zdc_lambda/Snakefile | 27 ++- .../zdc_lambda/analysis/lambda_plots.py | 11 +- benchmarks/zdc_lambda/config.yml | 3 +- benchmarks/zdc_photon/Snakefile | 31 ++-- .../zdc_photon/analysis/zdc_photon_plots.py | 13 +- benchmarks/zdc_photon/config.yml | 3 +- benchmarks/zdc_pi0/Snakefile | 17 +- benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py | 11 +- benchmarks/zdc_pi0/config.yml | 3 +- benchmarks/zdc_sigma/Snakefile | 25 ++- benchmarks/zdc_sigma/analysis/sigma_plots.py | 11 +- benchmarks/zdc_sigma/config.yml | 3 +- 18 files changed, 172 insertions(+), 201 deletions(-) diff --git a/benchmarks/insert_muon/Snakefile b/benchmarks/insert_muon/Snakefile index b36026b..622733c 100644 --- a/benchmarks/insert_muon/Snakefile +++ b/benchmarks/insert_muon/Snakefile @@ -18,13 +18,14 @@ rule insert_muon_simulate: 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 +34,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 4c45aac..052091e 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}_rec_mu-_{p}GeV_{index}.edm4eic.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 6f7681f..a2a6ff8 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_{1,2,3,4,5}.edm4hep.root retry: max: 2 when: diff --git a/benchmarks/insert_neutron/Snakefile b/benchmarks/insert_neutron/Snakefile index 55e3ff9..d14a7be 100644 --- a/benchmarks/insert_neutron/Snakefile +++ b/benchmarks/insert_neutron/Snakefile @@ -19,13 +19,14 @@ rule insert_neutron_simulate: 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 +35,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 9314a97..b9d2002 100644 --- a/benchmarks/insert_neutron/analysis/neutron_plots.py +++ b/benchmarks/insert_neutron/analysis/neutron_plots.py @@ -18,8 +18,10 @@ #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 aac0a73..ec75faa 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_{1,2,3,4,5}.edm4eic.root retry: max: 2 when: diff --git a/benchmarks/zdc_lambda/Snakefile b/benchmarks/zdc_lambda/Snakefile index 2156333..f7f934d 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: @@ -16,17 +16,14 @@ rule zdc_lambda_simulate: 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 +32,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/lambda_plots.py b/benchmarks/zdc_lambda/analysis/lambda_plots.py index 88aff7f..bf855d9 100644 --- a/benchmarks/zdc_lambda/analysis/lambda_plots.py +++ b/benchmarks/zdc_lambda/analysis/lambda_plots.py @@ -19,13 +19,10 @@ 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 396b192..be6c0bf 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_{1,2,3,4,5}.edm4eic.root retry: max: 2 when: diff --git a/benchmarks/zdc_photon/Snakefile b/benchmarks/zdc_photon/Snakefile index cfb5124..0c551ab 100644 --- a/benchmarks/zdc_photon/Snakefile +++ b/benchmarks/zdc_photon/Snakefile @@ -8,11 +8,7 @@ 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})' """ @@ -23,18 +19,15 @@ rule zdc_photon_simulate: 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 +35,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}.edm4hep.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}.edm4hep.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 b06cf43..94a9310 100644 --- a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py +++ b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py @@ -23,14 +23,11 @@ def gauss(x, A,mu, sigma): 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 946650a..094cdc0 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_{1,2,3,4,5}.edm4hep.root retry: max: 2 when: diff --git a/benchmarks/zdc_pi0/Snakefile b/benchmarks/zdc_pi0/Snakefile index 57fd2c9..8bc5c40 100644 --- a/benchmarks/zdc_pi0/Snakefile +++ b/benchmarks/zdc_pi0/Snakefile @@ -17,13 +17,14 @@ rule zdc_pi0_simulate: 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 +33,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}.edm4hep.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}.edm4hep.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 92adac5..1b54430 100644 --- a/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py +++ b/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py @@ -21,13 +21,10 @@ def gauss(x, A,mu, sigma): 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 a2fcb5d..d645f91 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_{1,2,3,4,5}.edm4hep.root retry: max: 2 when: diff --git a/benchmarks/zdc_sigma/Snakefile b/benchmarks/zdc_sigma/Snakefile index a92a378..6f6e4ab 100644 --- a/benchmarks/zdc_sigma/Snakefile +++ b/benchmarks/zdc_sigma/Snakefile @@ -16,17 +16,14 @@ rule zdc_sigma_simulate: 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 +32,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/sigma_plots.py b/benchmarks/zdc_sigma/analysis/sigma_plots.py index 6488b12..d1ce8b1 100644 --- a/benchmarks/zdc_sigma/analysis/sigma_plots.py +++ b/benchmarks/zdc_sigma/analysis/sigma_plots.py @@ -19,13 +19,10 @@ 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_zdc_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 fc81428..922804c 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_{1,2,3,4,5}.edm4eic.root retry: max: 2 when: