backgrounds: use updated beam gas files

.gitlab-ci.yml

@@ -128,7 +128,7 @@ get_data:
       - runner_system_failure
 
 include: 
-#  - local: 'benchmarks/backgrounds/config.yml'
+  - local: 'benchmarks/backgrounds/config.yml'
   - local: 'benchmarks/ecal_gaps/config.yml'
   - local: 'benchmarks/tracking_detectors/config.yml'
   - local: 'benchmarks/barrel_ecal/config.yml'

benchmarks/backgrounds/Snakefile

@@ -2,44 +2,13 @@ import os
 import shutil
 
 
-rule backgrounds_get_beam_gas_electron:
-    input:
-        provider.remote(remote_path("EPIC/EVGEN/BACKGROUNDS/BEAMGAS/electron/beam_gas_ep_10GeV_foam_emin10keV_10Mevt_vtx.hepmc")),
-    output:
-        "input/backgrounds/beam_gas_electron.hepmc",
-    run:
-        shutil.move(input[0], output[0])
-
-
-rule backgrounds_get_beam_gas_proton:
-    input:
-        provider.remote(remote_path("EPIC/EVGEN/BACKGROUNDS/BEAMGAS/proton/ProtonBeamGasEvents/100GeV/100GeV_1.hepmc")),
-    output:
-        "input/backgrounds/beam_gas_proton.hepmc",
-    run:
-        shutil.move(input[0], output[0])
-
-
-rule backgrounds_get_DIS:
-    input:
-        provider.remote(remote_path("EPIC/EVGEN/DIS/NC/{BEAM}/minQ2={MINQ2}/pythia8NCDIS_{BEAM}_minQ2={MINQ2}_{SUFFIX}.hepmc")),
-    wildcard_constraints:
-        BEAM="\d+x\d+",
-        MINQ2="\d+",
-    output:
-        "input/backgrounds/pythia8NCDIS_{BEAM}_minQ2={MINQ2}_{SUFFIX}.hepmc",
-    run:
-        shutil.move(input[0], output[0])
-
-
 rule backgrounds_sim:
     input:
-        hepmc="input/backgrounds/{NAME}.hepmc",
         warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
     output:
-        "sim_output/{DETECTOR_CONFIG}/backgrounds/{NAME}.edm4hep.root",
+        "sim_output/{DETECTOR_CONFIG}/backgrounds/{PATH}.edm4hep.root",
     log:
-        "sim_output/{DETECTOR_CONFIG}/backgrounds/{NAME}.edm4hep.root.log",
+        "sim_output/{DETECTOR_CONFIG}/backgrounds/{PATH}.edm4hep.root.log",
     params:
         N_EVENTS=100
     shell:
@@ -51,7 +20,7 @@ ddsim \
   -v WARNING \
   --numberOfEvents {params.N_EVENTS} \
   --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
-  --inputFiles {input.hepmc} \
+  --inputFiles root://dtn-eic.jlab.org//work/eic2/{wildcards.PATH}.hepmc3.tree.root \
   --outputFile {output}
 """
 
@@ -62,16 +31,21 @@ rule backgrounds_ecal_backwards:
     input:
         matplotlibrc=".matplotlibrc",
         script="benchmarks/backgrounds/ecal_backwards.py",
-        electron_beam_gas_gen="input/backgrounds/beam_gas_electron.hepmc",
-        electron_beam_gas_sim="sim_output/" + DETECTOR_CONFIG + "/backgrounds/beam_gas_electron.edm4hep.root",
-        physics_signal_sim="sim_output/" + DETECTOR_CONFIG + "/backgrounds/pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_1.edm4hep.root",
-        proton_beam_gas_gen="input/backgrounds/beam_gas_proton.hepmc",
-        proton_beam_gas_sim="sim_output/" + DETECTOR_CONFIG + "/backgrounds/beam_gas_proton.edm4hep.root",
+        electron_beam_gas_sim="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",
+        physics_signal_sim="sim_output/" + DETECTOR_CONFIG + "/backgrounds/EPIC/EVGEN/DIS/NC/10x100/minQ2=1/pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_1.edm4hep.root",
+        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")
     threads: workflow.cores
     shell:
         """
+set -m # monitor mode to prevent lingering shells
+cleanup() {{
+  echo Cleaning up
+  kill $WORKER_PID $SCHEDULER_PID
+}}
+trap cleanup EXIT
+
 PORT=$RANDOM
 dask scheduler --port $PORT &
 export DASK_SCHEDULER=localhost:$PORT
