Merge branch 'master' into tracking_benchmark

.github/workflows/mirror.yaml

@@ -55,5 +55,5 @@ jobs:
           -H "Accept: application/vnd.github+json" \
           -H "Authorization: Bearer $GITHUB_TOKEN" \
           -H "X-GitHub-Api-Version: 2022-11-28" \
-          -d '{"context": "eicweb/detector_benchmarks ('"$DETECTOR_CONFIG"')", "state": "pending", "description": "Waiting for response from the EICweb", "target_url": "${{ fromJson(steps.trigger_eicweb.outputs.json).web_url }}"}' \
+          -d '{"context": "eicweb/detector_benchmarks (nightly, '"$DETECTOR_CONFIG"')", "state": "pending", "description": "Waiting for response from the EICweb", "target_url": "${{ fromJson(steps.trigger_eicweb.outputs.json).web_url }}"}' \
           "https://api.github.com/repos/${{ github.repository }}/statuses/${{ github.event.pull_request.head.sha || github.sha }}"

.gitlab-ci.yml

@@ -35,6 +35,7 @@ default:
       - summary.txt
     reports:
       dotenv: .env
+    when: always
 
 stages:
   - status-pending
@@ -60,8 +61,8 @@ stages:
           "https://api.github.com/repos/${GITHUB_REPOSITORY}/statuses/${GITHUB_SHA}" \
           -d '{"state":"'"${STATE}"'",
                "target_url":"'"${CI_PIPELINE_URL}"'",
-               "description":"'"${DESCRIPTION} $(TZ=America/New_York date)"'",
-               "context":"eicweb/detector_benchmarks ('"$DETECTOR_CONFIG"')"
+               "description":"'"$(TZ=America/New_York date)"'",
+               "context":"eicweb/detector_benchmarks ('"${BENCHMARKS_TAG}"', '"$DETECTOR_CONFIG"')"
               }' ;
       fi
 
@@ -71,7 +72,6 @@ benchmarks:detector:pending:
   extends: .status
   variables:
     STATE: "pending"
-    DESCRIPTION: "Started..."
   when: always
 
 common:setup:
@@ -184,8 +184,11 @@ deploy_results:
   script:
     - snakemake $SNAKEMAKE_FLAGS --cores 1 results/metadata.json
     - 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
+    - wget https://dl.pelicanplatform.org/7.13.0/pelican_Linux_x86_64.tar.gz
+    - sha256sum -c <(echo '38ac8548c67302299e50a1b81c159ed418e90d84a6606ddd377fd2c8b164d114  pelican_Linux_x86_64.tar.gz')
+    - tar zxf pelican_Linux_x86_64.tar.gz
+    - mv results pipeline-$CI_PIPELINE_ID; tar cf pipeline-$CI_PIPELINE_ID.tar pipeline-$CI_PIPELINE_ID/; mv pipeline-$CI_PIPELINE_ID results
+    - ./pelican-*/pelican object copy pipeline-$CI_PIPELINE_ID.tar $OSDF_ENDPOINT$OSDF_OUTPUT_PREFIX/pipeline-$CI_PIPELINE_ID.tar
 
 benchmarks:detector:success:
   stage: status-report
@@ -193,7 +196,6 @@ benchmarks:detector:success:
   extends: .status
   variables:
     STATE: "success"
-    DESCRIPTION: "Succeeded!"
   after_script:
     # Cleanup scratch space
     - rm -rfv $LOCAL_DATA_PATH
@@ -205,7 +207,6 @@ benchmarks:detector:failure:
   extends: .status
   variables:
     STATE: "failure"
-    DESCRIPTION: "Failed!"
   when: on_failure
 
 pages:

Snakefile

@@ -60,7 +60,7 @@ def get_remote_path(path):
     if use_s3:
         return f"s3https://eics3.sdcc.bnl.gov:9000/eictest/{path}"
     elif use_xrootd:
-        return f"root://dtn-eic.jlab.org//work/eic2/{path}"
+        return f"root://dtn-eic.jlab.org//volatile/eic/{path}"
     else:
         raise runtime_exception('Unexpected value for config["remote_provider"]: {config["remote_provider"]}')
 
@@ -74,7 +74,7 @@ rule fetch_epic:
     cache: True
     retries: 3
     shell: """
-xrdcp --debug 2 root://dtn-eic.jlab.org//work/eic2/{output.filepath} {output.filepath}
+xrdcp --debug 2 root://dtn-eic.jlab.org//volatile/eic/{output.filepath} {output.filepath}
 """ if use_xrootd else """
 mc cp S3/eictest/{output.filepath} {output.filepath}
 """ if use_s3 else f"""
@@ -132,6 +132,9 @@ cat > {output} <<EOF
   "CI_PIPELINE_ID": "${{CI_PIPELINE_ID:-}}",
   "CI_PIPELINE_SOURCE": "${{CI_PIPELINE_SOURCE:-}}",
   "CI_PROJECT_ID": "${{CI_PROJECT_ID:-}}",
+  "GITHUB_REPOSITORY": "${{GITHUB_REPOSITORY:-}}",
+  "GITHUB_SHA": "${{GITHUB_SHA:-}}",
+  "GITHUB_PR": "${{GITHUB_PR:-}}",
   "PIPELINE_NAME": "${{PIPELINE_NAME:-}}"
 }}
 EOF

benchmarks/backgrounds/Snakefile

@@ -38,6 +38,9 @@ rule backgrounds_ecal_backwards:
         proton_beam_gas_sim="sim_output/" + DETECTOR_CONFIG + "/backgrounds/EPIC/EVGEN/BACKGROUNDS/BEAMGAS/proton/pythia8.306-1.0/100GeV/pythia8.306-1.0_ProtonBeamGas_100GeV_run001.edm4hep.root",
     output:
         directory("results/backgrounds/backwards_ecal")
+    log:
+        scheduler=".logs/results/backgrounds/backwards_ecal/scheduler.log",
+        worker=".logs/results/backgrounds/backwards_ecal/worker.log",
     threads: workflow.cores
     shell:
         """
