diff --git a/km3buu/output.py b/km3buu/output.py
index b69462047ee694fb44b2ebf00f52c5f946232507..1310d15b5dfd311576f9f2ae3f628badb8d1ab0e 100644
--- a/km3buu/output.py
+++ b/km3buu/output.py
@@ -258,33 +258,34 @@ class GiBUUOutput:
 
     @property
     def mean_xsec(self):
-        if self.flux_data is None:
-            return lambda energy: self.energy_max
-        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
+        if self.flux_data is None
+            return lambda energy: self.xsection["sum"]
+        else
+            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):
         """