@@ -80,12 +54,11 @@ dask worker tcp://$DASK_SCHEDULER --nworkers {threads} --nthreads 1 &
 WORKER_PID=$!
 env \
 MATPLOTLIBRC={input.matplotlibrc} \
-ELECTRON_BEAM_GAS_GEN=$(realpath {input.electron_beam_gas_gen}) \
+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_SIM=$(realpath {input.electron_beam_gas_sim}) \
 PHYSICS_PROCESS_SIM=$(realpath {input.physics_signal_sim}) \
-PROTON_BEAM_GAS_GEN=$(realpath {input.proton_beam_gas_gen}) \
+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_SIM=$(realpath {input.proton_beam_gas_sim}) \
 OUTPUT_DIR={output} \
 python {input.script}
-kill $WORKER_PID $SCHEDULER_PID
 """

benchmarks/backgrounds/config.yml

@@ -4,7 +4,11 @@ sim:backgrounds:
   script:
     - mkdir -p $LOCAL_DATA_PATH/input
     - ln -s $LOCAL_DATA_PATH/input input
-    - snakemake --cores 2 sim_output/$DETECTOR_CONFIG/beam_gas_{electron,proton}.edm4hep.root sim_output/$DETECTOR_CONFIG/pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_vtxfix_1.edm4hep.root
+    - |
+      snakemake --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
 
 bench:backgrounds_emcal_backwards:
   extends: .det_benchmark

benchmarks/backgrounds/ecal_backwards.org

@@ -15,7 +15,7 @@ import dask_awkward as dak
 import dask_histogram as dh
 import numpy as np
 import uproot
-import pyhepmc
+from pyHepMC3 import HepMC3
 #+end_src
 
 #+begin_src jupyter-python :results silent
@@ -66,47 +66,47 @@ output_dir.mkdir(parents=True, exist_ok=True)
 #+begin_src jupyter-python :results silent
 builder = ak.ArrayBuilder()
 n = 0
-with pyhepmc.open(PROTON_BEAM_GAS_GEN) as f:
-    it = iter(f)
-    for event in f:
-        builder.begin_list()
-        for vertex in event.vertices:
-            assert event.length_unit.name == "MM"
-            with builder.record("Vector4D"):
-                builder.field("x")
-                builder.real(vertex.position.x)
-                builder.field("y")
-                builder.real(vertex.position.y)
-                builder.field("z")
-                builder.real(vertex.position.z)
-                builder.field("t")
-                builder.real(vertex.position.t)
-        builder.end_list()
-        n += 1
-        if n > 10000: break
+f = HepMC3.ReaderPlugin(PROTON_BEAM_GAS_GEN, "libHepMC3rootIO.so", "newReaderRootTreefile")
+event = HepMC3.GenEvent()
+while f.read_event(event):
+    builder.begin_list()
+    assert event.length_unit().name == "MM"
+    for vertex in event.vertices():
+        with builder.record("Vector4D"):
+            builder.field("x")
+            builder.real(vertex.position().x())
+            builder.field("y")
+            builder.real(vertex.position().y())
+            builder.field("z")
+            builder.real(vertex.position().z())
+            builder.field("t")
+            builder.real(vertex.position().t())
+    builder.end_list()
+    n += 1
+    if n > 10000: break
 
 vertices_proton_beam_gas = builder.snapshot()
 
 builder = ak.ArrayBuilder()
 n = 0
-with pyhepmc.open(ELECTRON_BEAM_GAS_GEN) as f:
-    it = iter(f)
-    for event in f:
-        builder.begin_list()
-        for vertex in event.vertices:
-            assert event.length_unit.name == "MM"
-            with builder.record("Vector3D"):
-                builder.field("x")
-                builder.real(vertex.position.x)
-                builder.field("y")
-                builder.real(vertex.position.y)
-                builder.field("z")
-                builder.real(vertex.position.z)
-                builder.field("t")
-                builder.real(vertex.position.t)
-        builder.end_list()
-        n += 1
-        if n > 10000: break
+f = HepMC3.ReaderPlugin(ELECTRON_BEAM_GAS_GEN, "libHepMC3rootIO.so", "newReaderRootTreefile")
+event = HepMC3.GenEvent()
+while f.read_event(event):
+    builder.begin_list()
+    assert event.length_unit().name == "MM"
+    for vertex in event.vertices():
+        with builder.record("Vector4D"):
+            builder.field("x")
+            builder.real(vertex.position().x())
+            builder.field("y")
+            builder.real(vertex.position().y())
+            builder.field("z")
+            builder.real(vertex.position().z())
+            builder.field("t")
+            builder.real(vertex.position().t())
+    builder.end_list()
+    n += 1
+    if n > 10000: break
 
 vertices_electron_beam_gas = builder.snapshot()
 #+end_src
@@ -153,6 +153,7 @@ plt.xlabel("vertex[0].t, mm")
 plt.legend()
 plt.savefig(output_dir / "vertex_time_distribution.png", bbox_inches="tight")
 plt.show()
+plt.clf()
 
 for label, ds in datasets.items():
     if "vertices" not in ds: continue
@@ -164,6 +165,7 @@ plt.xlabel("vertex[0].z, mm")
 plt.legend()
 plt.savefig(output_dir / "vertex_z_distribution.png", bbox_inches="tight")
 plt.show()
+plt.clf()
 
 for label, ds in datasets.items():
     if "vertices" not in ds: continue
@@ -178,6 +180,7 @@ plt.ylabel("vertex[0].x, mm")
 plt.legend()
 plt.savefig(output_dir / "vertex_xz_distribution.png", bbox_inches="tight")
 plt.show()
+plt.clf()
 
 for ix, (label, ds) in enumerate(datasets.items()):
     if "vertices" not in ds: continue
@@ -192,6 +195,7 @@ for ix, (label, ds) in enumerate(datasets.items()):
     plt.title(label)
     plt.savefig(output_dir / f"vertex_yz_distribution_{ix}.png", bbox_inches="tight")
     plt.show()
+    plt.clf()
 #+end_src
 
 ** Simulation results
@@ -210,7 +214,8 @@ for collection_name in ["EcalEndcapNHits", "EcalEndcapPHits"]:
         for ix, vals in enumerate([xs, ys]):
             centers = np.unique(vals)
             diffs = centers[1:] - centers[:-1]
-            bin_widths[ix] = np.min(diffs[diffs > 0]) if np.sum(diffs > 0) > 0 else 1.
+            EPSILON = 1e-5
+            bin_widths[ix] = np.min(diffs[diffs > EPSILON]) if np.sum(diffs > EPSILON) > 0 else 1.
             print(f"bin_widths[{ix}]", bin_widths[ix])
 
         bins = {
@@ -231,6 +236,7 @@ for collection_name in ["EcalEndcapNHits", "EcalEndcapPHits"]:
         plt.ylabel("hit y, mm", loc="top")
         plt.savefig(output_dir / f"{collection_name}_event_display_{dataset_ix}.png", bbox_inches="tight")
         plt.show()
+        plt.clf()
 #+end_src
 
 ** Discovering number of cells
@@ -285,6 +291,7 @@ for collection_name in ["EcalEndcapNHits", "EcalEndcapPHits"]:
     plt.ylabel(f"Number of towers per {dr} mm slice in $r$", loc="top")
     plt.savefig(output_dir / f"{collection_name}_num_towers.png", bbox_inches="tight")
     plt.show()
+    plt.clf()
 #+end_src
 
 ** Plotting the rates
@@ -324,6 +331,7 @@ for collection_name in ["EcalEndcapNHits", "EcalEndcapPHits"]:
         plt.yscale("log")
         plt.savefig(output_dir / f"{collection_name}_hit_rate_vs_r_edep_min_{edep_min:.3f}.png", bbox_inches="tight")
         plt.show()
+        plt.clf()
 #+end_src
 
 #+begin_src jupyter-python
@@ -354,6 +362,7 @@ for collection_name in ["EcalEndcapNHits", "EcalEndcapPHits"]:
         plt.yscale("log")
         plt.savefig(output_dir / f"{collection_name}_hit_rate_vs_z_totedep_min_{totedep_min:.1f}.png", bbox_inches="tight")
         plt.show()
+        plt.clf()
 #+end_src
 
 #+begin_src jupyter-python
@@ -444,4 +453,5 @@ for edep_min in [0]: # GeV
     plt.xticks(rotation=90, ha='right')
     fig_cmb.savefig(f"rates_edep_min_{edep_min}.png", bbox_inches="tight")
     plt.show()
+    plt.clf()
 #+end_src

benchmarks/ecal_gaps/Snakefile

@@ -66,6 +66,7 @@ rule ecal_gaps:
     threads: workflow.cores
     shell:
         """
+set -m # monitor mode to prevent lingering shells
 cleanup() {{
   echo Cleaning up
   kill $WORKER_PID $SCHEDULER_PID