ML model for e/pi separation in the BIC

.gitlab-ci.yml

@@ -128,6 +128,7 @@ include:
   - local: 'benchmarks/backwards_ecal/config.yml'
   - local: 'benchmarks/beamline/config.yml'
   - local: 'benchmarks/calo_pid/config.yml'
+  - local: 'benchmarks/bic_pid/config.yml'
   - local: 'benchmarks/campaign/config.yml'
   - local: 'benchmarks/ecal_gaps/config.yml'
   - local: 'benchmarks/far_forward_dvcs/config.yml'

Snakefile

@@ -2,6 +2,7 @@ configfile: "snakemake.yml"
 
 import functools
 import os
+import subprocess
 from snakemake.logging import logger
 
 
@@ -48,6 +49,7 @@ include: "benchmarks/backwards_ecal/Snakefile"
 include: "benchmarks/barrel_ecal/Snakefile"
 include: "benchmarks/beamline/Snakefile"
 include: "benchmarks/calo_pid/Snakefile"
+include: "benchmarks/bic_pid/Snakefile"
 include: "benchmarks/campaign/Snakefile"
 include: "benchmarks/ecal_gaps/Snakefile"
 include: "benchmarks/far_forward_dvcs/Snakefile"

benchmarks/bic_pid/Snakefile

@@ -0,0 +1,129 @@
+from glob import glob
+import os
+
+def format_energy_for_dd4hep(s):
+    return s.rstrip("kMGeV") + "*" + s.lstrip("0123456789")
+
+
+def theta_min_from_phase_space(s):
+    return s.replace("deg", "").split("to")[0]
+
+
+def theta_max_from_phase_space(s):
+    return s.replace("deg", "").split("to")[1]
+
+
+rule bic_pid_sim:
+    input:
+        warmup="warmup.edm4hep.root",
+        geometry_lib=find_epic_libraries(),
+    output:
+        "sim_output/bic_pid/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root",
+    log:
+        "sim_output/bic_pid/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root.log",
+    wildcard_constraints:
+        PARTICLE="(e-|pi-)",
+        ENERGY="[0-9]+[kMG]eV",
+        PHASE_SPACE="[0-9]+to[0-9]+deg",
+        INDEX=r"\d{4}",
+    params:
+        N_EVENTS=1000,
+        SEED=lambda wildcards: "1" + wildcards.INDEX,
+        DETECTOR_PATH=os.environ["DETECTOR_PATH"],
+        DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
+        ENERGY=lambda wildcards: format_energy_for_dd4hep(wildcards.ENERGY),
+        THETA_MIN=lambda wildcards: theta_min_from_phase_space(wildcards.PHASE_SPACE),
+        THETA_MAX=lambda wildcards: theta_max_from_phase_space(wildcards.PHASE_SPACE),
+        DD4HEP_HASH=get_spack_package_hash("dd4hep"),
+        NPSIM_HASH=get_spack_package_hash("npsim"),
+    cache: True
+    shell:
+        r"""
+set -m 
+npsim \
+  --runType batch \
+  --enableGun \
+  --gun.momentumMin "{params.ENERGY}" \
+  --gun.momentumMax "{params.ENERGY}" \
+  --gun.thetaMin "{params.THETA_MIN}*deg" \
+  --gun.thetaMax "{params.THETA_MAX}*deg" \
+  --gun.particle {wildcards.PARTICLE} \
+  --gun.distribution eta \
+  --random.seed {params.SEED} \
+  --filter.tracker edep0 \
+  -v WARNING \
+  --numberOfEvents {params.N_EVENTS} \
+  --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \
+  --outputFile {output} > {log} 2>&1
+"""
+
+
+rule bic_pid_recon:
+    input:
+        sim="sim_output/bic_pid/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root",
+        warmup="warmup.edm4hep.root",
+    output:
+        "sim_output/bic_pid/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.edm4eic.root",
+    log:
+        "sim_output/bic_pid/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.edm4eic.root.log",
+    wildcard_constraints:
+        PARTICLE="(e-|pi-)",
+        ENERGY="[0-9]+[kMG]eV",
+        PHASE_SPACE="[0-9]+to[0-9]+deg",
+        INDEX=r"\d{4}",
+    params:
+        DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG,
+        EICRECON_HASH=get_spack_package_hash("eicrecon"),
+    cache: True
+    shell:
+        """
+DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \
+exec eicrecon {input.sim} -Ppodio:output_file={output} \
+  -Ppodio:output_collections=MCParticles,EcalBarrelScFiRecHits,EcalBarrelImagingRecHits
+"""
+
+
+rule bic_pid_input_list:
+    output:
+        "listing/bic_pid/{DETECTOR_CONFIG}/{PARTICLE}.lst",
+    params:
+        energy="1GeV",
+        phase_space="45to135deg",
+    run:
+        pattern = (
+            f"sim_output/bic_pid/{wildcards.DETECTOR_CONFIG}/{wildcards.PARTICLE}/"
+            f"{params.energy}/{params.phase_space}/"
+            f"{wildcards.PARTICLE}_{params.energy}_{params.phase_space}.*.eicrecon.edm4eic.root"
+        )
+        files = sorted(glob(pattern))
+
+        if len(files) != 100:
+            raise ValueError(
+                f"Expected 100 files for {wildcards.PARTICLE}, found {len(files)}.\n"
+                f"Pattern used: {pattern}"
+            )
+
+        os.makedirs(os.path.dirname(output[0]), exist_ok=True)
+        with open(output[0], "wt") as fp:
+            fp.write("\n".join(files))
+
+
+rule bic_pid:
+    input:
+        electrons="listing/bic_pid/{DETECTOR_CONFIG}/e-.lst",
+        pions="listing/bic_pid/{DETECTOR_CONFIG}/pi-.lst",
+        matplotlibrc=".matplotlibrc",
+        script="benchmarks/bic_pid/bic_pid.py",
+    output:
+        directory("results/{DETECTOR_CONFIG}/bic_pid")
+    shell:
+        """
+env \
+MATPLOTLIBRC={input.matplotlibrc} \
+DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \
+PLOT_TITLE={wildcards.DETECTOR_CONFIG} \
+INPUT_ELECTRONS="{input.electrons}" \
+INPUT_PIONS="{input.pions}" \
+OUTPUT_DIR={output} \
+python {input.script}
+"""