diff --git a/km3buu/output.py b/km3buu/output.py
index 1c0802ec679b3fad8f1680ddb5f0acc2c8501647..31766f0b30178cd19141ff6ca4e866d3393a516b 100644
--- a/km3buu/output.py
+++ b/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]
diff --git a/km3buu/tests/test_output.py b/km3buu/tests/test_output.py
index cb8de46664731163fc919035027aa4c3bb0cd3f4..e6ff083e86bf2be85e29326e899867684426625f 100644
--- a/km3buu/tests/test_output.py
+++ b/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(