From 6d808bbe219bdc243c63b3f0a9fa0835704ef99f Mon Sep 17 00:00:00 2001
From: Johannes Schumann <johannes.schumann@fau.de>
Date: Sun, 19 Dec 2021 02:57:20 +0100
Subject: [PATCH] Fix muon range wrt. start stop ranges

---
 km3buu/physics.py | 32 +++++++++++++++++---------------
 1 file changed, 17 insertions(+), 15 deletions(-)

diff --git a/km3buu/physics.py b/km3buu/physics.py
index 0aaaf53..52d931f 100644
--- a/km3buu/physics.py
+++ b/km3buu/physics.py
@@ -279,25 +279,27 @@ def muon_range_seawater(start_energy, stop_energy):
     muon_mass = Particle.from_string("mu").mass / 1e3
     if start_energy <= muon_mass:
         return 0
+    elif start_energy < stop_energy:
+        raise ValueError("Final energy must be smaller than initial energy.")
 
     etmp = start_energy
     dx = 0.0
 
-    params = np.array(
-        [(0., 2.3e-1 * DENSITY_SEA_WATER, 15.5e-4 * DENSITY_SEA_WATER),
-         (30., 2.67e-1 * DENSITY_SEA_WATER, 3.4e-4 * DENSITY_SEA_WATER),
-         (35.3e3, -6.5e-1 * DENSITY_SEA_WATER, 3.66e-4 * DENSITY_SEA_WATER)],
-        dtype=[("energy", "f8"), ("a", "f8"), ("b", "f8")])
-
-    mask_stages = (params["energy"] < stop_energy) & (np.append(
-        params["energy"][1:], np.inf) > start_energy)
-
-    energy_steps = np.append(params["energy"][mask_stages], stop_energy)
-    energy_steps[0] = start_energy
-
-    for i in np.where(mask_stages)[0]:
-        dx += muon_range(energy_steps[i], energy_steps[i + 1], params["a"][i],
-                         params["b"][i])
+    params = np.array([
+        (35.3e3, -6.5e-1 * DENSITY_SEA_WATER, 3.66e-4 * DENSITY_SEA_WATER),
+        (30., 2.67e-1 * DENSITY_SEA_WATER, 3.4e-4 * DENSITY_SEA_WATER),
+        (0., 2.3e-1 * DENSITY_SEA_WATER, 15.5e-4 * DENSITY_SEA_WATER),
+    ],
+                      dtype=[("energy", "f8"), ("a", "f8"), ("b", "f8")])
+
+    mask_stages = (params["energy"] < start_energy) & (np.append(
+        np.inf, params["energy"][:-1]) > stop_energy)
+    energy_steps = np.append(start_energy, params["energy"][mask_stages])
+    energy_steps[-1] = stop_energy
+
+    for i, stage in enumerate(np.where(mask_stages)[0]):
+        dx += muon_range(energy_steps[i], energy_steps[i + 1],
+                         params["a"][stage], params["b"][stage])
 
     return dx
 
-- 
GitLab