Refine flux determination method

km3buu/output.py

@@ -583,10 +583,13 @@ class GiBUUOutput:
         def fluxfunc(x, a, b):
             return a * x**b
 
-        mask = ~np.isclose(self.flux_data["events"], 0)
+        lower_limit = np.exp(np.log(np.max(self.flux_data["flux"])) * 0.2)
+        upper_limit = np.exp(np.log(np.max(self.flux_data["flux"])) * 0.8)
+        mask = (self.flux_data["flux"] > lower_limit) & (self.flux_data["flux"]
+                                                         < upper_limit)
         popt, pcov = curve_fit(fluxfunc,
                                self.flux_data["energy"][mask],
-                               self.flux_data["events"][mask],
+                               self.flux_data["flux"][mask],
                                p0=[1, -1])
 
         self._flux_index = popt[1]

km3buu/tests/test_output.py

@@ -73,6 +73,9 @@ class TestGiBUUOutput(unittest.TestCase):
         assert self.output.Z == 8
         assert self.output.A == 16
 
+    def test_flux_index(self):
+        assert np.isclose(self.output.flux_index, -0.904, rtol=1e-3)
+
     def test_w2weights(self):
         w2 = self.output.w2weights(123.0, 2.6e28, 4 * np.pi)
         np.testing.assert_array_almost_equal(