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

.gitlab-ci.yml

@@ -5,6 +5,7 @@ variables:
   DETECTOR_CONFIG: epic_craterlake
   GITHUB_SHA: ''
   GITHUB_REPOSITORY: ''
+  SNAKEMAKE_FLAGS: '--cache'
 
 workflow:
   name: '$PIPELINE_NAME'

.pre-commit-config.yaml

+repos:
+-   repo: https://github.com/pre-commit/pre-commit-hooks
+    rev: v5.0.0
+    hooks:
+    -   id: check-added-large-files
+        args: ['--maxkb=128']
+    -   id: check-yaml

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 $SNAKEMAKE_FLAGS --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
@@ -19,7 +19,7 @@ bench:backgrounds_emcal_backwards:
     - ln -s $LOCAL_DATA_PATH/input input
     - export PYTHONUSERBASE=$LOCAL_DATA_PATH/deps
     - pip install -r benchmarks/backgrounds/requirements.txt
-    - snakemake --cores 8 backgrounds_ecal_backwards
+    - snakemake $SNAKEMAKE_FLAGS --cores 8 backgrounds_ecal_backwards
 
 collect_results:backgrounds:
   extends: .det_benchmark
@@ -29,5 +29,5 @@ collect_results:backgrounds:
   script:
     - ls -lrht
     - mv results{,_save}/ # move results directory out of the way to preserve it
-    - snakemake --cores 1 --delete-all-output backgrounds_ecal_backwards
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output backgrounds_ecal_backwards
     - mv results{_save,}/

benchmarks/backwards_ecal/Snakefile

@@ -4,11 +4,29 @@ def get_n_events(wildcards):
     n_events = int(n_events // (energy ** 0.5))
     return n_events
 
+import functools
+import json
+import ctypes.util
+import warnings
+from snakemake.logging import logger
+
+@functools.cache
+def get_spack_package_hash(package_name):
+    try:
+        ver_info = json.loads(subprocess.check_output(["spack", "find", "--json", package_name]))
+        return ver_info[0]["package_hash"]
+    except FileNotFoundError as e:
+        logger.warning("Spack is not installed")
+        return ""
+    except subprocess.CalledProcessError as e:
+        print(e)
+        return ""
 
 rule backwards_ecal_sim:
     input:
         steering_file=ancient("EPIC/EVGEN/SINGLE/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.steer"),
         warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+        geometry_lib=os.environ["DETECTOR_PATH"] + "/../../lib/" + ctypes.util.find_library("epic"),
     output:
         "sim_output/backwards_ecal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root",
     log:
@@ -19,7 +37,13 @@ rule backwards_ecal_sim:
         PHASE_SPACE="(3to50|45to135|130to177)deg",
         INDEX="\d{4}",
     params:
-        N_EVENTS=get_n_events
+        N_EVENTS=get_n_events,
+        SEED=lambda wildcards: "1" + wildcards.INDEX,
+        DETECTOR_PATH=os.environ["DETECTOR_PATH"],
+        DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
+        DD4HEP_HASH=get_spack_package_hash("dd4hep"),
+        NPSIM_HASH=get_spack_package_hash("npsim"),
+    cache: True
     shell:
         """
 set -m # monitor mode to prevent lingering processes
@@ -27,11 +51,11 @@ exec ddsim \
   --runType batch \
   --enableGun \
   --steeringFile "{input.steering_file}" \
-  --random.seed 1{wildcards.INDEX} \
+  --random.seed {params.SEED} \
   --filter.tracker edep0 \
   -v WARNING \
   --numberOfEvents {params.N_EVENTS} \
-  --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+  --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \
   --outputFile {output}
 """
 
@@ -45,9 +69,13 @@ rule backwards_ecal_recon:
         "sim_output/backwards_ecal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root.log",
     wildcard_constraints:
         INDEX="\d{4}",
+    params:
+        DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
+        EICRECON_HASH=get_spack_package_hash("eicrecon"),
+    cache: True
     shell: """
 set -m # monitor mode to prevent lingering processes
-exec env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \
+exec env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \
   eicrecon {input} -Ppodio:output_file={output} \
   -Ppodio:output_collections=MCParticles,EcalEndcapNRecHits,EcalEndcapNClusters
 """

benchmarks/backwards_ecal/backwards_ecal.org

@@ -113,6 +113,124 @@ for energy in energies:
     ))
 #+end_src
 