@@ -49,17 +52,17 @@ cleanup() {{
 trap cleanup EXIT
 
 PORT=$RANDOM
-dask scheduler --port $PORT &
+dask scheduler --port $PORT 2>{log.scheduler} &
 export DASK_SCHEDULER=localhost:$PORT
 SCHEDULER_PID=$!
-dask worker tcp://$DASK_SCHEDULER --nworkers {threads} --nthreads 1 &
+dask worker tcp://$DASK_SCHEDULER --nworkers {threads} --nthreads 1 2>{log.worker} &
 WORKER_PID=$!
 env \
 MATPLOTLIBRC={input.matplotlibrc} \
-ELECTRON_BEAM_GAS_GEN=root://dtn-eic.jlab.org//work/eic2/EPIC/EVGEN/BACKGROUNDS/BEAMGAS/electron/GETaLM1.0.0-1.0/10GeV/GETaLM1.0.0-1.0_ElectronBeamGas_10GeV_foam_emin10keV_run001.hepmc3.tree.root \
+ELECTRON_BEAM_GAS_GEN=root://dtn-eic.jlab.org//volatile/eic/EPIC/EVGEN/BACKGROUNDS/BEAMGAS/electron/GETaLM1.0.0-1.0/10GeV/GETaLM1.0.0-1.0_ElectronBeamGas_10GeV_foam_emin10keV_run001.hepmc3.tree.root \
 ELECTRON_BEAM_GAS_SIM=$(realpath {input.electron_beam_gas_sim}) \
 PHYSICS_PROCESS_SIM=$(realpath {input.physics_signal_sim}) \
-PROTON_BEAM_GAS_GEN=root://dtn-eic.jlab.org//work/eic2/EPIC/EVGEN/BACKGROUNDS/BEAMGAS/proton/pythia8.306-1.0/100GeV/pythia8.306-1.0_ProtonBeamGas_100GeV_run001.hepmc3.tree.root \
+PROTON_BEAM_GAS_GEN=root://dtn-eic.jlab.org//volatile/eic/EPIC/EVGEN/BACKGROUNDS/BEAMGAS/proton/pythia8.306-1.0/100GeV/pythia8.306-1.0_ProtonBeamGas_100GeV_run001.hepmc3.tree.root \
 PROTON_BEAM_GAS_SIM=$(realpath {input.proton_beam_gas_sim}) \
 OUTPUT_DIR={output} \
 python {input.script}

benchmarks/backgrounds/requirements.txt

-awkward >= 2.4.0
-dask >= 2023
-dask_awkward
-dask_histogram
-distributed >= 2023
+awkward ~= 2.7.2
+dask ~= 2024.12.1
+dask_awkward ~= 2024.12.2
+dask_histogram ~= 2024.12.1
+distributed ~= 2024.12.1
 pyhepmc
-uproot >= 5.2.0
+uproot ~= 5.5.1

benchmarks/backwards_ecal/Snakefile

@@ -18,7 +18,7 @@ rule backwards_ecal_sim:
         PARTICLE="(e-|pi-)",
         ENERGY="[0-9]+[kMG]eV",
         PHASE_SPACE="(3to50|45to135|130to177)deg",
-        INDEX="\d{4}",
+        INDEX=r"\d{4}",
     params:
         N_EVENTS=get_n_events,
         SEED=lambda wildcards: "1" + wildcards.INDEX,
@@ -51,7 +51,7 @@ rule backwards_ecal_recon:
     log:
         "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}",
+        INDEX=r"\d{4}",
     params:
         DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
         EICRECON_HASH=get_spack_package_hash("eicrecon"),
@@ -64,14 +64,69 @@ exec env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \
 """
 
 
-rule backwards_ecal_recon_many:
+rule backwards_ecal_local_sim_list:
     input:
         expand(
             "sim_output/backwards_ecal/{{DETECTOR_CONFIG}}/{{PARTICLE}}/{{ENERGY}}/{{PHASE_SPACE}}/{{PARTICLE}}_{{ENERGY}}_{{PHASE_SPACE}}.{INDEX:04d}.eicrecon.tree.edm4eic.root",
             INDEX=range(20),
         ),
     output:
-        touch("sim_output/backwards_ecal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/flag"),
+        "listing/backwards_ecal/local/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst",
+    run:
+        with open(output[0], "wt") as fp:
+            fp.write("\n".join(input))
+
+
+if config.get("stream_from_xrootd", True) not in [False, "", "0", "false"]:
+    rule backwards_ecal_campaign_sim_list:
+        output:
+            "listing/backwards_ecal/{CAMPAIGN}/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst",
+        params:
+            search_path=lambda wildcards: f"EPIC/RECO/{wildcards.CAMPAIGN}/epic_craterlake/SINGLE/{wildcards.PARTICLE}/{wildcards.ENERGY}/{wildcards.PHASE_SPACE}/",
+        shell: """
+    xrdfs root://dtn-eic.jlab.org/ ls /volatile/eic/{params.search_path} \
+      | awk '{{ print "root://dtn-eic.jlab.org/"$1; }}' \
+      | sort \
+      > {output}
+    if [ ! -s {output} ]; then
+      echo "Got an empty file listing for path \"\""
+      exit 1
+    fi
+    """
+else:
+    checkpoint backwards_ecal_campaign_sim_list_checkpoint:
+        output:
+            "listing/backwards_ecal/{CAMPAIGN}/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst.orig",
+        params:
+            search_path=lambda wildcards: f"EPIC/RECO/{wildcards.CAMPAIGN}/epic_craterlake/SINGLE/{wildcards.PARTICLE}/{wildcards.ENERGY}/{wildcards.PHASE_SPACE}/",
+        shell: """
+    xrdfs root://dtn-eic.jlab.org/ ls /volatile/eic/{params.search_path} \
+      | sed -e 's#^/volatile/eic/##' \
+      | sort \
+      > {output}
+    if [ ! -s {output} ]; then
+      echo "Got an empty file listing for path \"\""
+      exit 1
+    fi
+    """
+
+    def get_backwards_ecal_campaign_sim_list(wildcards):
+        with checkpoints.backwards_ecal_campaign_sim_list_checkpoint.get(**wildcards).output[0].open() as fp:
+            return [line.rstrip() for line in fp.readlines()]
+
+    rule backwards_ecal_campaign_sim_list:
+        input:
+            # depend on paths from the file list
+            get_backwards_ecal_campaign_sim_list,
+            orig_list="listing/backwards_ecal/{CAMPAIGN}/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst.orig",
+        output:
+            "listing/backwards_ecal/{CAMPAIGN}/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst",
+        shell: """
+    cp {input.orig_list} {output}
+    """
+
+
+ruleorder: backwards_ecal_local_sim_list > backwards_ecal_campaign_sim_list
 
 
 DETECTOR_CONFIG=os.environ["DETECTOR_CONFIG"]
@@ -79,7 +134,7 @@ DETECTOR_CONFIG=os.environ["DETECTOR_CONFIG"]
 rule backwards_ecal:
     input:
         expand(
-            "sim_output/backwards_ecal/" + DETECTOR_CONFIG + "/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/flag",
+            "listing/backwards_ecal/{{CAMPAIGN}}/" + DETECTOR_CONFIG + "/{PARTICLE}/{ENERGY}/{PHASE_SPACE}.lst",
             PARTICLE=["pi-", "e-"],
             ENERGY=[
                 "100MeV",
@@ -96,14 +151,36 @@ rule backwards_ecal:
         matplotlibrc=".matplotlibrc",
         script="benchmarks/backwards_ecal/backwards_ecal.py",
     output:
-        directory("results/backwards_ecal")
+        directory("results/backwards_ecal/{CAMPAIGN}/")
+    log:
+        scheduler=".logs/results/backwards_ecal/{CAMPAIGN}/scheduler.log",
+        worker=".logs/results/backwards_ecal/{CAMPAIGN}/worker.log",
+    threads: workflow.cores
     shell:
         """
+if [[ "{wildcards.CAMPAIGN}" == "local" ]]; then
+export PLOT_TITLE="Benchmark simulation"
+else
+export PLOT_TITLE="\\textbf{{ePIC}} Simulation {wildcards.CAMPAIGN}"
+fi
+
+set -m # monitor mode to prevent lingering processes
+cleanup() {{
+  echo Cleaning up
+  kill $WORKER_PID $SCHEDULER_PID
+}}
+trap cleanup EXIT
+
+PORT=$RANDOM
+dask scheduler --port $PORT 2>{log.scheduler} &
+export DASK_SCHEDULER=localhost:$PORT
+SCHEDULER_PID=$!
+dask worker tcp://$DASK_SCHEDULER --nworkers {threads} --nthreads 1 2>{log.worker} &
+WORKER_PID=$!
 env \
 MATPLOTLIBRC={input.matplotlibrc} \
 DETECTOR_CONFIG=""" + DETECTOR_CONFIG + """ \
-PLOT_TITLE=""" + DETECTOR_CONFIG + """ \
-INPUT_PATH_FORMAT=sim_output/backwards_ecal/""" + DETECTOR_CONFIG + """/{{particle}}/{{energy}}/130to177deg/{{particle}}_{{energy}}_130to177deg.{{ix:04d}}.eicrecon.tree.edm4eic.root \
+INPUT_PATH_FORMAT=listing/backwards_ecal/{wildcards.CAMPAIGN}/""" + DETECTOR_CONFIG + """/{{particle}}/{{energy}}/130to177deg.lst \
 OUTPUT_DIR={output} \
 python {input.script}
 """

benchmarks/backwards_ecal/backwards_ecal.org

@@ -18,12 +18,32 @@ import os
 from pathlib import Path
 
 import awkward as ak
+import boost_histogram as bh
+import dask_histogram as dh
 import numpy as np
-import vector
 import uproot
-from sklearn.metrics import roc_curve
+import vector
+#+end_src
+
+#+begin_src jupyter-python :results silent
+from dask.distributed import Client
 
+client = Client(os.environ.get("DASK_SCHEDULER", "localhost:8786"))
+#+end_src
+
+#+begin_src jupyter-python
 vector.register_awkward()
+
+from dask.distributed import WorkerPlugin
+class VectorImporter(WorkerPlugin):
+    idempotent=True
+
+    def setup(self, worker):
+        import vector
+
+        vector.register_awkward()
+
+client.register_plugin(VectorImporter())
 #+end_src
 
 * Plotting setup
@@ -57,7 +77,6 @@ setup_presentation_style()
 DETECTOR_CONFIG=os.environ.get("DETECTOR_CONFIG")
 PLOT_TITLE=os.environ.get("PLOT_TITLE")
 INPUT_PATH_FORMAT=os.environ.get("INPUT_PATH_FORMAT", "EPIC/RECO/24.04.0/epic_craterlake/SINGLE/{particle}/{energy}/130to177deg/{particle}_{energy}_130to177deg.{ix:04d}.eicrecon.tree.edm4eic.root")
-INPUT_PATH_INDEX_RANGE=list(range(20))
 
 output_dir=Path(os.environ.get("OUTPUT_DIR", "./"))
 output_dir.mkdir(parents=True, exist_ok=True)
@@ -92,95 +111,114 @@ energies = [
     "10GeV",
     "20GeV",
 ]
-
 filter_name = [
     "MCParticles.*",
-    "*EcalEndcapNRecHits*",
-    "*EcalEndcapNClusters*",
+    "EcalEndcapNClusters.energy",
 ]
 
 pi_eval = {}
 e_eval = {}
 
+def readlist(path):
+    with open(path, "rt") as fp:
+        paths = [line.rstrip() for line in fp.readlines()]
+    return paths
+
 for energy in energies:
-    pi_eval[energy] = filter_pointing(uproot.concatenate(
-        {INPUT_PATH_FORMAT.format(particle="pi-", energy=energy, ix=ix): "events" for ix in INPUT_PATH_INDEX_RANGE},
+    pi_eval[energy] = uproot.dask(
+        {path: "events" for path in readlist(INPUT_PATH_FORMAT.format(particle="pi-", energy=energy))},
         filter_name=filter_name,
-    ))
-    e_eval[energy] = filter_pointing(uproot.concatenate(
-        {INPUT_PATH_FORMAT.format(particle="e-", energy=energy, ix=ix): "events" for ix in INPUT_PATH_INDEX_RANGE},
+        steps_per_file=1,
+        open_files=False,
+    ).map_partitions(filter_pointing)
+    e_eval[energy] = uproot.dask(
+        {path: "events" for path in readlist(INPUT_PATH_FORMAT.format(particle="e-", energy=energy))},
         filter_name=filter_name,
-    ))
+        steps_per_file=1,
+        open_files=False,
+    ).map_partitions(filter_pointing)
 #+end_src
 
 ** Energy resolution
 
+The first step in the analysis is plotting distributions of cluster energies.
+Only electron distribution is needed, but we also do pi- for the pion rejection
+study.
+
 #+begin_src jupyter-python
-fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
+e_over_p_hist = {}
+e_over_p_axis = bh.axis.Regular(101, 0., 1.01)
 
-fig.suptitle(PLOT_TITLE)
+for ix, energy in enumerate(energies):
+    for particle_name, dataset in [("e-", e_eval[energy]), ("pi-", pi_eval[energy])]:
+        energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
+        def clf(events):
+            return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
+        e_over_p_hist[(particle_name, energy)] = dh.factory(
+            clf(dataset),
+            axes=(e_over_p_axis,),
+        )
+
+e_over_p_hist = client.gather(client.compute(e_over_p_hist))
+#+end_src
+
+The next step is to plot the histograms together with fits and produce a plot
+for resolution vs energy that summarizes them.
+
+#+begin_src jupyter-python
+fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
 
 axs = np.ravel(np.array(axs))
 
 sigmas_rel_FWHM_cb = {}
 fractions_below = {}
 
+def norm_by_sum(arr):
+    return arr / np.sum(arr)
+
 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}")
+    plt.sca(axs[ix])
+    hist = e_over_p_hist[("e-", energy)]
+    patch = plt.stairs(norm_by_sum(hist.values()), hist.axes[0].edges, label=rf"$e^-$ {PLOT_TITLE}")
+    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
+    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(hist.values()[10:]), 2., 3., 0.95, 0.05)
 
-        plt.legend()
-        plt.xlabel("$E/p$", loc="right")
-        plt.ylabel("Event yield", loc="top")
+    try:
+        import scipy.optimize
+        par, pcov = scipy.optimize.curve_fit(f, hist.axes[0].centers[5:], hist.values()[5:], p0=p0, maxfev=10000)
+    except RuntimeError:
+        print(hist)
+        raise
+    plt.plot(hist.axes[0].centers, f(hist.axes[0].centers, *par), label=rf"Crystal Ball fit", color="tab: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
+        plt.axvline(x_minus, ls="--", lw=0.75, color=patch.get_facecolor(), label=r"$\mu - $FWHM")
+        plt.axvline(x_plus, ls=":", lw=0.75, color=patch.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(hist.values()[hist.axes[0].centers < cutoff_x]) / np.sum(hist.values())
+
+        return sigma_rel_FWHM_cb, fraction_below
+
+    sigma_rel_FWHM_cb, fraction_below = summarize_fit(par)
+    sigmas_rel_FWHM_cb.setdefault(PLOT_TITLE, {})[energy] = sigma_rel_FWHM_cb
+    fractions_below.setdefault(PLOT_TITLE, {})[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")
@@ -200,29 +238,39 @@ for clf_label, sigma_rel_FWHM_cb in sigmas_rel_FWHM_cb.items():
         import scipy.optimize
         par, pcov = scipy.optimize.curve_fit(f, energy_values[cond], sigma_over_e[cond], maxfev=10000)
     except RuntimeError:
-        par = None
+        print("energy_values", energy_values[cond])
+        print("sigma_over_e", sigma_over_e[cond])
+        raise
     stochastic, constant = par
 
     plt.plot(
         energy_values,
         sigma_over_e,
         marker=".",
+        ls="none",
         label=f"{clf_label}"
     )
+    xmin = np.min(energy_values[cond])
+    xmax = np.max(energy_values[cond])
+    xs = np.arange(xmin, xmax, 0.1)
     plt.plot(
-        energy_values[cond],
-        f(energy_values[cond], *par),
-        color="black",
+        xs,
+        f(xs, *par),
         ls="--",
         lw=0.5,
-        label=f"{clf_label}, ${np.ceil(stochastic * 10) / 10:.1f}\% / \sqrt{{E}} \oplus {np.ceil(constant * 10) / 10:.1f}\%$",
+        label=fr"Functional fit: ${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\%$",
+xmin = np.min(energy_values)
+xmax = np.max(energy_values) * 1.05
+xs = np.arange(xmin, xmax + 0.1, 0.1)
+plt.fill_between(
+    xs,
+    np.sqrt((2 / np.sqrt(xs)) ** 2 + 1 ** 2),
+    np.sqrt((2 / np.sqrt(xs)) ** 2 + 3 ** 2),
+    lw=0., alpha=0.2, color="black", label=r"YR requirement $2\% / \sqrt{E} \oplus (1-3)\%$",
 )
-plt.title(INPUT_PATH_FORMAT)
+plt.xlim(0., xmax)
+plt.ylim(top=6.)
 plt.legend()
 plt.xlabel("Energy, GeV", loc="right")
 plt.ylabel(r"$\sigma_{E} / E$ derived from FWHM, %", loc="top")
@@ -238,10 +286,6 @@ fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
 fig_log, axs_log = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
 fig_roc, axs_roc = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
 
-fig.suptitle(PLOT_TITLE)
-fig_log.suptitle(PLOT_TITLE)
-fig_roc.suptitle(PLOT_TITLE)
-
 axs = np.ravel(np.array(axs))
 axs_log = np.ravel(np.array(axs_log))
 axs_roc = np.ravel(np.array(axs_roc))
@@ -249,49 +293,41 @@ axs_roc = np.ravel(np.array(axs_roc))
 rocs = {}
 
 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])
-        pi_pred = clf(pi_eval[energy])
-
-        for do_log, ax in [(False, axs[ix]), (True, axs_log[ix])]:
-            plt.sca(ax)
-            plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0., 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.hist(pi_pred, weights=np.full_like(pi_pred, 1.0 / ak.num(pi_pred, axis=0)), bins=np.linspace(0., 1.01, 101), label=rf"$\pi^-$ {clf_label}", histtype="step")
-            plt.title(f"{energy}")
-            plt.legend()
-            plt.xlabel("Classifier output")
-            plt.ylabel("Event yield")
-            if do_log:
-                plt.yscale("log")
-
-        plt.sca(axs_roc[ix])
-        fpr, tpr, _ = roc_curve(
-            np.concatenate([np.ones_like(e_pred), np.zeros_like(pi_pred)]),
-            np.concatenate([e_pred, pi_pred]),
-        )
-        cond = fpr != 0 # avoid infinite rejection (region of large uncertainty)
-        cond &= tpr != 1 # avoid linear interpolation (region of large uncertainty)
-        def mk_interp(tpr, fpr):
-            def interp(eff):
-                return np.interp(eff, tpr, fpr)
-            return interp
-        rocs.setdefault(clf_label, {})[energy] = mk_interp(tpr, fpr)
-        plt.plot(tpr[cond] * 100, 1 / fpr[cond], label=f"{clf_label}")
-        plt.yscale("log")
+    energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
+    def clf(events):
+        return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
+    e_pred = clf(e_eval[energy])
+    pi_pred = clf(pi_eval[energy])
+
+    for do_log, ax in [(False, axs[ix]), (True, axs_log[ix])]:
+        plt.sca(ax)
+        patch = plt.stairs(norm_by_sum(e_over_p_hist[("e-", energy)].values()), e_over_p_hist[("e-", energy)].axes[0].edges, label=rf"$e^-$ {PLOT_TITLE}")
+        patch = plt.stairs(norm_by_sum(e_over_p_hist[("pi-", energy)].values()), e_over_p_hist[("pi-", energy)].axes[0].edges, label=rf"$\pi^-$ {PLOT_TITLE}")
         plt.title(f"{energy}")
-        plt.legend(loc="lower left")
-        plt.xlabel("Electron efficiency, %")
-        plt.ylabel("Pion rejection factor")
+        plt.legend()
+        plt.xlabel("Classifier output")
+        plt.ylabel("Event yield")
+        if do_log:
+            plt.yscale("log")
+
+    plt.sca(axs_roc[ix])
+    fpr = np.cumsum(e_over_p_hist[("pi-", energy)].values()[::-1])
+    fpr /= fpr[-1]
+    tpr = np.cumsum(e_over_p_hist[("e-", energy)].values()[::-1])
+    tpr /= tpr[-1]
+    cond = fpr != 0 # avoid infinite rejection (region of large uncertainty)
+    cond &= tpr != 1 # avoid linear interpolation (region of large uncertainty)
+    def mk_interp(tpr, fpr):
+        def interp(eff):
+            return np.interp(eff, tpr, fpr)
+        return interp
+    rocs.setdefault(PLOT_TITLE, {})[energy] = mk_interp(tpr, fpr)
+    plt.plot(tpr[cond] * 100, 1 / fpr[cond], label=f"{PLOT_TITLE}")
+    plt.yscale("log")
+    plt.title(f"{energy}")
+    plt.legend(loc="lower left")
+    plt.xlabel("Electron efficiency, %")
+    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")
@@ -311,8 +347,9 @@ for clf_label, roc in rocs.items():
         marker=".",
         label=f"{clf_label}",
     )
+xmax = np.max(energy_values) * 1.05
+plt.xlim(0., xmax)
 plt.yscale("log")
-plt.title(INPUT_PATH_FORMAT)
 plt.legend()
 plt.xlabel("Energy, GeV")
 plt.ylabel("Pion rejection at 95%")

benchmarks/backwards_ecal/config.yml

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

benchmarks/backwards_ecal/requirements.txt

-awkward >= 2.4.0
-scikit-learn
-uproot >= 5.2.0
-vector
+awkward ~= 2.7.2
+dask ~= 2024.12.1
+dask_awkward ~= 2024.12.2
+dask_histogram ~= 2024.12.1
+distributed ~= 2024.12.1
+pyhepmc
+uproot ~= 5.5.1
+vector ~= 1.5.2

benchmarks/calo_pid/Snakefile

@@ -15,7 +15,7 @@ rule calo_pid_sim:
         ENERGY_MAX="[0-9]+[kMG]eV",
         THETA_MIN="[0-9]+",
         THETA_MAX="[0-9]+",
-        INDEX="\d{4}",
+        INDEX=r"\d{4}",
     params:
         N_EVENTS=1000,
         SEED=lambda wildcards: "1" + wildcards.INDEX,
@@ -23,6 +23,8 @@ rule calo_pid_sim:
         DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
         ENERGY_MIN=lambda wildcards: format_energy_for_dd4hep(wildcards.ENERGY_MIN),
         ENERGY_MAX=lambda wildcards: format_energy_for_dd4hep(wildcards.ENERGY_MAX),
+        THETA_MIN=lambda wildcards: wildcards.THETA_MIN,
+        THETA_MAX=lambda wildcards: wildcards.THETA_MAX,
         DD4HEP_HASH=get_spack_package_hash("dd4hep"),
         NPSIM_HASH=get_spack_package_hash("npsim"),
     cache: True
@@ -55,31 +57,34 @@ rule calo_pid_recon:
     log:
         "sim_output/calo_pid/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root.log",
     wildcard_constraints:
-        INDEX="\d{4}",
+        INDEX=r"\d{4}",
     shell: """
 set -m # monitor mode to prevent lingering processes
 exec env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \
   eicrecon {input} -Ppodio:output_file={output} \
-  -Ppodio:output_collections=MCParticles,EcalEndcapNRecHits,EcalEndcapNClusters,EcalEndcapNParticleIDInput_features,EcalEndcapNParticleIDTarget
+  -Ppodio:output_collections=MCParticles,EcalEndcapNRecHits,EcalEndcapNClusters,EcalEndcapNParticleIDInput_features,EcalEndcapNParticleIDTarget,EcalEndcapNParticleIDOutput_probability_tensor
 """
 
 
-rule calo_pid:
+rule calo_pid_input_list:
     input:
         electrons=expand(
-            "sim_output/calo_pid/{{DETECTOR_CONFIG}}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root",
-            PARTICLE=["e-"],
-            ENERGY=["100MeVto20GeV"],
-            PHASE_SPACE=["130to177deg"],
-            INDEX=range(100),
-        ),
-        pions=expand(
-            "sim_output/calo_pid/{{DETECTOR_CONFIG}}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root",
-            PARTICLE=["pi-"],
+            "sim_output/calo_pid/{{DETECTOR_CONFIG}}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root",
             ENERGY=["100MeVto20GeV"],
             PHASE_SPACE=["130to177deg"],
             INDEX=range(100),
         ),
+    output:
+        "listing/calo_pid/{DETECTOR_CONFIG}/{PARTICLE}.lst",
+    run:
+        with open(output[0], "wt") as fp:
+            fp.write("\n".join(input))
+
+
+rule calo_pid:
+    input:
+        electrons="listing/calo_pid/{DETECTOR_CONFIG}/e-.lst",
+        pions="listing/calo_pid/{DETECTOR_CONFIG}/pi-.lst",
         matplotlibrc=".matplotlibrc",
         script="benchmarks/calo_pid/calo_pid.py",
     output:

benchmarks/calo_pid/calo_pid.org

@@ -32,8 +32,8 @@ vector.register_awkward()
 #+begin_src jupyter-python :results silent
 DETECTOR_CONFIG=os.environ.get("DETECTOR_CONFIG")
 PLOT_TITLE=os.environ.get("PLOT_TITLE")
-INPUT_PIONS=os.environ.get("INPUT_PIONS", "").split(" ")
-INPUT_ELECTRONS=os.environ.get("INPUT_ELECTRONS", "").split(" ")
+INPUT_PIONS=os.environ.get("INPUT_PIONS")
+INPUT_ELECTRONS=os.environ.get("INPUT_ELECTRONS")
 
 output_dir=Path(os.environ.get("OUTPUT_DIR", "./"))
 output_dir.mkdir(parents=True, exist_ok=True)
@@ -75,8 +75,13 @@ def filter_pointing(events):
     cond = (part_momentum.eta[:,0] > -3.5) & (part_momentum.eta[:,0] < -2.)
     return events[cond]
 
-e = filter_pointing(uproot.concatenate({filename: "events" for filename in INPUT_ELECTRONS}, filter_name=["MCParticles.*", "*EcalEndcapN*"]))
-pi = filter_pointing(uproot.concatenate({filename: "events" for filename in INPUT_PIONS}, filter_name=["MCParticles.*", "*EcalEndcapN*"]))
+def readlist(path):
+    with open(path, "rt") as fp:
+        paths = [line.rstrip() for line in fp.readlines()]
+    return paths
+
+e = filter_pointing(uproot.concatenate({filename: "events" for filename in readlist(INPUT_ELECTRONS)}, filter_name=["MCParticles.*", "*EcalEndcapN*"]))
+pi = filter_pointing(uproot.concatenate({filename: "events" for filename in readlist(INPUT_PIONS)}, filter_name=["MCParticles.*", "*EcalEndcapN*"]))
 
 e_train = e[:len(pi)//2]
 pi_train = pi[:len(pi)//2]
@@ -358,6 +363,7 @@ if "_EcalEndcapNParticleIDOutput_probability_tensor_floatData" in pi_train.field
     eval_proba = ak.concatenate([pi_eval_proba, e_eval_proba])
 
     plt.hist(clf.predict_proba(eval_x.to_numpy())[:,1] - eval_proba[:,1].to_numpy())
+    plt.savefig(output_dir / f"proba_diff.pdf", bbox_inches="tight")
     plt.show()
 else:
     print("EcalEndcapNParticleIDOutput not present")

benchmarks/ecal_gaps/Snakefile

@@ -14,7 +14,7 @@ rule ecal_gaps_sim:
         PARTICLE="e-",
         ENERGY="(500MeV|5GeV|20GeV)",
         PHASE_SPACE="(3to50|45to135|130to177)deg",
-        INDEX="\d{4}",
+        INDEX=r"\d{4}",
     params:
         N_EVENTS=1000,
         SEED=lambda wildcards: "1" + wildcards.INDEX,
@@ -47,7 +47,7 @@ rule ecal_gaps_recon:
     log:
         "sim_output/ecal_gaps/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root.log",
     wildcard_constraints:
-        INDEX="\d{4}",
+        INDEX=r"\d{4}",
     params:
         DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
         EICRECON_HASH=get_spack_package_hash("eicrecon"),
@@ -74,6 +74,9 @@ rule ecal_gaps:
         ),
     output:
         directory("results/{DETECTOR_CONFIG}/ecal_gaps"),
+    log:
+        scheduler=".logs/results/{DETECTOR_CONFIG}/ecal_gaps/scheduler.log",
+        worker=".logs/results/{DETECTOR_CONFIG}/ecal_gaps/worker.log",
     threads: workflow.cores
     shell:
         """
@@ -85,10 +88,10 @@ cleanup() {{
 trap cleanup EXIT
 
 PORT=$RANDOM
-dask scheduler --port $PORT &
+dask scheduler --port $PORT 2>{log.scheduler} &
 export DASK_SCHEDULER=localhost:$PORT
 SCHEDULER_PID=$!
-dask worker tcp://$DASK_SCHEDULER --nworkers {threads} --nthreads 1 &
+dask worker tcp://$DASK_SCHEDULER --nworkers {threads} --nthreads 1 2>{log.worker} &
 WORKER_PID=$!
 env \
 MATPLOTLIBRC={input.matplotlibrc} \

benchmarks/ecal_gaps/requirements.txt

-awkward >= 2.4.0
-dask >= 2023
-dask_awkward
-dask_histogram
-distributed >= 2023
+awkward ~= 2.7.2
+dask ~= 2024.12.1
+dask_awkward ~= 2024.12.2
+dask_histogram ~= 2024.12.1
+distributed ~= 2024.12.1
 pyhepmc
-uproot ~= 5.2.0
+uproot ~= 5.5.1

benchmarks/femc_electron/Snakefile

@@ -14,10 +14,11 @@ rule femc_electron_generate:
         N_EVENTS=get_n_events,
         th_max=28,
         th_min=2.0,
+        P=lambda wildcards: wildcards.P,
     shell:
         """
 mkdir -p sim_output/femc_electron
-root -l -b -q '{input.script}({params.N_EVENTS},"{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., {params.P})'
 """
 
 rule femc_electron_simulate:

benchmarks/femc_photon/Snakefile

@@ -14,10 +14,11 @@ rule femc_photon_generate:
         N_EVENTS=get_n_events,
         th_max=28,
         th_min=2.0,
+        P=lambda wildcards: wildcards.P,
     shell:
         """
 mkdir -p sim_output/femc_photon
-root -l -b -q '{input.script}({params.N_EVENTS},"{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., {params.P})'
 """
 
 rule femc_photon_simulate:

benchmarks/femc_pi0/Snakefile

@@ -14,10 +14,11 @@ rule femc_pi0_generate:
         N_EVENTS=get_n_events,
         th_max=28,
         th_min=2.0,
+        P=lambda wildcards: wildcards.P,
     shell:
         """
 mkdir -p sim_output/femc_pi0
-root -l -b -q '{input.script}({params.N_EVENTS},"{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., {params.P})'
 """
 
 rule femc_pi0_simulate:

benchmarks/insert_muon/Snakefile

@@ -7,9 +7,10 @@ rule insert_muon_generate:
         NEVENTS_GEN=5000,
         th_max=7.0,
         th_min=1.7,
+        P=lambda wildcards: wildcards.P,
     shell:
         """
-root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "mu-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "mu-", {params.th_min}, {params.th_max}, 0., 360., {params.P})'
 """
 
 rule insert_muon_simulate:
@@ -20,6 +21,7 @@ rule insert_muon_simulate:
     output:
         SIM_FILE="sim_output/insert_muon/{DETECTOR_CONFIG}_sim_mu-_{P}GeV_{INDEX}.edm4hep.root",
     params:
+        INDEX=lambda wildcards: wildcards.INDEX,
         PHYSICS_LIST="FTFP_BERT",
         DETECTOR_PATH=os.environ["DETECTOR_PATH"],
         DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
@@ -32,7 +34,7 @@ NEVENTS_SIM=1000
 # Running simulation
 npsim \
    --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \
-   --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \
+   --skipNEvents $(( $NEVENTS_SIM * {params.INDEX} )) \
    --numberOfEvents $NEVENTS_SIM \
    --physicsList {params.PHYSICS_LIST} \
    --inputFiles {input.GEN_FILE} \

benchmarks/insert_neutron/Snakefile

@@ -7,10 +7,11 @@ rule insert_neutron_generate:
         NEVENTS_GEN=1000,
         th_max=5.7,
         th_min=2.0,
+        P=lambda wildcards: wildcards.P,
     shell:
         """
 mkdir -p sim_output/insert_neutron
-root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "neutron", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+root -l -b -q '{input.script}({params.NEVENTS_GEN},"{output.GEN_FILE}", "neutron", {params.th_min}, {params.th_max}, 0., 360., {params.P})'
 """
 
 rule insert_neutron_simulate:
@@ -21,6 +22,7 @@ rule insert_neutron_simulate:
     output:
         SIM_FILE="sim_output/insert_neutron/{DETECTOR_CONFIG}_sim_neutron_{P}GeV_{INDEX}.edm4hep.root",
     params:
+        INDEX=lambda wildcards: wildcards.INDEX,
         PHYSICS_LIST="FTFP_BERT",
         DETECTOR_PATH=os.environ["DETECTOR_PATH"],
         DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
@@ -33,7 +35,7 @@ NEVENTS_SIM=200
 # Running simulation
 npsim \
    --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \
-   --skipNEvents $(( $NEVENTS_SIM * {wildcards.INDEX} )) \
+   --skipNEvents $(( $NEVENTS_SIM * {params.INDEX} )) \
    --numberOfEvents $NEVENTS_SIM \
    --physicsList {params.PHYSICS_LIST} \
    --inputFiles {input.GEN_FILE} \

benchmarks/insert_neutron/analysis/neutron_plots.py

@@ -79,11 +79,12 @@ slc=abs(bc)<3
 fnc=gauss
 sigma=np.sqrt(y[slc])+(y[slc]==0)
 p0=(100, 0, 5)
-coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
-xx=np.linspace(-5,5,100)
-plt.plot(xx,fnc(xx,*coeff))
-# except:
-#     pass
+try:
+    coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
+    xx=np.linspace(-5,5,100)
+    plt.plot(xx,fnc(xx,*coeff))
+except RuntimeError:
+    print("fit failed")
 plt.xlabel("$\\theta_{rec}-\\theta_{truth}$ [mrad]")
 plt.ylabel("events")
 plt.title(f"$p={p}$ GeV, ${eta_min}<\\eta<{eta_max}$")

benchmarks/insert_tau/Snakefile

@@ -13,9 +13,10 @@ rule insert_tau_generate:
         N_EVENTS=get_n_events,
         th_max=7.0,
         th_min=1.7,
+        P=lambda wildcards: wildcards.P,
     shell:
         """
-root -l -b -q '{input.script}({params.N_EVENTS},"{output.GEN_FILE}", "tau-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+root -l -b -q '{input.script}({params.N_EVENTS},"{output.GEN_FILE}", "tau-", {params.th_min}, {params.th_max}, 0., 360., {params.P})'
 """
 
 rule insert_tau_simulate:
@@ -27,6 +28,7 @@ rule insert_tau_simulate:
         SIM_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_sim_tau-_{P}GeV_{INDEX}.edm4hep.root",
     params:
         N_EVENTS=get_n_events,
+        INDEX=lambda wildcards: wildcards.INDEX,
         PHYSICS_LIST="FTFP_BERT",
         DETECTOR_PATH=os.environ["DETECTOR_PATH"],
         DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
@@ -39,7 +41,7 @@ rule insert_tau_simulate:
 npsim \
    --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \
    --numberOfEvents {params.N_EVENTS} \
-   --skipNEvents $(( {params.N_EVENTS} * {wildcards.INDEX} )) \
+   --skipNEvents $(( {params.N_EVENTS} * {params.INDEX} )) \
    --physicsList {params.PHYSICS_LIST} \
    --inputFiles {input.GEN_FILE} \
    --outputFile {output.SIM_FILE}