Flux Index Fit

CHANGELOG.rst

 Unreleased changes
 ------------------
+* Improved curve for flux index
 
 v1.1.2
 ----------------------------

km3buu/output.py

@@ -16,7 +16,7 @@ import re
 import numpy as np
 from io import StringIO
 from os import listdir, environ
-from os.path import isfile, join, abspath
+from os.path import isfile, join, abspath, exists
 from tempfile import TemporaryDirectory
 import awkward as ak
 import uproot
@@ -237,10 +237,9 @@ class GiBUUOutput:
         self._generated_events = -1
         self._flux_index = np.nan
 
-        try:
-            self._read_flux_file()
+        if self._read_flux_file():
             self._determine_flux_index()
-        except OSError:
+        else:
             self._read_single_energy()
 
     def _read_root_output(self):
@@ -285,6 +284,8 @@ class GiBUUOutput:
 
     def _read_flux_file(self):
         fpath = join(self._data_path, FLUXDESCR_FILENAME)
+        if not exists(fpath):
+            return False
         self.flux_data = np.loadtxt(fpath,
                                     dtype=FLUX_INFORMATION_DTYPE,
                                     usecols=(
@@ -297,6 +298,7 @@ class GiBUUOutput:
         self._energy_min = np.min(self.flux_data["energy"])
         self._energy_max = np.max(self.flux_data["energy"])
         self._generated_events = int(np.sum(self.flux_data["events"]))
+        return True
 
     def _event_xsec(self, root_tupledata):
         weights = np.array(root_tupledata.weight)
@@ -550,15 +552,19 @@ class GiBUUOutput:
         def fluxfunc(x, a, b):
             return a * x**b
 
-        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["flux"][mask],
-                               p0=[1, -1])
-        self._flux_index = popt[1]
+        energy_mask = self.flux_data["flux"] > 0
+        lower_limit = np.min(self.flux_data["energy"][energy_mask]) * 1.2
+        upper_limit = np.max(self.flux_data["energy"][energy_mask]) * 0.8
+        mask = (self.flux_data["energy"]
+                > lower_limit) & (self.flux_data["energy"] < upper_limit)
+        try:
+            popt, pcov = curve_fit(fluxfunc,
+                                   self.flux_data["energy"][mask],
+                                   self.flux_data["flux"][mask],
+                                   p0=[1, -1])
+            self._flux_index = popt[1]
+        except:
+            self._flux_index = np.nan
 
     @property
     def flux_index(self):

km3buu/tests/test_output.py

@@ -77,7 +77,7 @@ class TestGiBUUOutput(unittest.TestCase):
         assert self.output.A == 16
 
     def test_flux_index(self):
-        assert np.isclose(self.output.flux_index, -0.004812255, rtol=1e-3)
+        assert np.isclose(self.output.flux_index, -0.00576924, rtol=1e-3)
 
     def test_w2weights(self):
         w2 = self.output.w2weights(123.0, 2.6e28, 4 * np.pi)