Merge branch 'master' into pr/lfhcal_benchmarks

.gitlab-ci.yml

@@ -144,7 +144,9 @@ include:
   - local: 'benchmarks/zdc_sigma/config.yml'
   - local: 'benchmarks/zdc_lambda/config.yml'
   - local: 'benchmarks/insert_neutron/config.yml'
-
+  - local: 'benchmarks/femc_electron/config.yml'
+  - local: 'benchmarks/femc_photon/config.yml'
+  - local: 'benchmarks/femc_pi0/config.yml'
 deploy_results:
   allow_failure: true
   stage: deploy
@@ -167,6 +169,9 @@ deploy_results:
     - "collect_results:insert_muon"
     - "collect_results:zdc_photon"
     - "collect_results:zdc_pi0"
+    - "collect_results:femc_electron"
+    - "collect_results:femc_photon"
+    - "collect_results:femc_pi0"
   script:
     - echo "deploy results!"
     - find results -print | sort | tee summary.txt

Snakefile

@@ -15,6 +15,9 @@ include: "benchmarks/zdc_photon/Snakefile"
 include: "benchmarks/zdc_pi0/Snakefile"
 include: "benchmarks/zdc_sigma/Snakefile"
 include: "benchmarks/insert_neutron/Snakefile"
+include: "benchmarks/femc_electron/Snakefile"
+include: "benchmarks/femc_photon/Snakefile"
+include: "benchmarks/femc_pi0/Snakefile"
 
 use_s3 = config["remote_provider"].lower() == "s3"
 use_xrootd = config["remote_provider"].lower() == "xrootd"

benchmarks/backgrounds/config.yml

@@ -5,7 +5,7 @@ sim:backgrounds:
     - mkdir -p $LOCAL_DATA_PATH/input
     - ln -s $LOCAL_DATA_PATH/input input
     - |
-      snakemake -cache --cores 2 \
+      snakemake --cache --cores 2 \
         sim_output/$DETECTOR_CONFIG/backgrounds/EPIC/EVGEN/BACKGROUNDS/BEAMGAS/electron/GETaLM1.0.0-1.0/10GeV/GETaLM1.0.0-1.0_ElectronBeamGas_10GeV_foam_emin10keV_run001.edm4hep.root \
         sim_output/$DETECTOR_CONFIG/backgrounds/EPIC/EVGEN/DIS/NC/10x100/minQ2=1/pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_1.edm4hep.root \
         sim_output/$DETECTOR_CONFIG/backgrounds/EPIC/EVGEN/BACKGROUNDS/BEAMGAS/proton/pythia8.306-1.0/100GeV/pythia8.306-1.0_ProtonBeamGas_100GeV_run001.edm4hep.root

benchmarks/femc_electron/Snakefile

@@ -0,0 +1,58 @@
+rule femc_electron_generate:
+	input:
+                script="benchmarks/femc_electron/analysis/gen_particles.cxx",
+	params:
+		NEVENTS_GEN=1000,
+		th_max=28,
+		th_min=2.0
+	output:
+		GEN_FILE="sim_output/femc_electron/e-_{P}GeV.hepmc"
+	shell:
+		"""
+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 femc_electron_simulate:
+	input:
+		GEN_FILE="sim_output/femc_electron/e-_{P}GeV.hepmc"
+	params:
+		PHYSICS_LIST="FTFP_BERT"
+	output:
+		SIM_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_sim_e-_{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_electron_recon:
+        input:
+                SIM_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_sim_e-_{P}GeV.edm4hep.root"
+        output:
+                REC_FILE="sim_output/femc_electron/{DETECTOR_CONFIG}_rec_e-_{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,EcalEndcapPInsertRecHits,EcalEndcapPClusters  -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule femc_electron_analysis:
+	input:
+                expand("sim_output/femc_electron/{DETECTOR_CONFIG}_rec_e-_{P}GeV.edm4eic.root",
+		    P=[10, 20, 30, 40, 50, 60, 70, 80],
+		    DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+                script="benchmarks/femc_electron/analysis/femc_electron_plots.py",
+	output:
+		results_dir=directory("results/{DETECTOR_CONFIG}/femc_electron"),
+	shell:
+		"""
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""

benchmarks/femc_electron/analysis/femc_electron_plots.py

@@ -0,0 +1,222 @@
+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_electron/{config}_rec_e-_{p}GeV.edm4eic.root:events').arrays() for p in (10, 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")
+
+
+#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)):
+        if len(E[evt])==0:
+            continue
+        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(7, 85, 100)
+plt.plot(xx, 100*fnc(xx,*coeff), label=f'fit:{100*coeff[0]:.1f}%$\\oplus\\frac{{{100*coeff[1]:.0f}\\%}}{{\\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")
+
+#energy res in eta bins
+fig, axs=plt.subplots(1,2, figsize=(16,8))
+
+partitions=[1.5, 2.0, 2.5, 3.0]
+for eta_min, eta_max in zip(partitions[:-1], partitions[1:]):
+    pvals=[]
+    res=[]
+    dres=[]
+    scale=[]
+    dscale=[]
+    for p in arrays_sim:
+        bins=np.linspace(15*p/20,22*p/20, 50)
+        eta_truth=arrays_sim[p]['eta_truth']
+        y,x=np.histogram(ak.flatten(arrays_sim[p][(eta_truth<eta_max)&(eta_truth>eta_min)]['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)
+
+        try:
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+            if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
+                continue
+            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)
+        except:
+            pass
+    plt.sca(axs[0])
+    plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o', label=f'${eta_min:.1f}<\\eta<{eta_max:.1f}$')
+
+    plt.sca(axs[1])
+    plt.errorbar(pvals, 100*np.array(scale), 100*np.array(dscale), ls='', marker='o', label=f'${eta_min:.1f}<\\eta<{eta_max:.1f}$')
+
+plt.sca(axs[0])
+plt.legend()
+plt.ylim(0)
+plt.ylabel("E resolution [%]")
+plt.xlabel("E truth [GeV]")
+
+plt.sca(axs[1])
+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_eta_partitions.pdf")