diff --git a/skyllh/analyses/i3/publicdata_ps/aeff.py b/skyllh/analyses/i3/publicdata_ps/aeff.py
index 6dd4aa8141..2d185cf659 100644
--- a/skyllh/analyses/i3/publicdata_ps/aeff.py
+++ b/skyllh/analyses/i3/publicdata_ps/aeff.py
@@ -1,8 +1,5 @@
 import numpy as np
-from scipy import (
-    integrate,
-    interpolate,
-)
+from scipy import interpolate
 
 from skyllh.analyses.i3.publicdata_ps.utils import (
     FctSpline2D,
@@ -326,18 +323,13 @@ def get_detection_prob_for_decnu(self, decnu, enu_min, enu_max, enu_range_min, e
 
         spl = interpolate.splrep(x, y, k=1, s=0)
 
-        def _eval_spl_func(x):
-            return interpolate.splev(x, spl, der=0, ext=1)
-
-        norm = integrate.quad(_eval_spl_func, enu_range_min, enu_range_max, limit=200, full_output=1)[0]
+        norm = interpolate.splint(enu_range_min, enu_range_max, spl)
 
         enu_min = np.atleast_1d(enu_min)
         enu_max = np.atleast_1d(enu_max)
 
         det_prob = np.empty((len(enu_min),), dtype=np.double)
         for i in range(len(enu_min)):
-            integral = integrate.quad(_eval_spl_func, enu_min[i], enu_max[i], limit=200, full_output=1)[0]
-
-            det_prob[i] = integral / norm
+            det_prob[i] = interpolate.splint(enu_min[i], enu_max[i], spl) / norm
 
         return det_prob
diff --git a/skyllh/analyses/i3/publicdata_ps/signalpdf.py b/skyllh/analyses/i3/publicdata_ps/signalpdf.py
index c54af80341..b21725a97d 100644
--- a/skyllh/analyses/i3/publicdata_ps/signalpdf.py
+++ b/skyllh/analyses/i3/publicdata_ps/signalpdf.py
@@ -1,5 +1,4 @@
 import numpy as np
-from scipy import integrate
 
 from skyllh.analyses.i3.publicdata_ps.aeff import (
     PDAeff,
@@ -83,16 +82,6 @@ def __init__(
         # Add the PDF axes.
         self.add_axis(PDFAxis(name='log_energy', vmin=self.log10_reco_e_min, vmax=self.log10_reco_e_max))
 
-        # Check integrity.
-        integral = (
-            integrate.quad(
-                self.f_e_spl.evaluate, self.log10_reco_e_min, self.log10_reco_e_max, limit=200, full_output=1
-            )[0]
-            / self.f_e_spl.norm
-        )
-        if not np.isclose(integral, 1):
-            raise ValueError(f'The integral over log10_reco_e of the energy term must be unity! But it is {integral}!')
-
     def assert_is_valid_for_trial_data(self, tdm, tl=None):
         pass
 
diff --git a/skyllh/analyses/i3/publicdata_ps/utils.py b/skyllh/analyses/i3/publicdata_ps/utils.py
index 44ff53d08e..594d90e271 100644
--- a/skyllh/analyses/i3/publicdata_ps/utils.py
+++ b/skyllh/analyses/i3/publicdata_ps/utils.py
@@ -1,8 +1,5 @@
 import numpy as np
-from scipy import (
-    integrate,
-    interpolate,
-)
+from scipy import interpolate
 
 from skyllh.core.binning import (
     get_bincenters_from_binedges,
@@ -42,7 +39,9 @@ class from scipy.
 
         self.norm = None
         if norm:
-            self.norm = integrate.quad(self.__call__, self.x_min, self.x_max, limit=200, full_output=1)[0]
+            # The spline is defined only in the `x` (bincenters) interval by construction.
+            # We choose not to extrapolate out-of-range values.
+            self.norm = float(self.spl_f.integrate(x[0], x[-1]))
 
     def __call__(self, x, oor_value=0):
         """Evaluates the spline at the given x values. For x-values