Fix muon range wrt. start stop ranges

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