Merge remote-tracking branch 'origin/master' into beamline_benchmark

.github/workflows/mirror.yaml

@@ -6,7 +6,7 @@ on:
   workflow_dispatch:
 
 concurrency:
-  group: mirror
+  group: mirror-${{ github.event_name }}
   cancel-in-progress: false
 
 jobs:

.gitlab-ci.yml

@@ -3,7 +3,6 @@ image: ${BENCHMARKS_REGISTRY}/${BENCHMARKS_CONTAINER}:${BENCHMARKS_TAG}
 variables:
   DETECTOR: epic
   DETECTOR_CONFIG: epic_craterlake
-  DETECTOR_REPOSITORYURL: 'https://github.com/eic/epic.git'
   GITHUB_SHA: ''
   GITHUB_REPOSITORY: ''
 
@@ -89,13 +88,6 @@ common:setup:
       echo "BENCHMARKS_CONTAINER: ${BENCHMARKS_CONTAINER}"
       echo "BENCHMARKS_REGISTRY: ${BENCHMARKS_REGISTRY}"
     - source setup/bin/env.sh && ./setup/bin/install_common.sh
-
-
-common:detector:
-  stage: initialize
-  needs: ["common:setup"]
-  script:
-    - source .local/bin/env.sh && build_detector.sh
     - mkdir_local_data_link sim_output
     - mkdir -p results
     - mkdir -p config
@@ -103,7 +95,7 @@ common:detector:
 
 get_data:
   stage: data_init
-  needs: ["common:detector"]
+  needs: ["common:setup"]
   script:
     - source .local/bin/env.sh
     - ln -s "${LOCAL_DATA_PATH}/sim_output" sim_output
@@ -112,10 +104,11 @@ get_data:
 
 .det_benchmark:
   needs:
-    - ["get_data","common:detector"]
+    - ["get_data","common:setup"]
   before_script:
     - mc config host add S3 https://eics3.sdcc.bnl.gov:9000 ${S3_ACCESS_KEY} ${S3_SECRET_KEY}
     - source .local/bin/env.sh
+    - source /opt/detector/epic-main/setup.sh
     - ls -lrtha 
     - ln -s "${LOCAL_DATA_PATH}/sim_output" sim_output
     - ln -s "${LOCAL_DATA_PATH}/datasets/data" data
@@ -136,23 +129,55 @@ include:
   - local: 'benchmarks/ecal_gaps/config.yml'
   - local: 'benchmarks/tracking_detectors/config.yml'
   - local: 'benchmarks/tracking_performances/config.yml'
+  - local: 'benchmarks/tracking_performances_dis/config.yml'
   - local: 'benchmarks/barrel_ecal/config.yml'
   - local: 'benchmarks/barrel_hcal/config.yml'
+  - local: 'benchmarks/lfhcal/config.yml'
   - local: 'benchmarks/zdc/config.yml'
   - local: 'benchmarks/zdc_lyso/config.yml'
-  - local: 'benchmarks/material_maps/config.yml'
+  - local: 'benchmarks/zdc_photon/config.yml'
+  - local: 'benchmarks/zdc_pi0/config.yml'
   - local: 'benchmarks/material_scan/config.yml'
   - local: 'benchmarks/pid/config.yml'
   - local: 'benchmarks/timing/config.yml'
   - local: 'benchmarks/b0_tracker/config.yml'
-
+  - local: 'benchmarks/insert_muon/config.yml'
+  - 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
   needs:
-    - ["collect_results:zdc","collect_results:barrel_ecal","collect_results:barrel_hcal","collect_results:material_scan"]
+    - "collect_results:backgrounds"
+    - "collect_results:backwards_ecal"
+    - "collect_results:barrel_ecal"
+    - "collect_results:barrel_hcal"
+    - "collect_results:ecal_gaps"
+    - "collect_results:material_scan"
+    - "collect_results:pid"
+    - "collect_results:tracking_performance"
+    - "collect_results:tracking_performance_campaigns"
+    - "collect_results:zdc_sigma"
+    - "collect_results:zdc_lambda"
+    - "collect_results:insert_neutron"
+    - "collect_results:tracking_performances_dis"
+    - "collect_results:zdc"
+    - "collect_results:zdc_lyso"
+    - "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
+    - xrdfs $XROOTD_RW_ENDPOINT mkdir $XROOTD_OUTPUT_PREFIX/pipeline-$CI_PIPELINE_ID
+    - xrdcp -r results $XROOTD_RW_ENDPOINT/$XROOTD_OUTPUT_PREFIX/pipeline-$CI_PIPELINE_ID
 
 benchmarks:detector:success:
   stage: status-report

