Skip to content

Commit

Permalink
added benchmarks for pi0 and photon as well
Browse files Browse the repository at this point in the history
  • Loading branch information
sebouh137 committed Sep 26, 2024
1 parent a138a3b commit 6926300
Show file tree
Hide file tree
Showing 11 changed files with 985 additions and 12 deletions.
10 changes: 5 additions & 5 deletions benchmarks/femc_electron/Snakefile
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
rule generate:
rule femc_electron_generate:
input:
script="src/gen_particles.cxx",
script="benchmarks/femc_electron/analysis/gen_particles.cxx",
params:
NEVENTS_GEN=1000,
th_max=28,
Expand All @@ -13,7 +13,7 @@ mkdir -p sim_output/femc_electron
root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "e-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
"""

rule simulate:
rule femc_electron_simulate:
input:
GEN_FILE="sim_output/femc_electron/e-_{P}GeV.hepmc"
params:
Expand All @@ -32,7 +32,7 @@ npsim \
--outputFile {output.SIM_FILE}
"""

rule recon:
rule femc_electron_recon:
input:
SIM_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_sim_e-_{P}GeV.edm4hep.root"
output:
Expand All @@ -43,7 +43,7 @@ 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,EcalEndcapPInsertRecHits,EcalEndcapPClusters -Pjana:nevents=$NEVENTS_REC
"""

rule analysis:
rule femc_electron_analysis:
input:
expand("sim_output/femc_electron/{DETECTOR_CONFIG}_rec_e-_{P}GeV.edm4hep.root",
P=[20, 30, 40, 50, 60, 70, 80],
Expand Down
32 changes: 25 additions & 7 deletions benchmarks/femc_electron/analysis/femc_electron_plots.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
pass

import uproot as ur
arrays_sim={p:ur.open(f'/Users/spaul/EIC/detector_benchmarks/sim_output/femc_electron/epic_craterlake_rec_e-_{p}GeV.edm4hep.root:events').arrays() for p in (20, 30, 40, 50, 60,70,80)}
arrays_sim={p:ur.open(f'/sim_output/femc_electron/{config}_rec_e-_{p}GeV.edm4hep.root:events').arrays() for p in (20, 30, 40, 50, 60,70,80)}

for p in arrays_sim:
array=arrays_sim[p]
Expand Down Expand Up @@ -66,24 +66,42 @@
pvals=[]

#number of hits per cluster
fig, axs=plt.subplots(1,2, figsize=(16,8))
avgs=[]
stds=[]
pvals=[]

for p in arrays_sim:

a=arrays_sim[p]
n=[]
nn=-a['EcalEndcapPClusters.hits_begin']+a['EcalEndcapPClusters.hits_end']
E=a['EcalEndcapPClusters.energy']
for evt in range(len(array)):
maxE=np.max(E[evt])
found=False
for i in range(len(E[evt])):
if E[evt][i]==maxE:
n.append(nn[evt][i])
found=True
break
#if not found:
# n.append(0)

if p ==50:
plt.sca(axs[0])
y,x,_=plt.hist(ak.flatten(-a['EcalEndcapPClusters.hits_begin']+a['EcalEndcapPClusters.hits_end']),
range=(0,100), bins=21, histtype='step', label=f"E={p} GeV")
y,x,_=plt.hist(n, range=(0,100), bins=100, histtype='step', label=f"E={p} GeV")
plt.ylabel("events")
plt.xlabel("# hits per Ecal cluster")
plt.xlabel("# hits in cluster")
plt.title(f"e-, E={p} GeV")

n=ak.flatten(-a['EcalEndcapPClusters.hits_begin']+a['EcalEndcapPClusters.hits_end'])
pvals.append(p)
avgs.append(np.mean(n))
stds.append(np.std(n))

plt.sca(axs[1])
plt.errorbar(pvals, avgs, stds, marker='o',ls='')
plt.xlabel("E [GeV]")
plt.ylabel("# hits per Ecal cluster [mean$\\pm$std]")
plt.ylabel("# hits in cluster [mean$\\pm$std]")
plt.ylim(0)
plt.savefig(outdir+"/nhits_per_cluster.pdf")

Expand Down
59 changes: 59 additions & 0 deletions benchmarks/femc_photon/Snakefile
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
rule femc_photon_generate:
input:
script="benchmarks/femc_photon/analysis/gen_particles.cxx",
params:
NEVENTS_GEN=1000,
th_max=28,
th_min=2.0
output:
GEN_FILE="sim_output/femc_photon/photon_{P}GeV.hepmc"
shell:
"""
mkdir -p sim_output/femc_photon
root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "gamma", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
"""

rule femc_photon_simulate:
input:
GEN_FILE="sim_output/femc_photon/photon_{P}GeV.hepmc"
params:
PHYSICS_LIST="FTFP_BERT"
output:
SIM_FILE="sim_output/femc_photon/{DETECTOR_CONFIG}_sim_photon_{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 femc_photon_recon:
input:
SIM_FILE="sim_output/femc_photon/{DETECTOR_CONFIG}_sim_photon_{P}GeV.edm4hep.root"
output:
REC_FILE="sim_output/femc_photon/{DETECTOR_CONFIG}_rec_photon_{P}GeV.edm4hep.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,EcalEndcapPInsertRecHits,EcalEndcapPClusters -Pjana:nevents=$NEVENTS_REC
"""