+** Energy resolution
+
+#+begin_src jupyter-python
+fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
+
+fig.suptitle(PLOT_TITLE)
+
+axs = np.ravel(np.array(axs))
+
+sigmas_rel_FWHM_cb = {}
+fractions_below = {}
+
+for ix, energy in enumerate(energies):
+    for use_clusters in [False, True]:
+        energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
+        if use_clusters:
+            clf_label = "leading cluster"
+        else:
+            clf_label = "sum all hits"
+        def clf(events):
+            if use_clusters:
+                return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
+            else:
+                return ak.sum(events["EcalEndcapNRecHits.energy"], axis=-1) / energy_value
+        e_pred = clf(e_eval[energy])
+
+        plt.sca(axs[ix])
+        counts, bins, patches = plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0.01, 1.01, 101), label=rf"$e^-$ {clf_label}", hatch=None if use_clusters else r"xxx", alpha=0.8 if use_clusters else 1.)
+        plt.title(f"{energy}")
+
+        e_over_p = (bins[1:] + bins[:-1]) / 2
+        import scipy.stats
+        def f(x, n, beta, m, loc, scale):
+            return n * scipy.stats.crystalball.pdf(x, beta, m, loc, scale)
+        p0 = (np.sum(counts[10:]), 2., 3., 0.95, 0.05)
+
+        try:
+            import scipy.optimize
+            par, pcov = scipy.optimize.curve_fit(f, e_over_p[5:], counts[5:], p0=p0, maxfev=10000)
+        except RuntimeError:
+            par = None
+        plt.plot(e_over_p, f(e_over_p, *par), label=rf"Crystal Ball fit", color="tab:green" if use_clusters else "green", lw=0.8)
+
+        def summarize_fit(par):
+            _, _, _, loc_cb, scale_cb = par
+            # Calculate FWHM
+            y_max = np.max(f(np.linspace(0., 1., 100), *par))
+            f_prime = lambda x: f(x, *par) - y_max / 2
+            x_plus, = scipy.optimize.root(f_prime, loc_cb + scale_cb).x
+            x_minus, = scipy.optimize.root(f_prime, loc_cb - scale_cb).x
+            color = "cyan" if use_clusters else "orange"
+            plt.axvline(x_minus, ls="--", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu - $FWHM")
+            plt.axvline(x_plus, ls=":", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu + $FWHM")
+            fwhm = (x_plus - x_minus) / loc_cb
+            sigma_rel_FWHM_cb = fwhm / 2 / np.sqrt(2 * np.log(2))
+
+            cutoff_x = loc_cb - fwhm
+            fraction_below = np.sum(counts[e_over_p < cutoff_x]) / ak.num(e_pred, axis=0)
+
+            return sigma_rel_FWHM_cb, fraction_below
+
+        sigma_rel_FWHM_cb, fraction_below = summarize_fit(par)
+        sigmas_rel_FWHM_cb.setdefault(clf_label, {})[energy] = sigma_rel_FWHM_cb
+        fractions_below.setdefault(clf_label, {})[energy] = fraction_below
+
+        plt.legend()
+        plt.xlabel("$E/p$", loc="right")
+        plt.ylabel("Event yield", loc="top")
+
+fig.savefig(output_dir / f"resolution_plots.pdf", bbox_inches="tight")
+fig.savefig(output_dir / f"resolution_plots.png", bbox_inches="tight")
+plt.show()
+plt.close(fig)
+
+plt.figure()
+energy_values = np.array([float(energy.replace("GeV", "").replace("MeV", "e-3")) for energy in energies])
+
+for clf_label, sigma_rel_FWHM_cb in sigmas_rel_FWHM_cb.items():
+    sigma_over_e = np.array([sigma_rel_FWHM_cb[energy] for energy in energies]) * 100 # convert to %
+
+    def f(energy, stochastic, constant):
+        return np.sqrt((stochastic / np.sqrt(energy)) ** 2 + constant ** 2)
+    cond = energy_values >= 0.5
+    try:
+        import scipy.optimize
+        par, pcov = scipy.optimize.curve_fit(f, energy_values[cond], sigma_over_e[cond], maxfev=10000)
+    except RuntimeError:
+        par = None
+    stochastic, constant = par
+
+    plt.plot(
+        energy_values,
+        sigma_over_e,
+        marker=".",
+        label=f"{clf_label}"
+    )
+    plt.plot(
+        energy_values[cond],
+        f(energy_values[cond], *par),
+        color="black",
+        ls="--",
+        lw=0.5,
+        label=f"{clf_label}, ${np.ceil(stochastic * 10) / 10:.1f}\% / \sqrt{{E}} \oplus {np.ceil(constant * 10) / 10:.1f}\%$",
+    )
+plt.plot(
+    energy_values,
+    np.sqrt((1 / energy_values) ** 2 + (1 / np.sqrt(energy_values)) ** 2 + 1 ** 2),
+    color="black", label=r"YR requirement $1\% / E \oplus 2.5\% / \sqrt{E} \oplus 1\%$",
+)
+plt.title(INPUT_PATH_FORMAT)
+plt.legend()
+plt.xlabel("Energy, GeV", loc="right")
+plt.ylabel(r"$\sigma_{E} / E$ derived from FWHM, %", loc="top")
+plt.savefig(output_dir / f"resolution.pdf", bbox_inches="tight")
+plt.savefig(output_dir / f"resolution.png", bbox_inches="tight")
+plt.show()
+#+end_src
+
 ** Pion rejection
 
 #+begin_src jupyter-python
@@ -176,10 +294,13 @@ for ix, energy in enumerate(energies):
         plt.ylabel("Pion rejection factor")
 
 fig.savefig(output_dir / f"pred.pdf", bbox_inches="tight")
+fig.savefig(output_dir / f"pred.png", bbox_inches="tight")
 plt.close(fig)
 fig_log.savefig(output_dir / f"pred_log.pdf", bbox_inches="tight")
+fig_log.savefig(output_dir / f"pred_log.png", bbox_inches="tight")
 fig_log.show()
 fig_roc.savefig(output_dir / f"roc.pdf", bbox_inches="tight")
+fig_roc.savefig(output_dir / f"roc.png", bbox_inches="tight")
 fig_roc.show()
 
 plt.figure()
@@ -196,5 +317,6 @@ plt.legend()
 plt.xlabel("Energy, GeV")
 plt.ylabel("Pion rejection at 95%")
 plt.savefig(output_dir / f"pion_rej.pdf", bbox_inches="tight")
+plt.savefig(output_dir / f"pion_rej.png", bbox_inches="tight")
 plt.show()
 #+end_src

benchmarks/backwards_ecal/config.yml

@@ -16,7 +16,7 @@ sim:backwards_ecal:
         ]
   script:
     - |
-      snakemake --cache --cores 5 sim_output/backwards_ecal/${DETECTOR_CONFIG}/${PARTICLE}/${MOMENTUM}/130to177deg/flag
+      snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/backwards_ecal/${DETECTOR_CONFIG}/${PARTICLE}/${MOMENTUM}/130to177deg/flag
 
 bench:backwards_ecal:
   extends: .det_benchmark
@@ -26,7 +26,7 @@ bench:backwards_ecal:
   script:
     - export PYTHONUSERBASE=$LOCAL_DATA_PATH/deps
     - pip install -r benchmarks/backwards_ecal/requirements.txt
-    - snakemake --cores 1 backwards_ecal
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 backwards_ecal
 
 collect_results:backwards_ecal:
   extends: .det_benchmark
@@ -36,5 +36,5 @@ collect_results:backwards_ecal:
   script:
     - ls -lrht
     - mv results{,_save}/ # move results directory out of the way to preserve it
-    - snakemake --cores 1 --delete-all-output backwards_ecal
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output backwards_ecal
     - mv results{_save,}/

benchmarks/barrel_ecal/config.yml

@@ -2,13 +2,13 @@ sim:emcal_barrel_pions:
   extends: .det_benchmark
   stage: simulate
   script:
-    - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_{piplus,piminus}_energies5.0_5.0.edm4hep.root
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_{piplus,piminus}_energies5.0_5.0.edm4hep.root
 
 sim:emcal_barrel_pi0:
   extends: .det_benchmark
   stage: simulate
   script:
-    - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_pi0_energies5.0_5.0.edm4hep.root
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_pi0_energies5.0_5.0.edm4hep.root
 
 sim:emcal_barrel_electrons:
   extends: .det_benchmark
@@ -16,20 +16,20 @@ sim:emcal_barrel_electrons:
   script:
     - if [[ "$RUN_EXTENDED_BENCHMARK" == "true" ]] ; then snakemake --cores $DETECTOR_CONFIG/results/energy_scan/emcal_barrel_electron_fsam_scan.png; fi
     - export JUGGLER_N_EVENTS=400
-    - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_electron_energies5.0_5.0.edm4hep.root
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_electron_energies5.0_5.0.edm4hep.root
 
 sim:emcal_barrel_photons:
   extends: .det_benchmark
   stage: simulate
   script:
     - if [[ "$RUN_EXTENDED_BENCHMARK" == "true" ]] ; then snakemake --cores $DETECTOR_CONFIG/results/energy_scan/emcal_barrel_proton_fsam_scan.png; fi
-    - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_photon_energies5.0_5.0.edm4hep.root
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_photon_energies5.0_5.0.edm4hep.root
 
 sim:emcal_barrel_pion_rejection:
   extends: .det_benchmark
   stage: simulate
   script:
-    - snakemake --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_{piminus,electron}_energies1.0_18.0.edm4hep.root
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_{piminus,electron}_energies1.0_18.0.edm4hep.root
 
 calib:emcal_barrel_electrons:
   extends: .det_benchmark
@@ -39,7 +39,7 @@ calib:emcal_barrel_electrons:
   script:
     - ls -lhtR sim_output/
     - rootls -t sim_output/sim_emcal_barrel_electron_energies5.0_5.0.edm4hep.root
-    - snakemake --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_electron_calibration.json
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_electron_calibration.json
     - mv sim_output/sim_emcal_barrel_electron.edm4hep.root results/.
     - echo "JSON file(s) from analysis:" ; cat results/*.json
 
@@ -49,7 +49,7 @@ bench:emcal_barrel_pions:
   needs:
     - ["sim:emcal_barrel_pions"]
   script:
-    - snakemake --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_pions_Ethr.png
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_pions_Ethr.png
 
 bench:emcal_barrel_electrons_scan:
   extends: .det_benchmark
@@ -66,7 +66,7 @@ bench:emcal_barrel_pi0:
     - ["sim:emcal_barrel_pi0", "calib:emcal_barrel_electrons"]
   script:
     - echo "JSON file(s) from analysis:" ; cat results/*.json
-    - snakemake --cores 1 epic_craterlake/results/Barrel_emcal_pi0.json
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 epic_craterlake/results/Barrel_emcal_pi0.json
 
 bench:emcal_barrel_photons:
   extends: .det_benchmark
@@ -76,7 +76,7 @@ bench:emcal_barrel_photons:
   script:
     - ls -lhtR sim_output/
     - rootls -t sim_output/sim_emcal_barrel_photon_energies5.0_5.0.edm4hep.root
-    - snakemake --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_photon_calibration.json
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/results/emcal_barrel_photon_calibration.json
     - mv sim_output/sim_emcal_barrel_photon.edm4hep.root results/.
     - if [[ "$RUN_EXTENDED_BENCHMARK" == "true" ]] ; then snakemake --cores $DETECTOR_CONFIG/results/energy_scan/emcal_barrel_proton_fsam_scan.png; fi
 
@@ -89,7 +89,7 @@ bench:emcal_barrel_pion_rejection:
     - ls -lhtR sim_output/
     - rootls -t $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_piminus_energies1.0_18.0.edm4hep.root
     - rootls -t $DETECTOR_CONFIG/sim_output/sim_emcal_barrel_electron_energies1.0_18.0.edm4hep.root
-    - snakemake --cores 1 $DETECTOR_CONFIG/results/Barrel_emcal_pion_rej.json
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 $DETECTOR_CONFIG/results/Barrel_emcal_pion_rej.json
 
 collect_results:barrel_ecal:
   extends: .det_benchmark

benchmarks/ecal_gaps/config.yml

@@ -4,7 +4,7 @@ sim:ecal_gaps:
   script:
     - mkdir -p $LOCAL_DATA_PATH/input
     - ln -s $LOCAL_DATA_PATH/input input
-    - snakemake --cache --cores 10 ecal_gaps --omit-from ecal_gaps
+    - snakemake $SNAKEMAKE_FLAGS --cache --cores 10 ecal_gaps --omit-from ecal_gaps
 
 bench:ecal_gaps:
   extends: .det_benchmark
@@ -15,7 +15,7 @@ bench:ecal_gaps:
     - ln -s $LOCAL_DATA_PATH/input input
     - export PYTHONUSERBASE=$LOCAL_DATA_PATH/deps
     - pip install -r benchmarks/ecal_gaps/requirements.txt
-    - snakemake --cores 8 ecal_gaps
+    - snakemake $SNAKEMAKE_FLAGS --cores 8 ecal_gaps
 
 collect_results:ecal_gaps:
   extends: .det_benchmark
@@ -25,5 +25,5 @@ collect_results:ecal_gaps:
   script:
     - ls -lrht
     - mv results{,_save}/ # move results directory out of the way to preserve it
-    - snakemake --cores 1 --delete-all-output ecal_gaps
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output ecal_gaps
     - mv results{_save,}/

benchmarks/femc_electron/Snakefile

+def get_n_events(wildcards):
+    energy = float(wildcards.P)
+    n_events = 1000
+    n_events = int(n_events // ((energy / 20) ** 0.5))
+    return n_events
+
+
 rule femc_electron_generate:
     input:
         script="benchmarks/femc_electron/analysis/gen_particles.cxx",
     params:
-		NEVENTS_GEN=1000,
+        N_EVENTS=get_n_events,
         th_max=28,
         th_min=2.0
     output:
@@ -10,23 +17,23 @@ rule femc_electron_generate:
     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})'
+root -l -b -q '{input.script}({params.N_EVENTS},"{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:
+        N_EVENTS=get_n_events,
         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 \
+   --numberOfEvents {params.N_EVENTS} \
    --physicsList {params.PHYSICS_LIST} \
    --inputFiles {input.GEN_FILE} \
    --outputFile {output.SIM_FILE}
@@ -37,10 +44,11 @@ rule femc_electron_recon:
         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"
+    params:
+        N_EVENTS=get_n_events,
     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
+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={params.N_EVENTS}
 """
 
 rule femc_electron_analysis:

benchmarks/femc_electron/analysis/femc_electron_plots.py

@@ -132,7 +132,7 @@ for p in arrays_sim:
     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))
+    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
     #res=np.abs(coeff[2]/coeff[1])
     if p==50:
         xx=np.linspace(15*p/20,22*p/20, 100)
@@ -151,7 +151,7 @@ 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)
+coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0, sigma=dres, maxfev=10000)
 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()
@@ -190,7 +190,7 @@ for eta_min, eta_max in zip(partitions[:-1], partitions[1:]):
         p0=(100, p, 3)
 
         try:
-            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
             if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
                 continue
             pvals.append(p)

benchmarks/femc_electron/config.yml

@@ -13,7 +13,7 @@ sim:femc_electron:
       - P: 80
   timeout: 1 hours
   script:
-    - snakemake --cores 1 sim_output/femc_electron/epic_craterlake_rec_e-_${P}GeV.edm4eic.root
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 sim_output/femc_electron/epic_craterlake_rec_e-_${P}GeV.edm4eic.root
   retry:
     max: 2
     when:
@@ -24,7 +24,7 @@ bench:femc_electron:
   stage: benchmarks
   needs: ["sim:femc_electron"]
   script:
-    - snakemake --cores 1 results/epic_craterlake/femc_electron
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/femc_electron
 
 collect_results:femc_electron:
   extends: .det_benchmark
@@ -33,5 +33,5 @@ collect_results: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
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/femc_electron
     - mv results{_save,}/

benchmarks/femc_photon/Snakefile

+def get_n_events(wildcards):
+    energy = float(wildcards.P)
+    n_events = 1000
+    n_events = int(n_events // ((energy / 20) ** 0.5))
+    return n_events
+
+
 rule femc_photon_generate:
     input:
         script="benchmarks/femc_photon/analysis/gen_particles.cxx",
     params:
-		NEVENTS_GEN=1000,
+        N_EVENTS=get_n_events,
         th_max=28,
         th_min=2.0
     output:
@@ -10,23 +17,23 @@ rule femc_photon_generate:
     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})'
+root -l -b -q '{input.script}({params.N_EVENTS},"{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:
+        N_EVENTS=get_n_events,
         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 \
+   --numberOfEvents {params.N_EVENTS} \
    --physicsList {params.PHYSICS_LIST} \
    --inputFiles {input.GEN_FILE} \
    --outputFile {output.SIM_FILE}
@@ -37,10 +44,11 @@ rule femc_photon_recon:
         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"
+    params:
+        N_EVENTS=get_n_events,
     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
+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={params.N_EVENTS}
 """
 
 rule femc_photon_analysis:

benchmarks/femc_photon/analysis/femc_photon_plots.py

@@ -131,7 +131,7 @@ for p in arrays_sim:
     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))
+    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
     #res=np.abs(coeff[2]/coeff[1])
     if p==50:
         xx=np.linspace(15*p/20,22*p/20, 100)
@@ -150,7 +150,7 @@ 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)
+coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0, sigma=dres, maxfev=10000)
 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()
@@ -188,7 +188,7 @@ for eta_min, eta_max in zip(partitions[:-1], partitions[1:]):
         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))
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
             if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
                 continue
             pvals.append(p)

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/config.yml

@@ -13,7 +13,7 @@ sim:femc_photon:
       - P: 80
   timeout: 1 hours
   script:
-    - snakemake --cores 1 sim_output/femc_photon/epic_craterlake_rec_photon_${P}GeV.edm4eic.root
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 sim_output/femc_photon/epic_craterlake_rec_photon_${P}GeV.edm4eic.root
   retry:
     max: 2
     when:
@@ -24,7 +24,7 @@ bench:femc_photon:
   stage: benchmarks
   needs: ["sim:femc_photon"]
   script:
-    - snakemake --cores 1 results/epic_craterlake/femc_photon
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/femc_photon
 
 collect_results:femc_photon:
   extends: .det_benchmark
@@ -33,5 +33,5 @@ collect_results: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
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/femc_photon
     - mv results{_save,}/

benchmarks/femc_pi0/Snakefile

+def get_n_events(wildcards):
+    energy = float(wildcards.P)
+    n_events = 1000
+    n_events = int(n_events // ((energy / 20) ** 0.5))
+    return n_events
+
+
 rule femc_pi0_generate:
     input:
         script="benchmarks/femc_pi0/analysis/gen_particles.cxx",
     params:
-		NEVENTS_GEN=1000,
+        N_EVENTS=get_n_events,
         th_max=28,
         th_min=2.0
     output:
@@ -10,23 +17,23 @@ rule femc_pi0_generate:
     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})'
+root -l -b -q '{input.script}({params.N_EVENTS},"{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:
+        N_EVENTS=get_n_events,
         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 \
+   --numberOfEvents {params.N_EVENTS} \
    --physicsList {params.PHYSICS_LIST} \
    --inputFiles {input.GEN_FILE} \
    --outputFile {output.SIM_FILE}
@@ -37,10 +44,11 @@ rule femc_pi0_recon:
         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"
+    params:
+        N_EVENTS=get_n_events,
     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
+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={params.N_EVENTS}
 """
 
 rule femc_pi0_analysis:

benchmarks/femc_pi0/analysis/femc_pi0_plots.py

@@ -167,7 +167,7 @@ for p in arrays_sim:
     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))
+    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
     #res=np.abs(coeff[2]/coeff[1])
     if p==50:
         xx=np.linspace(15*p/20,22*p/20, 100)
@@ -186,7 +186,7 @@ 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)
+coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0, sigma=dres, maxfev=10000)
 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()
@@ -225,7 +225,7 @@ for eta_min, eta_max in zip(partitions[:-1], partitions[1:]):
         p0=(100, p, 3)
 
         try:
-            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
             if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
                 continue
             pvals.append(p)

benchmarks/femc_pi0/config.yml

@@ -12,7 +12,7 @@ sim:femc_pi0:
       - P: 70
       - P: 80
   script:
-    - snakemake --cores 1 sim_output/femc_pi0/epic_craterlake_rec_pi0_${P}GeV.edm4eic.root
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 sim_output/femc_pi0/epic_craterlake_rec_pi0_${P}GeV.edm4eic.root
   retry:
     max: 2
     when:
@@ -23,7 +23,7 @@ bench:femc_pi0:
   stage: benchmarks
   needs: ["sim:femc_pi0"]
   script:
-    - snakemake --cores 1 results/epic_craterlake/femc_pi0
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/femc_pi0
 
 collect_results:femc_pi0:
   extends: .det_benchmark
@@ -32,5 +32,5 @@ collect_results: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
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/femc_pi0
     - mv results{_save,}/

benchmarks/insert_muon/analysis/muon_plots.py

@@ -55,7 +55,7 @@ for p in 50,:
     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)))
+                                 sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000)
     print(coeff)
     xx=np.linspace(0,.7, 100)
     MIP=coeff[1]/1000

benchmarks/insert_muon/config.yml

@@ -5,7 +5,7 @@ sim:insert_muon:
     matrix:
       - P: 50
   script:
-    - snakemake --cores 5 sim_output/insert_muon/epic_craterlake_sim_mu-_${P}GeV_{0,1,2,3,4}.edm4hep.root
+    - snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/insert_muon/epic_craterlake_sim_mu-_${P}GeV_{0,1,2,3,4}.edm4hep.root
   retry:
     max: 2
     when:
@@ -16,7 +16,7 @@ bench:insert_muon:
   stage: benchmarks
   needs: ["sim:insert_muon"]
   script:
-    - snakemake --cores 1 results/epic_craterlake/insert_muon
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/insert_muon
 
 collect_results:insert_muon:
   extends: .det_benchmark
@@ -25,5 +25,5 @@ collect_results:insert_muon:
   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/insert_muon
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/insert_muon
     - mv results{_save,}/