Snakefile

@@ -6,8 +6,19 @@ include: "benchmarks/barrel_ecal/Snakefile"
 include: "benchmarks/ecal_gaps/Snakefile"
 include: "benchmarks/material_scan/Snakefile"
 include: "benchmarks/tracking_performances/Snakefile"
+include: "benchmarks/tracking_performances_dis/Snakefile"
+include: "benchmarks/lfhcal/Snakefile"
+include: "benchmarks/insert_muon/Snakefile"
+include: "benchmarks/zdc_lambda/Snakefile"
 include: "benchmarks/zdc_lyso/Snakefile"
 include: "benchmarks/beamline/Snakefile"
+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"
@@ -25,6 +36,9 @@ def get_remote_path(path):
 rule fetch_epic:
     output:
         filepath="EPIC/{PATH}"
+    params:
+        # wildcards are not included in hash for caching, we need to add them as params
+        PATH=lambda wildcards: wildcards.PATH
     cache: True
     retries: 3
     shell: """

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

+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

+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")

benchmarks/femc_electron/analysis/gen_particles.cxx

+#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;
+}

benchmarks/femc_electron/config.yml

+sim:femc_electron:
+  extends: .det_benchmark
+  stage: simulate
+  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:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs: ["sim:femc_electron"]
+  script:
+    - snakemake --cores 1 results/epic_craterlake/femc_electron
+
+collect_results:femc_electron:
+  extends: .det_benchmark
+  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,}/

benchmarks/femc_photon/Snakefile

+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}
+"""

benchmarks/femc_photon/analysis/femc_photon_plots.py

+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")

benchmarks/femc_photon/analysis/femc_photon_plots.py~

+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")

benchmarks/femc_photon/analysis/gen_particles.cxx

+#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;
+}

benchmarks/femc_photon/config.yml

+sim:femc_photon:
+  extends: .det_benchmark
+  stage: simulate
+  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:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs: ["sim:femc_photon"]
+  script:
+    - snakemake --cores 1 results/epic_craterlake/femc_photon
+
+collect_results:femc_photon:
+  extends: .det_benchmark
+  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,}/

benchmarks/femc_pi0/Snakefile

