diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index a1dc5d8c3af3903a3b73bd1a7b487ece7ada20ee..a56c79d26301ea6fc1afd87145333ffba52a8b8b 100644
--- a/.gitlab-ci.yml
+++ b/.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
diff --git a/Snakefile b/Snakefile
index bce37f8dcab0328b0bfa5aeb6a82a9b8afa19a20..5432328496cf9d827cdcd2f46ef6a1bb18065f6a 100644
--- a/Snakefile
+++ b/Snakefile
@@ -7,7 +7,6 @@ include: "benchmarks/ecal_gaps/Snakefile"
include: "benchmarks/material_scan/Snakefile"
include: "benchmarks/tracking_performances/Snakefile"
include: "benchmarks/tracking_performances_dis/Snakefile"
-include: "benchmarks/zdc_lyso/Snakefile"
include: "benchmarks/insert_muon/Snakefile"
include: "benchmarks/zdc_lambda/Snakefile"
include: "benchmarks/zdc_lyso/Snakefile"
@@ -15,7 +14,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"
diff --git a/benchmarks/femc_electron/Snakefile b/benchmarks/femc_electron/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..0e8dd0842b2a79fd67f59ccbae14bea41cd0ba22
--- /dev/null
+++ b/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}
+"""
diff --git a/benchmarks/femc_electron/analysis/femc_electron_plots.py b/benchmarks/femc_electron/analysis/femc_electron_plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..cc34633dc30887057bda2e8331f1ed509c6f33b5
--- /dev/null
+++ b/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")
diff --git a/benchmarks/femc_electron/analysis/gen_particles.cxx b/benchmarks/femc_electron/analysis/gen_particles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0877a85223bbc1af27af91e7dedfb212de7cb20b
--- /dev/null
+++ b/benchmarks/femc_electron/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/femc_electron/config.yml b/benchmarks/femc_electron/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..a600454d1f2a7b20b2450d22f4ed144717d2d901
--- /dev/null
+++ b/benchmarks/femc_electron/config.yml
@@ -0,0 +1,36 @@
+sim:femc_electron:
+ stage: simulate
+ extends: .det_benchmark
+ parallel:
+ matrix:
+ - P: 10
+ - P: 20
+ - P: 30
+ - P: 40
+ - P: 50
+ - P: 60
+ - P: 70
+ - P: 80
+ timeout: 1 hours
+ script:
+ - snakemake --cores 1 sim_output/femc_electron/epic_craterlake_rec_e-_${P}GeV.edm4eic.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:femc_electron:
+ stage: benchmarks
+ extends: .det_benchmark
+ needs: ["sim:femc_electron"]
+ script:
+ - snakemake --cores 1 results/epic_craterlake/femc_electron
+
+collect_results:femc_electron:
+ stage: collect
+ needs: ["bench:femc_electron"]
+ 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/femc_electron
+ - mv results{_save,}/
diff --git a/benchmarks/femc_photon/Snakefile b/benchmarks/femc_photon/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..40524365f7674bc92cb04c4d45160975644a858c
--- /dev/null
+++ b/benchmarks/femc_photon/Snakefile
@@ -0,0 +1,58 @@
+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.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_photon_analysis:
+ input:
+ expand("sim_output/femc_photon/{DETECTOR_CONFIG}_rec_photon_{P}GeV.edm4eic.root",
+ P=[10, 20, 30, 40, 50, 60, 70, 80],
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ 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}
+"""
diff --git a/benchmarks/femc_photon/analysis/femc_photon_plots.py b/benchmarks/femc_photon/analysis/femc_photon_plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..8170f4a1d0ecd006a5f511bab246cdd0124f9cdd
--- /dev/null
+++ b/benchmarks/femc_photon/analysis/femc_photon_plots.py
@@ -0,0 +1,220 @@
+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_photon_{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"photon, 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}_\\gamma$ [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")
diff --git a/benchmarks/femc_photon/analysis/femc_photon_plots.py~ b/benchmarks/femc_photon/analysis/femc_photon_plots.py~
new file mode 100644
index 0000000000000000000000000000000000000000..4227b4e3d910ec8dd65f16bf527661f273190fee
--- /dev/null
+++ b/benchmarks/femc_photon/analysis/femc_photon_plots.py~
@@ -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/epic_craterlake_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")
diff --git a/benchmarks/femc_photon/analysis/gen_particles.cxx b/benchmarks/femc_photon/analysis/gen_particles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0877a85223bbc1af27af91e7dedfb212de7cb20b
--- /dev/null
+++ b/benchmarks/femc_photon/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/femc_photon/config.yml b/benchmarks/femc_photon/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..a367e5093dccb1dc5b9642f2fede452a5c46f85a
--- /dev/null
+++ b/benchmarks/femc_photon/config.yml
@@ -0,0 +1,36 @@
+sim:femc_photon:
+ stage: simulate
+ extends: .det_benchmark
+ parallel:
+ matrix:
+ - P: 10
+ - P: 20
+ - P: 30
+ - P: 40
+ - P: 50
+ - P: 60
+ - P: 70
+ - P: 80
+ timeout: 1 hours
+ script:
+ - snakemake --cores 1 sim_output/femc_photon/epic_craterlake_rec_photon_${P}GeV.edm4eic.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:femc_photon:
+ stage: benchmarks
+ extends: .det_benchmark
+ needs: ["sim:femc_photon"]
+ script:
+ - snakemake --cores 1 results/epic_craterlake/femc_photon
+
+collect_results:femc_photon:
+ stage: collect
+ needs: ["bench:femc_photon"]
+ 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/femc_photon
+ - mv results{_save,}/
diff --git a/benchmarks/femc_pi0/Snakefile b/benchmarks/femc_pi0/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..4936040cb832a72cfef4a43fbb4dfc02296695b2
--- /dev/null
+++ b/benchmarks/femc_pi0/Snakefile
@@ -0,0 +1,58 @@
+rule femc_pi0_generate:
+ input:
+ script="benchmarks/femc_pi0/analysis/gen_particles.cxx",
+ params:
+ NEVENTS_GEN=1000,
+ th_max=28,
+ th_min=2.0
+ output:
+ GEN_FILE="sim_output/femc_pi0/pi0_{P}GeV.hepmc"
+ shell:
+ """
+mkdir -p sim_output/femc_pi0
+root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "pi0", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+"""
+
+rule femc_pi0_simulate:
+ input:
+ GEN_FILE="sim_output/femc_pi0/pi0_{P}GeV.hepmc"
+ params:
+ PHYSICS_LIST="FTFP_BERT"
+ output:
+ SIM_FILE="sim_output/femc_pi0/{DETECTOR_CONFIG}_sim_pi0_{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_pi0_recon:
+ input:
+ SIM_FILE="sim_output/femc_pi0/{DETECTOR_CONFIG}_sim_pi0_{P}GeV.edm4hep.root"
+ output:
+ REC_FILE="sim_output/femc_pi0/{DETECTOR_CONFIG}_rec_pi0_{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_pi0_analysis:
+ input:
+ expand("sim_output/femc_pi0/{DETECTOR_CONFIG}_rec_pi0_{P}GeV.edm4eic.root",
+ P=[10, 20, 30, 40, 50, 60, 70, 80],
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ script="benchmarks/femc_pi0/analysis/femc_pi0_plots.py",
+ output:
+ results_dir=directory("results/{DETECTOR_CONFIG}/femc_pi0"),
+ shell:
+ """
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""
diff --git a/benchmarks/femc_pi0/analysis/femc_pi0_plots.py b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..cdf9d2319c943721bb8fdf281dd34d7bda7fd069
--- /dev/null
+++ b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py
@@ -0,0 +1,257 @@
+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_pi0/{config}_rec_pi0_{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 clusters
+plt.figure()
+for eta_min, eta_max, field in (1.5, 2.8, 'nclust_endcap'),:
+ pvals=[]
+ f0=[]
+ f1=[]
+ f2=[]
+ f3=[]
+ f4=[]
+ f5p=[]
+ for p in arrays_sim:
+ array=arrays_sim[p]
+ n=array[field][(array['eta_truth']>eta_min)&(array['eta_truth']<eta_max)]
+ pvals.append(p)
+ tot=len(n)
+ f0.append(np.sum(1*(n==0))/tot)
+ f1.append(np.sum(1*(n==1))/tot)
+ f2.append(np.sum(1*(n==2))/tot)
+ f3.append(np.sum(1*(n==3))/tot)
+ f4.append(np.sum(1*(n==4))/tot)
+ f5p.append(np.sum(1*(n>=5))/tot)
+ plt.errorbar(pvals, f0, label="$n_{cl}$=0")
+ plt.errorbar(pvals, f1, label="$n_{cl}$=1")
+ plt.errorbar(pvals, f2, label="$n_{cl}$=2")
+ plt.errorbar(pvals, f3, label="$n_{cl}$=3")
+ plt.errorbar(pvals, f4, label="$n_{cl}$=4")
+ plt.errorbar(pvals, f5p, label="$n_{cl}\\geq$5")
+ plt.legend()
+ plt.ylabel("fraction of events")
+ plt.xlabel("$E_{\\pi^0}$ [GeV]")
+ plt.ylim(0, 1)
+ plt.legend(fontsize=20, ncol=2)
+ plt.savefig(outdir+f"/nclust.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"$\\pi^0$, 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(np.sum(arrays_sim[p]['EcalEndcapPClusters.energy'], axis=-1), bins=bins, histtype='step')
+ plt.ylabel("events")
+ plt.xlabel("$E^{rec}_{\\pi^0}$ [GeV]")
+ else:
+ y,x=np.histogram(np.sum(arrays_sim[p]['EcalEndcapPClusters.energy'], axis=-1), 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]:.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(np.sum(arrays_sim[p][(eta_truth<eta_max)&(eta_truth>eta_min)]['EcalEndcapPClusters.energy'], axis=-1), 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")
diff --git a/benchmarks/femc_pi0/analysis/gen_particles.cxx b/benchmarks/femc_pi0/analysis/gen_particles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0877a85223bbc1af27af91e7dedfb212de7cb20b
--- /dev/null
+++ b/benchmarks/femc_pi0/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/femc_pi0/config.yml b/benchmarks/femc_pi0/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..bcbd90afd07241217dc64b712fc209a3cf0bf84d
--- /dev/null
+++ b/benchmarks/femc_pi0/config.yml
@@ -0,0 +1,36 @@
+sim:femc_pi0:
+ stage: simulate
+ extends: .det_benchmark
+ parallel:
+ matrix:
+ - P: 10
+ - P: 20
+ - P: 30
+ - P: 40
+ - P: 50
+ - P: 60
+ - P: 70
+ - P: 80
+ timeout: 1 hours
+ script:
+ - snakemake --cores 1 sim_output/femc_pi0/epic_craterlake_rec_pi0_${P}GeV.edm4eic.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:femc_pi0:
+ stage: benchmarks
+ extends: .det_benchmark
+ needs: ["sim:femc_pi0"]
+ script:
+ - snakemake --cores 1 results/epic_craterlake/femc_pi0
+
+collect_results:femc_pi0:
+ stage: collect
+ needs: ["bench:femc_pi0"]
+ 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/femc_pi0
+ - mv results{_save,}/