rule femc_photon_analysis:
input:
expand("sim_output/femc_photon/{DETECTOR_CONFIG}_rec_photon_{P}GeV.edm4hep.root",
P=[20, 30, 40, 50, 60, 70, 80],
DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
particle=["{particle}"]),
script="benchmarks/femc_photon/analysis/femc_photon_plots.py",
output:
results_dir=directory("results/{DETECTOR_CONFIG}/femc_photon"),
shell:
"""
mkdir -p {output.results_dir}
python {input.script} {output.results_dir}
"""
171 changes: 171 additions & 0 deletions benchmarks/femc_photon/analysis/femc_photon_plots.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,171 @@
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}/{benchmark_name}
outdir=sys.argv[1]+"/"
try:
import os
os.mkdir(outdir[:-1])
except:
pass

import uproot as ur
arrays_sim={p:ur.open(f'sim_output/femc_photon/{config}_rec_e-_{p}GeV.edm4hep.root:events').arrays() for p in (20, 30, 40, 50, 60,70,80)}

for p in arrays_sim:
array=arrays_sim[p]
tilt=-.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['nclust_endcap']=[len(array['EcalEndcapPClusters.energy'][i]) for i in range(len(array))]

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)

#number of clusters
plt.figure()
for eta_min, eta_max, field in (1.5, 2.8, 'nclust_endcap'),:
for p in arrays_sim:
array=arrays_sim[p]
plt.hist(array[field][(array['eta_truth']>eta_min)&(array['eta_truth']<eta_max)],
bins=np.linspace(0,10,11), histtype='step', label=f'{p} GeV', density=True)
plt.ylabel("events")
plt.xlabel("# of Ecal clusters")
plt.legend()
plt.savefig(outdir+f"/{field}.pdf")

fig, axs=plt.subplots(1,2, figsize=(16,8))
avgs=[]
stds=[]
pvals=[]

#number of hits per cluster
fig, axs=plt.subplots(1,2, figsize=(16,8))
avgs=[]
stds=[]
pvals=[]

for p in arrays_sim:

a=arrays_sim[p]
n=[]
nn=-a['EcalEndcapPClusters.hits_begin']+a['EcalEndcapPClusters.hits_end']
E=a['EcalEndcapPClusters.energy']
for evt in range(len(array)):
maxE=np.max(E[evt])
found=False
for i in range(len(E[evt])):
if E[evt][i]==maxE:
n.append(nn[evt][i])
found=True
break
#if not found:
# n.append(0)

if p ==50:
plt.sca(axs[0])
y,x,_=plt.hist(n, range=(0,100), bins=100, histtype='step', label=f"E={p} GeV")
plt.ylabel("events")
plt.xlabel("# hits in cluster")
plt.title(f"e-, E={p} GeV")
pvals.append(p)
avgs.append(np.mean(n))
stds.append(np.std(n))

plt.sca(axs[1])
plt.errorbar(pvals, avgs, stds, marker='o',ls='')
plt.xlabel("E [GeV]")
plt.ylabel("# hits in cluster [mean$\\pm$std]")
plt.ylim(0)
plt.savefig(outdir+"/nhits_per_cluster.pdf")


#energy resolution
def gauss(x, A,mu, sigma):
return A * np.exp(-(x-mu)**2/(2*sigma**2))
from scipy.optimize import curve_fit

fig, axs=plt.subplots(1,3, figsize=(24,8))
pvals=[]
res=[]
dres=[]
scale=[]
dscale=[]
for p in arrays_sim:
bins=np.linspace(15*p/20,22*p/20, 50)
if p==50:
plt.sca(axs[0])
plt.title(f"E={p} GeV")
y,x,_=plt.hist(ak.flatten(arrays_sim[p]['EcalEndcapPClusters.energy']), bins=bins, histtype='step')
plt.ylabel("events")
plt.xlabel("$E^{rec}_e$ [GeV]")
else:
y,x=np.histogram(ak.flatten(arrays_sim[p]['EcalEndcapPClusters.energy']), bins=bins)
bcs=(x[1:]+x[:-1])/2

fnc=gauss
slc=abs(bcs-p)<3
sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
p0=(100, p, 3)

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 p==50:
xx=np.linspace(15*p/20,22*p/20, 100)

plt.plot(xx, fnc(xx,*coeff), label=f"$\\sigma_E/E={abs(coeff[2])/coeff[1]*100:.1f}\%$")
plt.axvline(p, color='g', ls='--', alpha=0.7)
plt.legend()
#plt.xlim(0,60)
#plt.show()
pvals.append(p)
res.append(abs(coeff[2])/coeff[1])
dres.append(np.sqrt(var_matrix[2][2])/coeff[1])
scale.append(abs(coeff[1])/p)
dscale.append(np.sqrt(var_matrix[1][1])/p)
plt.sca(axs[1])
plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o')
fnc = lambda E, a, b: np.hypot(a,b/np.sqrt(E))
p0=(.05, .12)
coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres)
xx=np.linspace(15, 85, 100)
plt.plot(xx, 100*fnc(xx,*coeff), label=f'fit:{100*coeff[0]:.0f}%$\\oplus\\frac{{{100*coeff[1]:.1f}\\%}}{{\\sqrt{{E}}}}$')
plt.legend()
plt.ylim(0)
plt.ylabel("E resolution [%]")
plt.xlabel("E truth [GeV]")
plt.sca(axs[2])

plt.errorbar(pvals, 100*np.array(scale), 100*np.array(dscale), ls='', marker='o')
plt.ylabel("energy scale [%]")
plt.xlabel("E truth [GeV]")
plt.axhline(100, color='0.5', alpha=0.5, ls='--')
plt.ylim(0, 110)
plt.tight_layout()
plt.savefig(outdir+"/energy_res.pdf")
Loading

0 comments on commit 6926300

Please sign in to comment.