+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}
+"""

benchmarks/femc_pi0/analysis/femc_pi0_plots.py

+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")

benchmarks/femc_pi0/analysis/gen_particles.cxx

+#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;
+}

benchmarks/femc_pi0/config.yml

+sim:femc_pi0:
+  extends: .det_benchmark
+  stage: simulate
+  parallel:
+    matrix:
+      - P: 10
+      - P: 20
+      - P: 30
+      - P: 40
+      - P: 50
+      - P: 60
+      - P: 70
+      - P: 80
+  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:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs: ["sim:femc_pi0"]
+  script:
+    - snakemake --cores 1 results/epic_craterlake/femc_pi0
+
+collect_results:femc_pi0:
+  extends: .det_benchmark
+  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,}/

benchmarks/insert_muon/Snakefile

+rule insert_muon_generate:
+	input:
+                script="benchmarks/insert_muon/analysis/gen_particles.cxx",
+	params:
+		NEVENTS_GEN=5000,
+		th_max=7.0,
+		th_min=1.7
+	output:
+		GEN_FILE="sim_output/insert_muon/mu-_{P}GeV.hepmc"
+	shell:
+		"""
+root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "mu-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+"""
+
+rule insert_muon_simulate:
+	input:
+                GEN_FILE="sim_output/insert_muon/mu-_{P}GeV.hepmc",
+                warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+	params:
+		PHYSICS_LIST="FTFP_BERT"
+	output:
+		SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV_{INDEX}.edm4hep.root"
+	shell:
+		"""
+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} \
+   --outputFile {output.SIM_FILE}
+"""
+
+rule insert_muon_recon:
+        input:
+                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_{INDEX}.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  -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule insert_muon_analysis:
+	input:
+                expand("sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV_{INDEX}.edm4hep.root",
+		    P=[50],
+		    DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
+                    INDEX=range(5),
+                ),
+                script="benchmarks/insert_muon/analysis/muon_plots.py",
+	output:
+		results_dir=directory("results/{DETECTOR_CONFIG}/insert_muon"),
+	shell:
+		"""
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""

benchmarks/insert_muon/analysis/gen_particles.cxx

+#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;
+}

benchmarks/insert_muon/analysis/muon_plots.py

+import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys
+import mplhep as hep
+hep.style.use("CMS")
+
+plt.rcParams['figure.facecolor']='white'
+plt.rcParams['savefig.facecolor']='white'
+plt.rcParams['savefig.bbox']='tight'
+
+plt.rcParams["figure.figsize"] = (7, 7)
+
+config=sys.argv[1].split("/")[1]  #results/{config}/{benchmark_name}
+outdir=sys.argv[1]+"/"
+try:
+    import os
+    os.mkdir(outdir[:-1])
+except:
+    pass
+
+def Landau(x, normalization,location,stdev):
+    #print(type(x))
+    u=(x-location)*3.591/stdev/2.355
+    renormalization = 1.64872*normalization
+    return renormalization * np.exp(-u/2 - np.exp(-u)/2)
+
+import uproot as ur
+arrays_sim={}
+momenta=50,
+for p in momenta:
+    arrays_sim[p] = ur.concatenate({
+        f'sim_output/insert_muon/{config}_sim_mu-_{p}GeV_{index}.edm4hep.root': 'events'
+        for index in range(5)
+    })
+
+for array in arrays_sim.values():
+    tilt=-0.025
+    px=array['MCParticles.momentum.x'][:,2]
+    py=array['MCParticles.momentum.y'][:,2]
+    pz=array['MCParticles.momentum.z'][:,2]
+    p=np.sqrt(px**2+py**2+pz**2)
+    
+    pxp=px*np.cos(tilt)-pz*np.sin(tilt)
+    pyp=py
+    pzp=pz*np.cos(tilt)+px*np.sin(tilt)
+    
+    array['eta_truth']=1/2*np.log((p+pzp)/(p-pzp))
+    array['phi_truth']=np.arctan2(pyp,pxp)
+    
+for p in 50,:
+    E=arrays_sim[p]["HcalEndcapPInsertHits.energy"]
+    y, x,_=plt.hist(1e3*ak.flatten(E),bins=100, range=(0, 1.2), histtype='step')
+    bc=(x[1:]+x[:-1])/2
+    from scipy.optimize import curve_fit
+    slc=abs(bc-.48)<.15
+    fnc=Landau
+    p0=[100, .5, .05]
+    #print(list(y), list(x))
+    coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+                                 sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+    print(coeff)
+    xx=np.linspace(0,.7, 100)
+    MIP=coeff[1]/1000
+    plt.plot(xx, fnc(xx,*coeff), label=f'Landau fit:\nMIP={coeff[1]*1e3:.0f}$\\pm${1e3*np.sqrt(var_matrix[1][1]):.0f} keV')
+    plt.xlabel("hit energy [MeV]")
+    plt.ylabel("hits")
+    plt.title(f"$E_{{\\mu^-}}=${p} GeV")
+    plt.legend(fontsize=20)
+    plt.savefig(outdir+"/MIP.pdf")
+
+    plt.figure(figsize=(10,7))
+    array=arrays_sim[p]
+    bins=30; r=((-np.pi, np.pi),(2.8, 4.2))
+    selection=np.sum(array["HcalEndcapPInsertHits.energy"]>0.5*MIP,axis=-1)>0
+    h1, xedges, yedges = np.histogram2d(list(array[selection]['phi_truth']),list(array[selection]['eta_truth']), bins=bins, range=r)
+    h2, xedges, yedges = np.histogram2d(list(array['phi_truth']),list(array['eta_truth']), bins=bins, range=r)
+
+    h = h1 / h2
+    pc=plt.pcolor(xedges, yedges, h.T,linewidth=0)
+    plt.xlabel("$\\phi^*$ [rad]")
+    plt.ylabel("$\\eta^*$")
+    cb = plt.colorbar(pc)
+    cb.set_label("acceptance")
+    plt.title(f"$E_{{\\mu^-}}=${p} GeV")
+    plt.savefig(outdir+"/acceptance.pdf")