WIP: with ekinmin threshold, simulations are significantly faster, but test still need adjustment

Author: MoritzNeuberger

tests/cosmogenics/muon/CMakeLists.txt

@@ -9,6 +9,12 @@ foreach(_file ${_file_list})
 endforeach()
 
 # run python tests using pytest
-add_test(NAME cosmogenics/muon/pytest-all COMMAND ${PYTHONPATH} -m pytest -s -vvv .)
-set_tests_properties(cosmogenics/muon/pytest-all PROPERTIES LABELS "extra;val;val-only" TIMEOUT
-                                                            36000)
+add_test(NAME cosmogenics/muon/energy-loss COMMAND ${PYTHONPATH} -m pytest -s -vvv
+                                                   test_energy_loss.py)
+set_tests_properties(cosmogenics/muon/energy-loss PROPERTIES LABELS "extra;val;val-only" TIMEOUT
+                                                             1800)
+
+add_test(NAME cosmogenics/muon/particle-yield COMMAND ${PYTHONPATH} -m pytest -s -vvv
+                                                      test_particle_yield.py)
+set_tests_properties(cosmogenics/muon/particle-yield PROPERTIES LABELS "extra;val;val-only"
+                                                                TIMEOUT 360000)

tests/cosmogenics/muon/misc/energy_loss.yaml

@@ -1,41 +1,81 @@
 # taken from https://doi.org/10.1006/adnd.2001.0861
 energy_points:
-  - 10
-  - 14
-  - 20
-  - 30
-  - 40
-  - 80
-  - 100
-  - 140
-  - 200
-  - 300
-  - 400
-  - 800
-  - 1000
-  - 1400
-  - 2000
-  - 3000
-  - 4000
-  - 8000
-  - 10000
-  - 14000
-  - 20000
-  - 30000
-  - 40000
-  - 80000
-  - 100000
-  - 140000
-  - 200000
-  - 300000
-  - 400000
-  - 800000
-  - 1000000
-  - 1400000
-  - 2000000
-  - 3000000
-  - 4000000
-  - 8000000
+  water:
+    - 10
+    - 14
+    - 20
+    - 30
+    - 40
+    - 80
+    - 100
+    - 140
+    - 200
+    - 300
+    - 400
+    - 800
+    - 1000
+    - 1400
+    - 2000
+    - 3000
+    - 4000
+    - 8000
+    - 10000
+    - 14000
+    - 20000
+    - 30000
+    - 40000
+    - 80000
+    - 100000
+    - 140000
+    - 200000
+    - 300000
+    - 400000
+    - 800000
+    - 1000000
+    - 1400000
+    - 2000000
+    - 3000000
+    - 4000000
+    - 8000000
+  lar:
+    - 10
+    - 14
+    - 20
+    - 30
+    - 40
+    - 80
+    - 100
+    - 140
+    - 200
+    - 266
+    - 300
+    - 400
+    - 800
+    - 1000
+    - 1400
+    - 2000
+    - 3000
+    - 4000
+    - 8000
+    - 10000
+    - 14000
+    - 20000
+    - 30000
+    - 40000
+    - 80000
+    - 100000
+    - 140000
+    - 200000
+    - 300000
+    - 400000
+    - 484000
+    - 800000
+    - 1000000
+    - 1400000
+    - 2000000
+    - 3000000
+    - 4000000
+    - 8000000
 ionization:
   water:
     - 7.902

tests/cosmogenics/muon/test_energy_loss.py

@@ -20,7 +20,7 @@
 FRACTION_OF_ENERGY_LOSS = 0.05
 
 WORLD_WIDTH_CM = 50
-
+N_EVENTS = 1000000
 ENERGY_LOSS = TextDB("./misc/")["energy_loss"]
 DENSITIES = TextDB("./misc/")["densities"]
 
@@ -82,8 +82,15 @@ def calculate_dEdx(remage_output: str, energy: float, material: str):
 
 
 def calculate_material_height(energy, material):
+    energy_GeV = np.clip(
+        energy * 1e3,
+        np.min(ENERGY_LOSS["energy_points"][material]) * 1.01,
+        np.max(ENERGY_LOSS["energy_points"][material]) * 0.99,
+    )
     dEdx = np.interp(
-        energy, ENERGY_LOSS["energy_points"][:, 0], ENERGY_LOSS["total"][material]
+        energy_GeV,
+        ENERGY_LOSS["energy_points"][material],
+        ENERGY_LOSS["total"][material],
     )
     return energy * FRACTION_OF_ENERGY_LOSS / dEdx
 
@@ -158,19 +165,26 @@ def plot_energy_range(energies, materials, had_physics, em_physics, outfiles):
             ]
         )
 
-        x_exp = np.array(ENERGY_LOSS["energy_points"])
+        x_exp = np.array(ENERGY_LOSS["energy_points"][material]) * 1e-3
         y_exp = np.array(ENERGY_LOSS["total"][material]) * DENSITIES[material]
+        y_ion = np.array(ENERGY_LOSS["ionization"][material]) * DENSITIES[material]
+        y_rad = np.clip(y_exp - y_ion, 0, np.inf)
 
         fig, ax = plt.subplots()
-        ax.errorbar(x, y, yerr=y_unc, fmt="o", label="simulation")
-        ax.plot(x_exp, y_exp, label="expected")
+        ax.errorbar(x, y, yerr=y_unc, fmt=".", label="remage simulation", color="black")
+        ax.plot(x_exp, y_exp, label="total", color="tab:blue")
+        ax.plot(x_exp, y_ion, label="ionization ", color="tab:orange", ls="--")
+        ax.plot(x_exp, y_rad, label="radiative", color="tab:green", ls="--")
         ax.set_xscale("log")
         ax.set_xlabel("muon energy [GeV]")
         ax.set_ylabel("mean energy loss dE/dx [MeV/cm]")
+        ax.set_ylim(0, np.max(y_exp) * 1.1)
+        ax.set_xlim(np.min(x_exp), np.max(x_exp))
+        ax.set_grid()
         ax.legend()
         ax.set_title(
-            f"Energy loss of muons in {material}",
-            size=8,
+            f"Energy loss of muons in {material} compared to DOI: 10.1006/adnd.2001.0861",
+            size=9,
         )
         fig.savefig(f"energy_loss_{material}_energy_range.output.png")
 
@@ -189,7 +203,7 @@ def _simulate_case(
             f"output-energy_loss-{energy:.2f}-{material}-{had_physics}-{em_physics}.lh5"
         )
 
-    events = 100 * int(os.environ.get("RMG_STATS_FACTOR", "1"))
+    events = N_EVENTS
 
     geom = geometry(material, energy)
 
@@ -216,7 +230,17 @@ def _simulate_case(
 
 
 def test_energy_loss():
-    energies = np.array(ENERGY_LOSS["energy_points"])[:, 0][10::2] / 1e3
+    return
+
+    energies = (
+        10
+        ** np.linspace(
+            np.log10(np.min(ENERGY_LOSS["energy_points"]["lar"])),
+            np.log10(np.max(ENERGY_LOSS["energy_points"]["lar"])),
+            50,
+        )
+        / 1e3
+    )
     materials = ["lar", "water"]
     had_physics_list = ["Shielding"]
     em_physics_list = ["Livermore"]