km3buu/output.py

@@ -20,7 +20,7 @@ from os.path import isfile, join, abspath
 from tempfile import TemporaryDirectory
 import awkward as ak
 import uproot
-from scipy.interpolate import UnivariateSpline
+from scipy.interpolate import UnivariateSpline, interp1d
 from scipy.spatial.transform import Rotation
 import scipy.constants as constants
 import mendeleev
@@ -208,6 +208,34 @@ class GiBUUOutput:
         xsec = np.divide(total_events * weights, n_files)
         return xsec
 
+    @property
+    def mean_xsec(self):
+        root_tupledata = self.arrays
+        energies = np.array(root_tupledata.lepIn_E)
+        weights = self._event_xsec(root_tupledata)
+        Emin = np.min(energies)
+        Emax = np.max(energies)
+        xsec, energy_bins = np.histogram(energies,
+                                         weights=weights,
+                                         bins=np.logspace(
+                                             np.log10(Emin), np.log10(Emax),
+                                             15))
+        deltaE = np.mean(self.flux_data["energy"][1:] -
+                         self.flux_data["energy"][:-1])
+        bin_events = np.array([
+            self.flux_interpolation.integral(energy_bins[i],
+                                             energy_bins[i + 1]) / deltaE
+            for i in range(len(energy_bins) - 1)
+        ])
+        x = (energy_bins[1:] + energy_bins[:-1]) / 2
+        y = xsec / bin_events / x
+        xsec_interp = interp1d(x,
+                               y,
+                               kind="linear",
+                               fill_value=(y[0], y[-1]),
+                               bounds_error=False)
+        return lambda e: xsec_interp(e) * e
+
     def w2weights(self, volume, target_density, solid_angle):
         """
         Calculate w2weights

km3buu/tests/test_jobcard.py

@@ -48,8 +48,8 @@ class TestJobcard(unittest.TestCase):
 class TestNeutrinoJobcard(unittest.TestCase):
     def setUp(self):
         self.test_fluxfile = TemporaryFile()
-        self.test_Z = np.random.randint(0, 100)
-        self.test_A = np.random.randint(0, 100)
+        self.test_Z = np.random.randint(1, 100)
+        self.test_A = np.random.randint(self.test_Z, 100)
         self.test_energy_min = np.random.uniform(0.0, 100.0)
         self.test_energy_max = np.random.uniform(self.test_energy_min, 100.0)
         self.photon_propagation_flag = np.random.choice([True, False])