diff --git a/docs/changes/250.api.rst b/docs/changes/250.api.rst
new file mode 100644
index 000000000..97a5b1d44
--- /dev/null
+++ b/docs/changes/250.api.rst
@@ -0,0 +1,6 @@
+In prior versions of pyirf, the energy dispersion matrix was normalized to a
+sum of 1 over the migration axis.
+This is wrong, the correct normalization is to an integral of 1, which is fixed now.
+
+The internal API of the interpolation functions had to be adapted to take in additional
+keywords, mainly the bin edges and the kind of normalization (standard or solid angle cone sections).
diff --git a/pyirf/interpolation/__init__.py b/pyirf/interpolation/__init__.py
index 47d71cb84..2863650b5 100644
--- a/pyirf/interpolation/__init__.py
+++ b/pyirf/interpolation/__init__.py
@@ -11,6 +11,7 @@
     BaseInterpolator,
     DiscretePDFInterpolator,
     ParametrizedInterpolator,
+    PDFNormalization,
 )
 from .component_estimators import (
     BaseComponentEstimator,
@@ -39,6 +40,7 @@
     "BaseInterpolator",
     "BaseNearestNeighborSearcher",
     "BaseExtrapolator",
+    "PDFNormalization",
     "DiscretePDFExtrapolator",
     "ParametrizedExtrapolator",
     "DiscretePDFComponentEstimator",
diff --git a/pyirf/interpolation/base_extrapolators.py b/pyirf/interpolation/base_extrapolators.py
index 79ff7a79c..28b3b38a8 100644
--- a/pyirf/interpolation/base_extrapolators.py
+++ b/pyirf/interpolation/base_extrapolators.py
@@ -3,6 +3,7 @@
 
 import numpy as np
 from pyirf.binning import bin_center
+from pyirf.interpolation.base_interpolators import PDFNormalization
 
 __all__ = ["BaseExtrapolator", "ParametrizedExtrapolator", "DiscretePDFExtrapolator"]
 
@@ -83,21 +84,22 @@ class DiscretePDFExtrapolator(BaseExtrapolator):
     Derived from pyirf.interpolation.BaseExtrapolator
     """
 
-    def __init__(self, grid_points, bin_edges, bin_contents):
+    def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA):
         """DiscretePDFExtrapolator
 
         Parameters
         ----------
-        grid_points: np.ndarray, shape=(n_points, n_dims)
+        grid_points : np.ndarray, shape=(n_points, n_dims)
             Grid points at which templates exist
-        bin_edges: np.ndarray, shape=(n_bins+1)
+        bin_edges : np.ndarray, shape=(n_bins+1)
             Edges of the data binning
-        bin_content: np.ndarray, shape=(n_points, ..., n_bins)
+        binned_pdf : np.ndarray, shape=(n_points, ..., n_bins)
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points, last dimension has to correspond to number of bins for
             the quantity that should be extrapolated (e.g. the Migra axis for EDisp)
-
+        normalization : PDFNormalization
+            How the PDF is normalized
 
         Note
         ----
@@ -105,6 +107,7 @@ def __init__(self, grid_points, bin_edges, bin_contents):
         """
         super().__init__(grid_points)
 
+        self.normalization = normalization
         self.bin_edges = bin_edges
         self.bin_mids = bin_center(self.bin_edges)
-        self.bin_contents = bin_contents
+        self.binned_pdf = binned_pdf
diff --git a/pyirf/interpolation/base_interpolators.py b/pyirf/interpolation/base_interpolators.py
index 0c53ed2c2..9f47f0c31 100644
--- a/pyirf/interpolation/base_interpolators.py
+++ b/pyirf/interpolation/base_interpolators.py
@@ -1,10 +1,27 @@
 """Base classes for interpolators"""
 from abc import ABCMeta, abstractmethod
+import enum
 
 import numpy as np
-from pyirf.binning import bin_center
 
-__all__ = ["BaseInterpolator", "ParametrizedInterpolator", "DiscretePDFInterpolator"]
+from ..binning import bin_center
+
+__all__ = [
+    "BaseInterpolator",
+    "ParametrizedInterpolator",
+    "DiscretePDFInterpolator",
+    "PDFNormalization",
+]
+
+
+class PDFNormalization(enum.Enum):
+    """How a discrete PDF is normalized"""
+
+    #: PDF is normalized to a "normal" area integral of 1
+    AREA = enum.auto()
+    #: PDF is normalized to 1 over the solid angle integral where the bin
+    #: edges represent the opening angles of cones in radian.
+    CONE_SOLID_ANGLE = enum.auto()
 
 
 class BaseInterpolator(metaclass=ABCMeta):
@@ -84,21 +101,24 @@ class DiscretePDFInterpolator(BaseInterpolator):
     Derived from pyirf.interpolation.BaseInterpolator
     """
 
-    def __init__(self, grid_points, bin_edges, bin_contents):
+    def __init__(
+        self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA
+    ):
         """DiscretePDFInterpolator
 
         Parameters
         ----------
-        grid_points: np.ndarray, shape=(n_points, n_dims)
+        grid_points : np.ndarray, shape=(n_points, n_dims)
             Grid points at which interpolation templates exist
-        bin_edges: np.ndarray, shape=(n_bins+1)
+        bin_edges : np.ndarray, shape=(n_bins+1)
             Edges of the data binning
-        bin_content: np.ndarray, shape=(n_points, ..., n_bins)
+        binned_pdf : np.ndarray, shape=(n_points, ..., n_bins)
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points, last dimension has to correspond to number of bins for
             the quantity that should be interpolated (e.g. the Migra axis for EDisp)
-
+        normalization : PDFNormalization
+            How the PDF is normalized
 
         Note
         ----
@@ -108,4 +128,5 @@ def __init__(self, grid_points, bin_edges, bin_contents):
 
         self.bin_edges = bin_edges
         self.bin_mids = bin_center(self.bin_edges)
-        self.bin_contents = bin_contents
+        self.binned_pdf = binned_pdf
+        self.normalization = normalization
diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py
index c34f0a1f7..dc859c868 100644
--- a/pyirf/interpolation/component_estimators.py
+++ b/pyirf/interpolation/component_estimators.py
@@ -7,7 +7,11 @@
 from scipy.spatial import Delaunay
 
 from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator
-from .base_interpolators import DiscretePDFInterpolator, ParametrizedInterpolator
+from .base_interpolators import (
+    DiscretePDFInterpolator,
+    PDFNormalization,
+    ParametrizedInterpolator,
+)
 from .griddata_interpolator import GridDataInterpolator
 from .quantile_interpolator import QuantileInterpolator
 
@@ -152,7 +156,7 @@ def __init__(
         self,
         grid_points,
         bin_edges,
-        bin_contents,
+        binned_pdf,
         interpolator_cls=QuantileInterpolator,
         interpolator_kwargs=None,
         extrapolator_cls=None,
@@ -168,7 +172,7 @@ def __init__(
             Grid points at which interpolation templates exist
         bin_edges: np.ndarray, shape=(n_bins+1)
             Common set of bin-edges for all discretized PDFs.
-        bin_contents: np.ndarray, shape=(n_points, ..., n_bins)
+        binned_pdf: np.ndarray, shape=(n_points, ..., n_bins)
             Discretized PDFs for all grid points and arbitrary further dimensions
             (in IRF term e.g. field-of-view offset bins). Actual interpolation dimension,
             meaning the dimensions that contains actual histograms, has to be along
@@ -191,16 +195,16 @@ def __init__(
         TypeError:
             When bin_edges is not a np.ndarray.
         TypeError:
-            When bin_content is not a np.ndarray..
+            When binned_pdf is not a np.ndarray..
         TypeError:
             When interpolator_cls is not a DiscretePDFInterpolator subclass.
         TypeError:
             When extrapolator_cls is not a DiscretePDFExtrapolator subclass.
         ValueError:
-            When number of bins in bin_edges and contents in bin_contents is
+            When number of bins in bin_edges and contents in binned_pdf is
             not matching.
         ValueError:
-            When number of histograms in bin_contents and points in grid_points
+            When number of histograms in binned_pdf and points in grid_points
             is not matching.
 
         Note
@@ -212,28 +216,28 @@ def __init__(
             grid_points,
         )
 
-        if not isinstance(bin_contents, np.ndarray):
-            raise TypeError("Input bin_contents is not a numpy array.")
-        elif self.n_points != bin_contents.shape[0]:
+        if not isinstance(binned_pdf, np.ndarray):
+            raise TypeError("Input binned_pdf is not a numpy array.")
+        elif self.n_points != binned_pdf.shape[0]:
             raise ValueError(
                 f"Shape missmatch, number of grid_points ({self.n_points}) and "
-                f"number of histograms in bin_contents ({bin_contents.shape[0]}) "
+                f"number of histograms in binned_pdf ({binned_pdf.shape[0]}) "
                 "not matching."
             )
         elif not isinstance(bin_edges, np.ndarray):
             raise TypeError("Input bin_edges is not a numpy array.")
-        elif bin_contents.shape[-1] != (bin_edges.shape[0] - 1):
+        elif binned_pdf.shape[-1] != (bin_edges.shape[0] - 1):
             raise ValueError(
                 f"Shape missmatch, bin_edges ({bin_edges.shape[0] - 1} bins) "
-                f"and bin_contents ({bin_contents.shape[-1]} bins) not matching."
+                f"and binned_pdf ({binned_pdf.shape[-1]} bins) not matching."
             )
 
         # Make sure that 1D input is sorted in increasing order
         if self.grid_dim == 1:
             sorting_inds = np.argsort(self.grid_points.squeeze())
-            
+
             self.grid_points = self.grid_points[sorting_inds]
-            bin_contents = bin_contents[sorting_inds]
+            binned_pdf = binned_pdf[sorting_inds]
 
         if interpolator_kwargs is None:
             interpolator_kwargs = {}
@@ -247,7 +251,7 @@ def __init__(
             )
 
         self.interpolator = interpolator_cls(
-            self.grid_points, bin_edges, bin_contents, **interpolator_kwargs
+            self.grid_points, bin_edges, binned_pdf, **interpolator_kwargs
         )
 
         if extrapolator_cls is None:
@@ -258,7 +262,7 @@ def __init__(
             )
         else:
             self.extrapolator = extrapolator_cls(
-                self.grid_points, bin_edges, bin_contents, **extrapolator_kwargs
+                self.grid_points, bin_edges, binned_pdf, **extrapolator_kwargs
             )
 
 
@@ -334,7 +338,7 @@ def __init__(
         # Make sure that 1D input is sorted in increasing order
         if self.grid_dim == 1:
             sorting_inds = np.argsort(self.grid_points.squeeze())
-            
+
             self.grid_points = self.grid_points[sorting_inds]
             params = params[sorting_inds]
 
@@ -349,7 +353,9 @@ def __init__(
                 f"interpolator_cls must be a ParametrizedInterpolator subclass, got {interpolator_cls}"
             )
 
-        self.interpolator = interpolator_cls(self.grid_points, params, **interpolator_kwargs)
+        self.interpolator = interpolator_cls(
+            self.grid_points, params, **interpolator_kwargs
+        )
 
         if extrapolator_cls is None:
             self.extrapolator = None
@@ -596,7 +602,7 @@ def __init__(
         super().__init__(
             grid_points=grid_points,
             bin_edges=migra_bins,
-            bin_contents=np.swapaxes(energy_dispersion, axis, -1),
+            binned_pdf=np.swapaxes(energy_dispersion, axis, -1),
             interpolator_cls=interpolator_cls,
             interpolator_kwargs=interpolator_kwargs,
             extrapolator_cls=extrapolator_cls,
@@ -681,14 +687,23 @@ def __init__(
 
         psf = np.swapaxes(psf, axis, -1)
 
-        # Renormalize along the source offset axis to have a proper PDF
-        self.omegas = np.diff(cone_solid_angle(source_offset_bins))
-        psf_normed = psf * self.omegas
+        if interpolator_kwargs is None:
+            interpolator_kwargs = {}
+
+        if extrapolator_kwargs is None:
+            extrapolator_kwargs = {}
+
+        interpolator_kwargs.setdefault(
+            "normalization", PDFNormalization.CONE_SOLID_ANGLE
+        )
+        extrapolator_kwargs.setdefault(
+            "normalization", PDFNormalization.CONE_SOLID_ANGLE
+        )
 
         super().__init__(
             grid_points=grid_points,
-            bin_edges=source_offset_bins,
-            bin_contents=psf_normed,
+            bin_edges=source_offset_bins.to_value(u.rad),
+            binned_pdf=psf,
             interpolator_cls=interpolator_cls,
             interpolator_kwargs=interpolator_kwargs,
             extrapolator_cls=extrapolator_cls,
@@ -707,7 +722,7 @@ def __call__(self, target_point):
 
         Returns
         -------
-        psf_interp: np.ndarray, shape=(n_points, ..., n_source_offset_bins)
+        psf_interp: u.Quantity[sr-1], shape=(n_points, ..., n_source_offset_bins)
             Interpolated psf table with same shape as input matrices. For PSF_TABLE
             of shape (n_points, n_energy_bins, n_fov_offset_bins, n_source_offset_bins)
 
@@ -716,4 +731,4 @@ def __call__(self, target_point):
         interpolated_psf_normed = super().__call__(target_point)
 
         # Undo normalisation to get a proper PSF and return
-        return np.swapaxes(interpolated_psf_normed / self.omegas, -1, self.axis)
+        return u.Quantity(np.swapaxes(interpolated_psf_normed, -1, self.axis), u.sr**-1, copy=False)
diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py
index 0d69bc802..f561e050d 100644
--- a/pyirf/interpolation/moment_morph_interpolator.py
+++ b/pyirf/interpolation/moment_morph_interpolator.py
@@ -2,23 +2,26 @@
 from pyirf.binning import bin_center, calculate_bin_indices
 from scipy.spatial import Delaunay
 
-from .base_interpolators import DiscretePDFInterpolator
+from .base_interpolators import DiscretePDFInterpolator, PDFNormalization
+from .utils import get_bin_width
 
 __all__ = [
     "MomentMorphInterpolator",
 ]
 
 
-def _estimate_mean_std(bin_edges, bin_contents):
+def _estimate_mean_std(bin_edges, binned_pdf, normalization):
     """
     Function to roughly estimate mean and standard deviation from a histogram.
 
     Parameters
     ----------
     bin_edges: np.ndarray, shape=(M+1)
-        Array of common bin-edges for bin_contents
-    bin_contents: np.ndarray, shape=(N, ..., M)
-        Array of bin-entries, actual
+        Array of common bin-edges for binned_pdf
+    binned_pdf: np.ndarray, shape=(N, ..., M)
+        PDF values from which to compute mean and std
+    normalization : PDFNormalization
+        How the PDF is normalized
 
     Returns
     -------
@@ -29,35 +32,40 @@ def _estimate_mean_std(bin_edges, bin_contents):
         the input template is =/= 0
     """
     # Create an 2darray where the 1darray mids is repeated n_template times
-    mids = np.broadcast_to(bin_center(bin_edges), bin_contents.shape)
+    mids = np.broadcast_to(bin_center(bin_edges), binned_pdf.shape)
+
+    width = get_bin_width(bin_edges, normalization)
+
+    # integrate pdf to get probability in each bin
+    probability = binned_pdf * width
     # Weighted averages to compute mean and std
-    mean = np.average(mids, weights=bin_contents, axis=-1)
+    mean = np.average(mids, weights=probability, axis=-1)
     std = np.sqrt(
-        np.average((mids - mean[..., np.newaxis]) ** 2, weights=bin_contents, axis=-1)
+        np.average((mids - mean[..., np.newaxis]) ** 2, weights=probability, axis=-1)
     )
 
     # Set std to 0.5*width for all those templates that have only one bin =/= 0. In those
     # cases mids-mean = 0 and therefore std = 0. Uses the width of the one bin with
-    # bin_content!=0
+    # binned_pdf!=0
     mask = std == 0
     if np.any(mask):
-        # Create array of bin_widths and clone for each template
-        width = np.diff(bin_edges) / 2
-        width = np.repeat(width[np.newaxis, :], bin_contents.shape[0], axis=0)
-        std[mask] = width[bin_contents[mask, :] != 0]
+        width = np.diff(bin_edges)
+        # std of a uniform distribution inside the bin
+        uniform_std = np.broadcast_to(np.sqrt(1/12) * width, binned_pdf[mask].shape)
+        std[mask] = uniform_std[binned_pdf[mask, :] != 0]
 
     return mean, std
 
 
-def _lookup(bin_edges, bin_contents, x):
+def _lookup(bin_edges, binned_pdf, x):
     """
     Function to return the bin-height at a desired point.
 
     Parameters
     ----------
     bin_edges: np.ndarray, shape=(M+1)
-        Array of common bin-edges for bin_contents
-    bin_contents: np.ndarray, shape=(N, ..., M)
+        Array of common bin-edges for binned_pdf
+    binned_pdf: np.ndarray, shape=(N, ..., M)
         Array of bin-entries, actual
     x: numpy.ndarray, shape=(N, ..., M)
         Array of M points for each input template, where the histogram-value (bin-height) should be found
@@ -79,11 +87,11 @@ def _lookup(bin_edges, bin_contents, x):
         [
             cont[binnr_row]
             for cont, binnr_row in zip(
-                bin_contents.reshape(-1, bin_contents.shape[-1]),
+                binned_pdf.reshape(-1, binned_pdf.shape[-1]),
                 binnr.reshape(-1, binnr.shape[-1]),
             )
         ]
-    ).reshape(bin_contents.shape)
+    ).reshape(binned_pdf.shape)
 
     # Set all invalid bins to 0
     lu[~valid] = 0
@@ -172,18 +180,20 @@ def barycentric_2D_interpolation_coefficients(grid_points, target_point):
     return coefficients
 
 
-def moment_morph_estimation(bin_edges, bin_contents, coefficients):
+def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization):
     """
     Function that wraps up the moment morph procedure [1] adopted for histograms.
 
     Parameters
     ----------
     bin_edges: np.ndarray, shape=(M+1)
-        Array of common bin-edges for bin_contents
-    bin_contents: np.ndarray, shape=(N, ..., M)
+        Array of common bin-edges for binned_pdf
+    binned_pdf: np.ndarray, shape=(N, ..., M)
         Array of bin-entries, actual
     coefficients: np.ndarray, shape=(N)
-        Estimation coefficients for each entry in bin_contents
+        Estimation coefficients for each entry in binned_pdf
+    normalization : PDFNormalization
+        How the PDF is normalized
 
     Returns
     -------
@@ -199,16 +209,16 @@ def moment_morph_estimation(bin_edges, bin_contents, coefficients):
     bin_mids = bin_center(bin_edges)
 
     # Catch all those templates, where at least one template histogram is all zeros.
-    zero_templates = ~np.all(~np.isclose(np.sum(bin_contents, axis=-1), 0), 0)
+    zero_templates = ~np.all(~np.isclose(np.sum(binned_pdf, axis=-1), 0), 0)
 
     # Manipulate those templates so that computations pass without error
-    bin_contents[:, zero_templates] = np.full(len(bin_mids), 1 / len(bin_mids))
+    binned_pdf[:, zero_templates] = np.full(len(bin_mids), 1 / len(bin_mids))
 
     # Estimate mean and std for each input template histogram. First adaption needed to extend
     # the moment morph procedure to histograms
-    mus, sigs = _estimate_mean_std(bin_edges=bin_edges, bin_contents=bin_contents)
+    mus, sigs = _estimate_mean_std(bin_edges=bin_edges, binned_pdf=binned_pdf, normalization=normalization)
     coefficients = coefficients.reshape(
-        bin_contents.shape[0], *np.ones(mus.ndim - 1, "int")
+        binned_pdf.shape[0], *np.ones(mus.ndim - 1, "int")
     )
 
     # Transform mean and std as in eq. (11) and (12) in [1]
@@ -221,7 +231,7 @@ def moment_morph_estimation(bin_edges, bin_contents, coefficients):
     bij = mus - mu_prime * aij
 
     # Transformation as in eq. (13) in [1]
-    mids = np.broadcast_to(bin_mids, bin_contents.shape)
+    mids = np.broadcast_to(bin_mids, binned_pdf.shape)
     transf_mids = aij[..., np.newaxis] * mids + bij[..., np.newaxis]
 
     # Compute the morphed historgram according to eq. (18) in [1]. The function "lookup" "resamples"
@@ -229,7 +239,7 @@ def moment_morph_estimation(bin_edges, bin_contents, coefficients):
     # bin-mid as new value for a whole transformed bin. Second adaption needed to extend
     # the moment morph procedure to histograms, adaptes the behaviour of eq. (16)
 
-    transf_hist = _lookup(bin_edges=bin_edges, bin_contents=bin_contents, x=transf_mids)
+    transf_hist = _lookup(bin_edges=bin_edges, binned_pdf=binned_pdf, x=transf_mids)
 
     f_new = np.sum(
         np.expand_dims(coefficients, -1) * transf_hist * np.expand_dims(aij, -1), axis=0
@@ -240,15 +250,18 @@ def moment_morph_estimation(bin_edges, bin_contents, coefficients):
 
     # Re-Normalize, needed, as the estimation of the std used above is not exact but the result is scaled with
     # the estimated std
-    norm = np.expand_dims(np.sum(f_new, axis=-1), -1)
+    bin_widths = get_bin_width(bin_edges, normalization)
+    norm = np.expand_dims(np.sum(f_new * bin_widths, axis=-1), -1)
 
     return np.divide(f_new, norm, out=np.zeros_like(f_new), where=norm != 0).reshape(
-        1, *bin_contents.shape[1:]
+        1, *binned_pdf.shape[1:]
     )
 
 
 class MomentMorphInterpolator(DiscretePDFInterpolator):
-    def __init__(self, grid_points, bin_edges, bin_contents):
+    def __init__(
+        self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA
+    ):
         """
         Interpolator class using moment morphing.
 
@@ -258,19 +271,21 @@ def __init__(self, grid_points, bin_edges, bin_contents):
             Grid points at which interpolation templates exist. May be one ot two dimensional.
         bin_edges: np.ndarray, shape=(M+1)
             Edges of the data binning
-        bin_content: np.ndarray, shape=(N, ..., M)
+        binned_pdf: np.ndarray, shape=(N, ..., M)
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points. Interpolation dimension, meaning the
             the quantity that should be interpolated (e.g. the Migra axis for EDisp)
             has to be at axis specified by axis-keyword as well as having entries
             corresponding to the number of bins given through bin_edges keyword.
+        normalization : PDFNormalization
+            How the PDF is normalized
 
         Note
         ----
             Also calls pyirf.interpolation.DiscretePDFInterpolator.__init__.
         """
-        super().__init__(grid_points, bin_edges, bin_contents)
+        super().__init__(grid_points, bin_edges, binned_pdf, normalization)
 
         if self.grid_dim == 2:
             self.triangulation = Delaunay(self.grid_points)
@@ -293,8 +308,9 @@ def _interpolate1D(self, target_point):
 
         return moment_morph_estimation(
             bin_edges=self.bin_edges,
-            bin_contents=self.bin_contents[segment_inds],
+            binned_pdf=self.binned_pdf[segment_inds],
             coefficients=coefficients,
+            normalization=self.normalization,
         )
 
     def _interpolate2D(self, target_point):
@@ -312,8 +328,9 @@ def _interpolate2D(self, target_point):
 
         return moment_morph_estimation(
             bin_edges=self.bin_edges,
-            bin_contents=self.bin_contents[simplex_inds],
+            binned_pdf=self.binned_pdf[simplex_inds],
             coefficients=coefficients,
+            normalization=self.normalization,
         )
 
     def interpolate(self, target_point):
diff --git a/pyirf/interpolation/nearest_neighbor_searcher.py b/pyirf/interpolation/nearest_neighbor_searcher.py
index 8bcaaa108..6ddcf5df7 100644
--- a/pyirf/interpolation/nearest_neighbor_searcher.py
+++ b/pyirf/interpolation/nearest_neighbor_searcher.py
@@ -20,7 +20,7 @@ class BaseNearestNeighborSearcher(BaseInterpolator):
     actual Interpolation/Extrapolation
     """
 
-    def __init__(self, grid_points, contents, norm_ord=2):
+    def __init__(self, grid_points, values, norm_ord=2):
         """
         BaseNearestNeighborSearcher
 
@@ -28,8 +28,8 @@ def __init__(self, grid_points, contents, norm_ord=2):
         ----------
         grid_points: np.ndarray, shape=(n_points, n_dims)
             Grid points at which templates exist
-        contents: np.ndarray, shape=(n_points, ...)
-            Corresponding IRF contents at grid_points
+        values: np.ndarray, shape=(n_points, ...)
+            Corresponding IRF values at grid_points
         norm_ord: non-zero int
             Order of the norm which is used to compute the distances,
             passed to numpy.linalg.norm [1]. Defaults to 2,
@@ -44,10 +44,9 @@ def __init__(self, grid_points, contents, norm_ord=2):
         ----
             Also calls pyirf.interpolation.BaseInterpolators.__init__
         """
-
         super().__init__(grid_points)
 
-        self.contents = contents
+        self.values = values
 
         # Test wether norm_ord is a number
         try:
@@ -67,9 +66,8 @@ def __init__(self, grid_points, contents, norm_ord=2):
 
     def interpolate(self, target_point):
         """
-        Takes a grid of IRF contents for a bunch of different parameters
-        and returns the contents at the nearest grid point
-        as seen from the target point.
+        Takes a grid of IRF values for a bunch of different parameters and returns
+        the values at the nearest grid point as seen from the target point.
 
         Parameters
         ----------
@@ -79,11 +77,11 @@ def interpolate(self, target_point):
         Returns
         -------
         content_new: numpy.ndarray, shape=(1,...,M,...)
-            Contents at nearest neighbor
+            values at nearest neighbor
 
         Note
         ----
-            In case of multiple nearest neighbors, the contents corresponding
+            In case of multiple nearest neighbors, the values corresponding
             to the first one are returned.
         """
 
@@ -96,7 +94,7 @@ def interpolate(self, target_point):
 
         index = np.argmin(distances)
 
-        return self.contents[index, :]
+        return self.values[index, :]
 
 
 class DiscretePDFNearestNeighborSearcher(BaseNearestNeighborSearcher):
@@ -106,7 +104,7 @@ class DiscretePDFNearestNeighborSearcher(BaseNearestNeighborSearcher):
     Compatible with discretized PDF IRF component API.
     """
 
-    def __init__(self, grid_points, bin_edges, bin_contents, norm_ord=2):
+    def __init__(self, grid_points, bin_edges, binned_pdf, norm_ord=2):
         """
         NearestNeighborSearcher compatible with discretized PDF IRF components API
 
@@ -116,7 +114,7 @@ def __init__(self, grid_points, bin_edges, bin_contents, norm_ord=2):
             Grid points at which templates exist
         bin_edges: np.ndarray, shape=(n_bins+1)
             Edges of the data binning. Ignored for nearest neighbor searching.
-        bin_content: np.ndarray, shape=(n_points, ..., n_bins)
+        binned_pdf: np.ndarray, shape=(n_points, ..., n_bins)
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points, last dimension has to correspond to number of bins for
@@ -131,9 +129,7 @@ def __init__(self, grid_points, bin_edges, bin_contents, norm_ord=2):
             Also calls pyirf.interpolation.BaseNearestNeighborSearcher.__init__
         """
 
-        super().__init__(
-            grid_points=grid_points, contents=bin_contents, norm_ord=norm_ord
-        )
+        super().__init__(grid_points=grid_points, values=binned_pdf, norm_ord=norm_ord)
 
 
 DiscretePDFInterpolator.register(DiscretePDFNearestNeighborSearcher)
@@ -168,7 +164,7 @@ def __init__(self, grid_points, params, norm_ord=2):
             Also calls pyirf.interpolation.BaseNearestNeighborSearcher.__init__
         """
 
-        super().__init__(grid_points=grid_points, contents=params, norm_ord=norm_ord)
+        super().__init__(grid_points=grid_points, values=params, norm_ord=norm_ord)
 
 
 ParametrizedInterpolator.register(ParametrizedNearestNeighborSearcher)
diff --git a/pyirf/interpolation/nearest_simplex_extrapolator.py b/pyirf/interpolation/nearest_simplex_extrapolator.py
index 321ac6b11..3dc8c7aa8 100644
--- a/pyirf/interpolation/nearest_simplex_extrapolator.py
+++ b/pyirf/interpolation/nearest_simplex_extrapolator.py
@@ -7,13 +7,14 @@
 import numpy as np
 from scipy.spatial import Delaunay
 
-from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator
+
+from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator, PDFNormalization
 from .moment_morph_interpolator import (
     barycentric_2D_interpolation_coefficients,
     linesegment_1D_interpolation_coefficients,
     moment_morph_estimation,
 )
-from .utils import find_nearest_simplex
+from .utils import find_nearest_simplex, get_bin_width
 
 __all__ = [
     "MomentMorphNearestSimplexExtrapolator",
@@ -122,7 +123,7 @@ def extrapolate(self, target_point):
 
 
 class MomentMorphNearestSimplexExtrapolator(DiscretePDFExtrapolator):
-    def __init__(self, grid_points, bin_edges, bin_contents):
+    def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA):
         """
         Extrapolator class extending/reusing parts of Moment Morphing
         by allowing for negative extrapolation coefficients computed
@@ -136,19 +137,21 @@ def __init__(self, grid_points, bin_edges, bin_contents):
             Have to be sorted in accending order for 1D.
         bin_edges: np.ndarray, shape=(M+1)
             Edges of the data binning
-        bin_content: np.ndarray, shape=(N, ..., M)
+        binned_pdf: np.ndarray, shape=(N, ..., M)
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points. Extrapolation dimension, meaning the
             the quantity that should be interpolated (e.g. the Migra axis for EDisp)
             has to be at the last axis as well as having entries
             corresponding to the number of bins given through bin_edges keyword.
+        normalization : PDFNormalization
+            How the PDF is normalized
 
         Note
         ----
             Also calls pyirf.interpolation.DiscretePDFExtrapolator.__init__.
         """
-        super().__init__(grid_points, bin_edges, bin_contents)
+        super().__init__(grid_points, bin_edges, binned_pdf, normalization)
 
         if self.grid_dim == 2:
             self.triangulation = Delaunay(self.grid_points)
@@ -169,8 +172,9 @@ def _extrapolate1D(self, segment_inds, target_point):
 
         return moment_morph_estimation(
             bin_edges=self.bin_edges,
-            bin_contents=self.bin_contents[segment_inds],
+            binned_pdf=self.binned_pdf[segment_inds],
             coefficients=coefficients,
+            normalization=self.normalization,
         )
 
     def _extrapolate2D(self, simplex_inds, target_point):
@@ -185,8 +189,9 @@ def _extrapolate2D(self, simplex_inds, target_point):
 
         return moment_morph_estimation(
             bin_edges=self.bin_edges,
-            bin_contents=self.bin_contents[simplex_inds],
+            binned_pdf=self.binned_pdf[simplex_inds],
             coefficients=coefficients,
+            normalization=self.normalization,
         )
 
     def extrapolate(self, target_point):
@@ -228,9 +233,10 @@ def extrapolate(self, target_point):
 
             # cut off values below 0 and re-normalize
             extrapolant[extrapolant < 0] = 0
-            norm = np.expand_dims(np.sum(extrapolant, axis=-1), -1)
+            bin_width = get_bin_width(self.bin_edges, self.normalization)
+            norm = np.expand_dims(np.sum(extrapolant * bin_width, axis=-1), -1)
             return np.divide(
                 extrapolant, norm, out=np.zeros_like(extrapolant), where=norm != 0
-            ).reshape(1, *self.bin_contents.shape[1:])
+            ).reshape(1, *self.binned_pdf.shape[1:])
         else:
             return extrapolant
diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py
index 7faeb5cb5..890624191 100644
--- a/pyirf/interpolation/quantile_interpolator.py
+++ b/pyirf/interpolation/quantile_interpolator.py
@@ -1,16 +1,18 @@
 import numpy as np
 from scipy.interpolate import griddata, interp1d
 
-from .base_interpolators import DiscretePDFInterpolator
+from .base_interpolators import DiscretePDFInterpolator, PDFNormalization
+from .utils import get_bin_width
 
 __all__ = ["QuantileInterpolator"]
 
 
-def cdf_values(bin_contents):
+def cdf_values(binned_pdf, bin_edges, normalization):
     """
     compute cdf values and assure they are normed to unity
     """
-    cdfs = np.cumsum(bin_contents, axis=-1)
+    bin_widths = get_bin_width(bin_edges, normalization)
+    cdfs = np.cumsum(binned_pdf * bin_widths, axis=-1)
 
     # assure the last cdf value is 1, ignore errors for empty pdfs as they are reset to 0 by nan_to_num
     with np.errstate(invalid="ignore"):
@@ -95,7 +97,7 @@ def pdf_from_ppf(bin_edges, interp_ppfs, quantiles):
         Edges of the bins in which the final pdf should be binned
 
     interp_ppfs: numpy.ndarray, shape=(1,...,L)
-        Corresponding ppf-values for all self.quantiles at the target_point,
+        Corresponding ppf-values for all quantiles at the target_point,
         not to be confused with QunatileInterpolators self.ppfs, the ppfs
         computed from the input distributions.
 
@@ -122,7 +124,7 @@ def interpolate_ppf(xy):
         pdf = xy[len(xy) // 2 :]
         interpolate = interp1d(ppf, pdf, bounds_error=False, fill_value=(0, 0))
         result = np.nan_to_num(interpolate(bin_edges[:-1]))
-        return np.diff(bin_edges) * result
+        return result
 
     # Interpolate pdf samples and evaluate at bin edges, weight with the bin_width to estimate
     # correct bin height via the midpoint rule formulation of the trapezoidal rule
@@ -131,16 +133,17 @@ def interpolate_ppf(xy):
     return pdf_values
 
 
-def norm_pdf(pdf_values):
+def norm_pdf(pdf_values, bin_edges, normalization):
     """
     Normalize binned_pdf to a sum of 1
     """
     norm = np.sum(pdf_values, axis=-1)
+    width = get_bin_width(bin_edges, normalization)
 
     # Norm all binned_pdfs to unity that are not empty
     normed_pdf_values = np.divide(
         pdf_values,
-        norm[..., np.newaxis],
+        norm[..., np.newaxis] * width,
         out=np.zeros_like(pdf_values),
         where=norm[..., np.newaxis] != 0,
     )
@@ -149,28 +152,37 @@ def norm_pdf(pdf_values):
 
 
 class QuantileInterpolator(DiscretePDFInterpolator):
-    def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e-3):
+    def __init__(
+        self,
+        grid_points,
+        bin_edges,
+        binned_pdf,
+        quantile_resolution=1e-3,
+        normalization=PDFNormalization.AREA,
+    ):
         """BinnedInterpolator constructor
 
         Parameters
         ----------
-        grid_points: np.ndarray
+        grid_points : np.ndarray
             Grid points at which interpolation templates exist
-        bin_edges: np.ndarray
+        bin_edges : np.ndarray
             Edges of the data binning
-        bin_contents: np.ndarray
+        binned_pdf : np.ndarray
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points, the last axis has to correspond to number
             of bins for the quantity that should be interpolated
             (e.g. the Migra axis for EDisp)
-        quantile_resolution: float
+        quantile_resolution : float
             Spacing between quantiles
+        normalization : PDFNormalization
+            How the discrete PDF is normalized
 
         Raises
         ------
         ValueError:
-            When last axis in bin_contents and number of bins are not equal.
+            When last axis in binned_pdf and number of bins are not equal.
 
         Note
         ----
@@ -178,7 +190,7 @@ def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e-
             DiscretePDFInterpolator
         """
         # Remember input shape
-        self.input_shape = bin_contents.shape
+        self.input_shape = binned_pdf.shape
 
         if self.input_shape[-1] != len(bin_edges) - 1:
             raise ValueError(
@@ -187,8 +199,8 @@ def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e-
 
         # include 0-bin at first position in each pdf to avoid edge-effects where the CDF would otherwise
         # start at a value != 0, also extend edges with one bin to the left
-        fill_zeros = np.zeros(shape=bin_contents.shape[:-1])[..., np.newaxis]
-        bin_contents = np.concatenate((fill_zeros, bin_contents), axis=-1)
+        fill_zeros = np.zeros(shape=binned_pdf.shape[:-1])[..., np.newaxis]
+        binned_pdf = np.concatenate((fill_zeros, binned_pdf), axis=-1)
         bin_edges = np.append(bin_edges[0] - np.diff(bin_edges)[0], bin_edges)
 
         # compute quantiles from quantile_resolution
@@ -197,11 +209,14 @@ def __init__(self, grid_points, bin_edges, bin_contents, quantile_resolution=1e-
         )
 
         super().__init__(
-            grid_points=grid_points, bin_edges=bin_edges, bin_contents=bin_contents
+            grid_points=grid_points,
+            bin_edges=bin_edges,
+            binned_pdf=binned_pdf,
+            normalization=normalization,
         )
 
         # Compute CDF values
-        self.cdfs = cdf_values(self.bin_contents)
+        self.cdfs = cdf_values(self.binned_pdf, self.bin_edges, self.normalization)
 
         # compute ppf values at quantiles, determine quantile step of [1]
         self.ppfs = ppf_values(self.bin_mids, self.cdfs, self.quantiles)
@@ -242,7 +257,7 @@ def interpolate(self, target_point):
         )[..., 1:]
 
         # Renormalize pdf to sum of 1
-        normed_interpolated_pdfs = norm_pdf(interpolated_pdfs)
+        normed_interpolated_pdfs = norm_pdf(interpolated_pdfs, self.bin_edges[1:], self.normalization)
 
         # Re-swap axes and set all nans to zero
         return np.nan_to_num(normed_interpolated_pdfs).reshape(1, *self.input_shape[1:])
diff --git a/pyirf/interpolation/tests/__init__.py b/pyirf/interpolation/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/pyirf/interpolation/tests/test_base_extrapolators.py b/pyirf/interpolation/tests/test_base_extrapolators.py
index b7997ea4a..9331a26cf 100644
--- a/pyirf/interpolation/tests/test_base_extrapolators.py
+++ b/pyirf/interpolation/tests/test_base_extrapolators.py
@@ -67,19 +67,19 @@ def test_DiscretePDFExtrapolator_instantiation():
     target2D = np.array([[0.25, 0.25]])
 
     bin_edges = np.linspace(-1, 1, 11)
-    bin_content = np.ones(shape=(len(grid_points1D), len(bin_edges) - 1))
+    binned_pdf = np.ones(shape=(len(grid_points1D), len(bin_edges) - 1))
 
     with pytest.raises(TypeError):  # Abstract class, cannot be instantiated
-        DiscretePDFExtrapolator(grid_points1D, bin_edges, bin_content)
+        DiscretePDFExtrapolator(grid_points1D, bin_edges, binned_pdf)
 
     class DummyBinnedExtrapolator(DiscretePDFExtrapolator):
         def extrapolate(self, target_point):
             return 42
 
-    interp1D = DummyBinnedExtrapolator(grid_points1D, bin_edges, bin_content)
+    interp1D = DummyBinnedExtrapolator(grid_points1D, bin_edges, binned_pdf)
     assert interp1D(target1D) == 42
 
-    interp2D = DummyBinnedExtrapolator(grid_points2D, bin_edges, bin_content)
+    interp2D = DummyBinnedExtrapolator(grid_points2D, bin_edges, binned_pdf)
     assert interp2D(target2D) == 42
 
 
diff --git a/pyirf/interpolation/tests/test_base_interpolators.py b/pyirf/interpolation/tests/test_base_interpolators.py
index c9bb63203..b2b319c3a 100644
--- a/pyirf/interpolation/tests/test_base_interpolators.py
+++ b/pyirf/interpolation/tests/test_base_interpolators.py
@@ -62,19 +62,19 @@ def test_DiscretePDFInterpolator_instantiation():
     target2D = np.array([[1.25, 1.25]])
 
     bin_edges = np.linspace(-1, 1, 11)
-    bin_content = np.ones(shape=(len(grid_points1D), len(bin_edges) - 1))
+    binned_pdf = np.ones(shape=(len(grid_points1D), len(bin_edges) - 1))
 
     with pytest.raises(TypeError):  # Abstract class, cannot be instantiated
-        DiscretePDFInterpolator(grid_points1D, bin_edges, bin_content)
+        DiscretePDFInterpolator(grid_points1D, bin_edges, binned_pdf)
 
     class DummyBinnedInterpolator(DiscretePDFInterpolator):
         def interpolate(self, target_point, **kwargs):
             return 42
 
-    interp1D = DummyBinnedInterpolator(grid_points1D, bin_edges, bin_content)
+    interp1D = DummyBinnedInterpolator(grid_points1D, bin_edges, binned_pdf)
     assert interp1D(target1D) == 42
 
-    interp2D = DummyBinnedInterpolator(grid_points2D, bin_edges, bin_content)
+    interp2D = DummyBinnedInterpolator(grid_points2D, bin_edges, binned_pdf)
     assert interp2D(target2D) == 42
 
 
diff --git a/pyirf/interpolation/tests/test_component_estimator_base_classes.py b/pyirf/interpolation/tests/test_component_estimator_base_classes.py
index 083a67fdb..816c2ed35 100644
--- a/pyirf/interpolation/tests/test_component_estimator_base_classes.py
+++ b/pyirf/interpolation/tests/test_component_estimator_base_classes.py
@@ -127,10 +127,12 @@ def __init__(self, grid_points):
     assert estim2D(target2D_inGrid) == 42
 
     estim1D = DummyEstimator(grid_points1D_good)
-    assert estim1D(target1D_outofGrid) == 43
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert estim1D(target1D_outofGrid) == 43
 
     estim2D = DummyEstimator(grid_points2D_good)
-    assert estim2D(target2D_outofGrid) == 43
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert estim2D(target2D_outofGrid) == 43
 
 
 def test_ParametrizedComponentEstimator_checks():
@@ -208,7 +210,8 @@ def extrapolate(self, target_point):
         extrapolator_cls=DummyExtrapolator,
     )
     assert estim(np.array([[1.5]])) == 42
-    assert estim(np.array([[0]])) == 43
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert estim(np.array([[0]])) == 43
 
 
 def test_DiscretePDFComponentEstimator_checks():
@@ -239,9 +242,9 @@ def extrapolate(self, target_point):
 
     grid_points = np.array([1, 2, 3])
     bin_edges = np.linspace(-1, 1, 11)
-    bin_content_good_shape = np.ones(shape=(len(grid_points), len(bin_edges) - 1))
-    bin_content_bad_nhist = np.ones(shape=(len(grid_points) - 1, len(bin_edges) - 1))
-    bin_content_bad_nbins = np.ones(shape=(len(grid_points), len(bin_edges)))
+    binned_pdf_good_shape = np.ones(shape=(len(grid_points), len(bin_edges) - 1))
+    binned_pdf_bad_nhist = np.ones(shape=(len(grid_points) - 1, len(bin_edges) - 1))
+    binned_pdf_bad_nbins = np.ones(shape=(len(grid_points), len(bin_edges)))
 
     with pytest.raises(
         TypeError,
@@ -249,7 +252,7 @@ def extrapolate(self, target_point):
     ):
         DiscretePDFComponentEstimator(
             grid_points=grid_points,
-            bin_contents=bin_content_good_shape,
+            binned_pdf=binned_pdf_good_shape,
             bin_edges=bin_edges,
             interpolator_cls=WrongInterpolator,
         )
@@ -257,15 +260,15 @@ def extrapolate(self, target_point):
     with pytest.raises(TypeError, match="Input bin_edges is not a numpy array."):
         DiscretePDFComponentEstimator(
             grid_points=grid_points,
-            bin_contents=bin_content_good_shape,
+            binned_pdf=binned_pdf_good_shape,
             bin_edges=bin_edges.tolist(),
             interpolator_cls=DummyInterpolator,
         )
 
-    with pytest.raises(TypeError, match="Input bin_contents is not a numpy array."):
+    with pytest.raises(TypeError, match="Input binned_pdf is not a numpy array."):
         DiscretePDFComponentEstimator(
             grid_points=grid_points,
-            bin_contents=bin_content_good_shape.tolist(),
+            binned_pdf=binned_pdf_good_shape.tolist(),
             bin_edges=bin_edges,
             interpolator_cls=DummyInterpolator,
         )
@@ -274,12 +277,12 @@ def extrapolate(self, target_point):
         ValueError,
         match=re.escape(
             "Shape missmatch, number of grid_points (3) and "
-            "number of histograms in bin_contents (2) not matching."
+            "number of histograms in binned_pdf (2) not matching."
         ),
     ):
         DiscretePDFComponentEstimator(
             grid_points=grid_points,
-            bin_contents=bin_content_bad_nhist,
+            binned_pdf=binned_pdf_bad_nhist,
             bin_edges=bin_edges,
             interpolator_cls=DummyInterpolator,
         )
@@ -288,12 +291,12 @@ def extrapolate(self, target_point):
         ValueError,
         match=re.escape(
             "Shape missmatch, bin_edges (10 bins) "
-            "and bin_contents (11 bins) not matching."
+            "and binned_pdf (11 bins) not matching."
         ),
     ):
         DiscretePDFComponentEstimator(
             grid_points=grid_points,
-            bin_contents=bin_content_bad_nbins,
+            binned_pdf=binned_pdf_bad_nbins,
             bin_edges=bin_edges,
             interpolator_cls=DummyInterpolator,
         )
@@ -305,7 +308,7 @@ def extrapolate(self, target_point):
         DiscretePDFComponentEstimator(
             grid_points=grid_points,
             bin_edges=bin_edges,
-            bin_contents=bin_content_good_shape,
+            binned_pdf=binned_pdf_good_shape,
             interpolator_cls=DummyInterpolator,
             extrapolator_cls=WrongExtrapolator,
         )
@@ -313,12 +316,13 @@ def extrapolate(self, target_point):
     estim = DiscretePDFComponentEstimator(
         grid_points=grid_points,
         bin_edges=bin_edges,
-        bin_contents=bin_content_good_shape,
+        binned_pdf=binned_pdf_good_shape,
         interpolator_cls=DummyInterpolator,
         extrapolator_cls=DummyExtrapolator,
     )
     assert estim(np.array([[1.5]])) == 42
-    assert estim(np.array([[0]])) == 43
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert estim(np.array([[0]])) == 43
 
 
 def test_DiscretePDFComponentEstimator_NearestNeighbors():
@@ -331,11 +335,11 @@ def test_DiscretePDFComponentEstimator_NearestNeighbors():
 
     grid_points = np.array([1, 2, 3])
     bin_edges = np.linspace(0, 1, 11)
-    bin_contents = np.array([np.full(10, x) for x in grid_points])
+    binned_pdf = np.array([np.full(10, x) for x in grid_points])
 
     estim = DiscretePDFComponentEstimator(
         grid_points=grid_points,
-        bin_contents=bin_contents,
+        binned_pdf=binned_pdf,
         bin_edges=bin_edges,
         interpolator_cls=DiscretePDFNearestNeighborSearcher,
         interpolator_kwargs={"norm_ord": 2},
@@ -343,8 +347,9 @@ def test_DiscretePDFComponentEstimator_NearestNeighbors():
         extrapolator_kwargs={"norm_ord": 1},
     )
 
-    assert np.allclose(estim(target_point=np.array([1.1])), bin_contents[0, :])
-    assert np.allclose(estim(target_point=np.array([4.1])), bin_contents[2, :])
+    assert np.allclose(estim(target_point=np.array([1.1])), binned_pdf[0, :])
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert np.allclose(estim(target_point=np.array([4.1])), binned_pdf[2, :])
 
     with pytest.raises(
         TypeError,
@@ -352,7 +357,7 @@ def test_DiscretePDFComponentEstimator_NearestNeighbors():
     ):
         DiscretePDFComponentEstimator(
             grid_points=grid_points,
-            bin_contents=bin_contents,
+            binned_pdf=binned_pdf,
             bin_edges=bin_edges,
             interpolator_cls=ParametrizedNearestNeighborSearcher,
             interpolator_kwargs={"norm_ord": 3},
@@ -382,7 +387,8 @@ def test_ParametrizedComponentEstimator_NearestNeighbors():
     )
 
     assert np.allclose(estim(target_point=np.array([1.1])), params[0, :])
-    assert np.allclose(estim(target_point=np.array([4.1])), params[2, :])
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert np.allclose(estim(target_point=np.array([4.1])), params[2, :])
 
     with pytest.raises(
         TypeError,
@@ -406,7 +412,7 @@ def test_DiscretePDFComponentEstimator_1Dsorting():
 
     grid_points = np.array([[3], [1], [2]])
     bin_edges = np.linspace(0, 1, 11)
-    bin_contents = np.array([np.full(10, x) for x in grid_points])
+    binned_pdf = np.array([np.full(10, x) for x in grid_points])
 
     class DummyInterpolator(DiscretePDFInterpolator):
         def interpolate(self, target_point):
@@ -415,22 +421,24 @@ def interpolate(self, target_point):
     class DummyExtrapolator(DiscretePDFExtrapolator):
         def extrapolate(self, target_point):
             if target_point < self.grid_points.min():
-                return self.bin_contents[0]
+                return self.binned_pdf[0]
             elif target_point > self.grid_points.max():
-                return self.bin_contents[-1]
+                return self.binned_pdf[-1]
 
     estim = DiscretePDFComponentEstimator(
         grid_points=grid_points,
-        bin_contents=bin_contents,
+        binned_pdf=binned_pdf,
         bin_edges=bin_edges,
         interpolator_cls=DummyInterpolator,
         extrapolator_cls=DummyExtrapolator,
     )
 
-    # Nearest neighbor is grid_point 1 at the index 1 of the original bin_contents
-    assert np.allclose(estim(target_point=np.array([0])), bin_contents[1, :])
-    # Nearest neighbor is grid_point 3 at the index 0 of the original bin_contents
-    assert np.allclose(estim(target_point=np.array([4])), bin_contents[0, :])
+    # Nearest neighbor is grid_point 1 at the index 1 of the original binned_pdf
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert np.allclose(estim(target_point=np.array([0])), binned_pdf[1, :])
+    # Nearest neighbor is grid_point 3 at the index 0 of the original binned_pdf
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert np.allclose(estim(target_point=np.array([4])), binned_pdf[0, :])
 
 
 def test_ParametrizedComponentEstimator_1Dsorting():
@@ -460,7 +468,9 @@ def extrapolate(self, target_point):
         extrapolator_cls=DummyExtrapolator,
     )
 
-    # Nearest neighbor is grid_point 1 at the index 1 of the original bin_contents
-    assert np.allclose(estim(target_point=np.array([0])), params[1, :])
-    # Nearest neighbor is grid_point 3 at the index 0 of the original bin_contents
-    assert np.allclose(estim(target_point=np.array([4])), params[0, :])
+    # Nearest neighbor is grid_point 1 at the index 1 of the original binned_pdf
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert np.allclose(estim(target_point=np.array([0])), params[1, :])
+    # Nearest neighbor is grid_point 3 at the index 0 of the original binned_pdf
+    with pytest.warns(UserWarning, match="has to be extrapolated"):
+        assert np.allclose(estim(target_point=np.array([4])), params[0, :])
diff --git a/pyirf/interpolation/tests/test_component_estimator_specific_classes.py b/pyirf/interpolation/tests/test_component_estimator_specific_classes.py
index 49531e4ce..c1f77913b 100644
--- a/pyirf/interpolation/tests/test_component_estimator_specific_classes.py
+++ b/pyirf/interpolation/tests/test_component_estimator_specific_classes.py
@@ -10,6 +10,7 @@ def test_EnergyDispersionEstimator(prod5_irfs):
     zen_pnt = np.array([key.value for key in prod5_irfs.keys()])
     edisps = np.array([irf["edisp"].data for irf in prod5_irfs.values()])
     bin_edges = list(prod5_irfs.values())[0]["edisp"].axes["migra"].edges
+    bin_width = np.diff(bin_edges)
 
     estimator = EnergyDispersionEstimator(
         grid_points=zen_pnt[[0, 2]],
@@ -24,13 +25,12 @@ def test_EnergyDispersionEstimator(prod5_irfs):
 
     interp = estimator(target_point=zen_pnt[[1]])
 
-    assert np.max(interp) <= 1
     assert np.min(interp) >= 0
     assert np.all(np.isfinite(interp))
     assert np.all(
         np.logical_or(
-            np.isclose(np.sum(interp, axis=-2), 1),
-            np.isclose(np.sum(interp, axis=-2), 0),
+            np.isclose(np.sum(interp * bin_width[:, np.newaxis], axis=-2), 1),
+            np.isclose(np.sum(interp * bin_width[:, np.newaxis], axis=-2), 0),
         )
     )
     assert interp.shape == edisps[[1]].shape
@@ -74,15 +74,15 @@ def hist(pnt):
 
     interp = estimator(target_point=zen_pnt[[1]])
 
-    interp *= omegas[np.newaxis, np.newaxis, np.newaxis, ...]
+    probability = (interp * omegas[np.newaxis, np.newaxis, np.newaxis, ...]).to_value(u.one)
 
-    assert np.max(interp) <= 1
-    assert np.min(interp) >= 0
+    assert np.max(probability) <= 1
+    assert np.min(probability) >= 0
     assert np.all(np.isfinite(interp))
     assert np.all(
         np.logical_or(
-            np.isclose(np.sum(interp, axis=-1), 1),
-            np.isclose(np.sum(interp, axis=-1), 0),
+            np.isclose(np.sum(probability, axis=-1), 1),
+            np.isclose(np.sum(probability, axis=-1), 0),
         )
     )
     assert interp.shape == dummy_psfs[[1]].shape
diff --git a/pyirf/interpolation/tests/test_moment_morph_interpolator.py b/pyirf/interpolation/tests/test_moment_morph_interpolator.py
index 605f4310a..77915ccb5 100644
--- a/pyirf/interpolation/tests/test_moment_morph_interpolator.py
+++ b/pyirf/interpolation/tests/test_moment_morph_interpolator.py
@@ -2,6 +2,8 @@
 import pytest
 from scipy.stats import norm
 
+from pyirf.interpolation.base_interpolators import PDFNormalization
+
 
 def expected_mean(a, b):
     return 5 + (a / 5) + (b / 15)
@@ -11,13 +13,13 @@ def expected_std(a, b):
     return 1 + 0.05 * (a + b)
 
 
-def binned_normal_pdf(mean_std_args, bins):
-    pdf = np.diff(
-        norm(loc=expected_mean(*mean_std_args), scale=expected_std(*mean_std_args)).cdf(
-            bins
-        )
+def binned_normal_pdf(a, b, bins):
+    dist = norm(
+        loc=expected_mean(a, b),
+        scale=expected_std(a, b),
     )
-    return pdf / pdf.sum()
+    pdf = np.diff(dist.cdf(bins))
+    return pdf / np.diff(bins)
 
 
 @pytest.fixture
@@ -29,14 +31,14 @@ def bins():
 def simple_1D_data(bins):
     grid = np.array([[20], [40]])
     target = np.array([30])
-    bin_contents = np.array([binned_normal_pdf([x, 0], bins) for x in grid])
+    binned_pdf = np.array([binned_normal_pdf(x, 0, bins) for x in grid])
 
-    truth = binned_normal_pdf([target, 0], bins)
+    truth = binned_normal_pdf(target, 0, bins)
 
     return {
         "grid": grid,
         "target": target,
-        "bin_contents": bin_contents,
+        "binned_pdf": binned_pdf,
         "truth": truth,
     }
 
@@ -45,14 +47,14 @@ def simple_1D_data(bins):
 def simple_2D_data(bins):
     grid = np.array([[20, 20], [60, 20], [40, 60]])
     target = np.array([25, 25])
-    bin_contents = np.array([binned_normal_pdf(x, bins) for x in grid])
+    binned_pdf = np.array([binned_normal_pdf(*x, bins) for x in grid])
 
-    truth = binned_normal_pdf(target, bins)
+    truth = binned_normal_pdf(*target, bins)
 
     return {
         "grid": grid,
         "target": target,
-        "bin_contents": bin_contents,
+        "binned_pdf": binned_pdf,
         "truth": truth,
     }
 
@@ -61,13 +63,13 @@ def test_estimate_mean_std(bins):
     from pyirf.interpolation.moment_morph_interpolator import _estimate_mean_std
 
     grid = np.array([[20], [40]])
-    bin_contents = np.array(
+    binned_pdf = np.array(
         [
             [
-                [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)],
+                [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 1, 0, bins)],
                 [
-                    binned_normal_pdf([x + 1.5, 0], bins),
-                    binned_normal_pdf([x + 10, 0], bins),
+                    binned_normal_pdf(x + 1.5, 0, bins),
+                    binned_normal_pdf(x + 10, 0, bins),
                 ],
             ]
             for x in grid
@@ -94,18 +96,18 @@ def test_estimate_mean_std(bins):
         ]
     ).squeeze()
 
-    mean, std = _estimate_mean_std(bins, bin_contents)
+    mean, std = _estimate_mean_std(bins, binned_pdf, PDFNormalization.AREA)
 
     # Assert estimation and truth within one bin
-    assert np.allclose(mean, true_mean, atol=np.diff(bins)[0] / 2)
-    assert np.allclose(std, true_std, atol=np.diff(bins)[0] / 2)
+    np.testing.assert_allclose(mean, true_mean, atol=np.diff(bins)[0] / 2)
+    np.testing.assert_allclose(std, true_std, atol=np.diff(bins)[0] / 2)
 
 
 def test_lookup():
     from pyirf.interpolation.moment_morph_interpolator import _lookup
 
     bins = np.array([0, 0.1, 0.2, 0.3, 0.4, 0.5])
-    bin_contents = np.array(
+    binned_pdf = np.array(
         [
             [[1, 2, 3, 4, 5], [6, 7, 8, 9, 10]],
             [[6, 7, 8, 9, 10], [1, 2, 3, 4, 5]],
@@ -126,7 +128,7 @@ def test_lookup():
         ]
     )
 
-    assert np.allclose(_lookup(bins, bin_contents, x), truth)
+    np.testing.assert_allclose(_lookup(bins, binned_pdf, x), truth)
 
 
 def test_linesegment_1D_interpolation_coefficients():
@@ -137,16 +139,16 @@ def test_linesegment_1D_interpolation_coefficients():
     grid_points = np.array([[20], [40]])
     target_point = np.array([[20]])
 
-    assert np.allclose(
+    np.testing.assert_allclose(
         linesegment_1D_interpolation_coefficients(grid_points, target_point),
-        np.array([1, 0]),
+        np.array([[1, 0]]),
     )
 
     target_point = np.array([[40]])
 
-    assert np.allclose(
+    np.testing.assert_allclose(
         linesegment_1D_interpolation_coefficients(grid_points, target_point),
-        np.array([0, 1]),
+        np.array([[0, 1]]),
     )
 
     target_point = np.array([[25]])
@@ -156,11 +158,11 @@ def test_linesegment_1D_interpolation_coefficients():
     mfrac = (target_point[0, 0] - grid_points[0, 0]) / (
         grid_points[1, 0] - grid_points[0, 0]
     )
-    cs = np.array([1 - mfrac, mfrac])
+    cs = np.array([[1 - mfrac, mfrac]])
 
     res = linesegment_1D_interpolation_coefficients(grid_points, target_point)
 
-    assert np.allclose(res, cs)
+    np.testing.assert_allclose(res, cs)
     assert np.sum(res) == 1
     assert np.all(np.logical_and(res < 1, res > 0))
 
@@ -173,28 +175,28 @@ def test_barycentric_2D_interpolation_coefficients():
     grid_points = np.array([[20, 20], [60, 20], [40, 60]])
     target_point = np.array([20, 20])
 
-    assert np.allclose(
+    np.testing.assert_allclose(
         barycentric_2D_interpolation_coefficients(grid_points, target_point),
         np.array([1, 0, 0]),
     )
 
     target_point = np.array([60, 20])
 
-    assert np.allclose(
+    np.testing.assert_allclose(
         barycentric_2D_interpolation_coefficients(grid_points, target_point),
         np.array([0, 1, 0]),
     )
 
     target_point = np.array([40, 60])
 
-    assert np.allclose(
+    np.testing.assert_allclose(
         barycentric_2D_interpolation_coefficients(grid_points, target_point),
         np.array([0, 0, 1]),
     )
 
     target_point = np.array([40, 20])
 
-    assert np.allclose(
+    np.testing.assert_allclose(
         barycentric_2D_interpolation_coefficients(grid_points, target_point),
         np.array([0.5, 0.5, 0]),
     )
@@ -203,7 +205,7 @@ def test_barycentric_2D_interpolation_coefficients():
     target_point = np.array([40, 100 / 3])
     res = barycentric_2D_interpolation_coefficients(grid_points, target_point)
 
-    assert np.allclose(res, np.array([1, 1, 1]) / 3)
+    np.testing.assert_allclose(res, np.array([1, 1, 1]) / 3)
 
 
 def test_moment_morph_estimation1D(bins, simple_1D_data):
@@ -212,16 +214,16 @@ def test_moment_morph_estimation1D(bins, simple_1D_data):
         moment_morph_estimation,
     )
 
-    grid, target, bin_contents, truth = simple_1D_data.values()
+    grid, target, binned_pdf, truth = simple_1D_data.values()
 
     coeffs = linesegment_1D_interpolation_coefficients(grid, target)
-    res = moment_morph_estimation(bins, bin_contents, coeffs)
+    res = moment_morph_estimation(bins, binned_pdf, coeffs, PDFNormalization.AREA)
 
-    assert np.isclose(np.sum(res), 1)
+    np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1)
     assert np.all(np.isfinite(res))
     assert res.shape == (1, len(bins) - 1)
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_moment_morph_estimation2D(bins, simple_2D_data):
@@ -230,34 +232,34 @@ def test_moment_morph_estimation2D(bins, simple_2D_data):
         moment_morph_estimation,
     )
 
-    grid, target, bin_contents, truth = simple_2D_data.values()
+    grid, target, binned_pdf, truth = simple_2D_data.values()
 
     coeffs = barycentric_2D_interpolation_coefficients(grid, target)
-    res = moment_morph_estimation(bins, bin_contents, coeffs)
+    res = moment_morph_estimation(bins, binned_pdf, coeffs, PDFNormalization.AREA)
 
-    assert np.isclose(np.sum(res), 1)
+    np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1)
     assert np.all(np.isfinite(res))
     assert res.shape == (1, len(bins) - 1)
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphInterpolator1D(bins, simple_1D_data):
     from pyirf.interpolation import MomentMorphInterpolator
 
-    grid, target, bin_contents, truth = simple_1D_data.values()
+    grid, target, binned_pdf, truth = simple_1D_data.values()
 
     interp = MomentMorphInterpolator(
-        grid_points=grid, bin_edges=bins, bin_contents=bin_contents
+        grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf
     )
 
     res = interp(target)
 
-    assert np.isclose(np.sum(res), 1)
+    np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1)
     assert np.all(np.isfinite(res))
     assert res.shape == (1, len(bins) - 1)
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphInterpolator1D_dirac_delta_input():
@@ -265,7 +267,7 @@ def test_MomentMorphInterpolator1D_dirac_delta_input():
 
     grid = np.array([[1], [3]])
     bin_edges = np.array([0, 1, 2, 3, 4])
-    bin_contents = np.array(
+    binned_pdf = np.array(
         [
             [[0, 1, 0, 0], [0.25, 0.25, 0.25, 0.25]],
             [[0, 0, 0, 1], [0.25, 0.25, 0.25, 0.25]],
@@ -273,41 +275,42 @@ def test_MomentMorphInterpolator1D_dirac_delta_input():
     )
     target = np.array([2])
 
-    interp = MomentMorphInterpolator(grid, bin_edges, bin_contents)
+    interp = MomentMorphInterpolator(grid, bin_edges, binned_pdf)
     res = interp(target)
 
-    assert np.allclose(res, np.array([[0, 0, 1, 0], [0.25, 0.25, 0.25, 0.25]]))
+    expected = np.array([[[0, 0, 1, 0], [0.25, 0.25, 0.25, 0.25]]])
+    np.testing.assert_allclose(res, expected)
 
 
 def test_MomentMorphInterpolator1D_all_empty(bins, simple_1D_data):
     from pyirf.interpolation import MomentMorphInterpolator
 
     grid, target, _, _ = simple_1D_data.values()
-    bin_contents = np.array([np.zeros(len(bins) - 1) for _ in grid])
+    binned_pdf = np.array([np.zeros(len(bins) - 1) for _ in grid])
 
     interp = MomentMorphInterpolator(
-        grid_points=grid, bin_edges=bins, bin_contents=bin_contents
+        grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf
     )
 
     res = interp(target)
 
-    assert np.allclose(res, 0)
+    np.testing.assert_allclose(res, 0)
 
 
 def test_MomentMorphInterpolator1D_partially_empty(bins, simple_1D_data):
     from pyirf.interpolation import MomentMorphInterpolator
 
-    grid, target, bin_contents, _ = simple_1D_data.values()
+    grid, target, binned_pdf, _ = simple_1D_data.values()
 
-    bin_contents[0, :] = np.zeros(len(bins) - 1)
+    binned_pdf[0, :] = np.zeros(len(bins) - 1)
 
     interp = MomentMorphInterpolator(
-        grid_points=grid, bin_edges=bins, bin_contents=bin_contents
+        grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf
     )
 
     res = interp(target)
 
-    assert np.allclose(res, 0)
+    np.testing.assert_allclose(res, 0)
 
 
 def test_MomentMorphInterpolator1D_mixed_data(bins):
@@ -317,34 +320,31 @@ def test_MomentMorphInterpolator1D_mixed_data(bins):
     target = np.array([30])
 
     # Create template histograms
-    bin_contents = np.array(
+    binned_pdf = np.array(
         [
             [
-                [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)],
-                [
-                    binned_normal_pdf([x + 10, 0], bins),
-                    binned_normal_pdf([x + 2, 0], bins),
-                ],
+                [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 1, 0, bins)],
+                [binned_normal_pdf(x + 10, 0, bins), binned_normal_pdf(x + 2, 0, bins)],
             ]
             for x in grid
         ]
     )
 
     # Make template histograms at indizes [:, 1, 1, :] all zeroed
-    bin_contents[:, 1, 1, :] = np.zeros(len(bins) - 1)
+    binned_pdf[:, 1, 1, :] = np.zeros(len(bins) - 1)
 
     # Zero template histogram at index [1, 0, 0, :]
-    bin_contents[1, 0, 0, :] = np.zeros(len(bins) - 1)
+    binned_pdf[1, 0, 0, :] = np.zeros(len(bins) - 1)
 
     truth = np.array(
         [
             [
-                binned_normal_pdf([target, 0], bins),
-                binned_normal_pdf([target + 1, 0], bins),
+                binned_normal_pdf(target, 0, bins),
+                binned_normal_pdf(target + 1, 0, bins),
             ],
             [
-                binned_normal_pdf([target + 10, 0], bins),
-                binned_normal_pdf([target + 2, 0], bins),
+                binned_normal_pdf(target + 10, 0, bins),
+                binned_normal_pdf(target + 2, 0, bins),
             ],
         ]
     )
@@ -354,17 +354,17 @@ def test_MomentMorphInterpolator1D_mixed_data(bins):
     truth[1, 1, :] = np.zeros(len(bins) - 1)
 
     interp = MomentMorphInterpolator(
-        grid_points=grid, bin_edges=bins, bin_contents=bin_contents
+        grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf
     )
 
     res = interp(target)
 
-    expected_norms = np.array([[0, 1], [1, 0]])
-    assert np.allclose(np.sum(res, axis=-1), expected_norms)
+    expected_norms = np.array([[[0, 1], [1, 0]]])
+    np.testing.assert_allclose(np.sum(res * np.diff(bins), axis=-1), expected_norms)
     assert np.all(np.isfinite(res))
-    assert res.shape == (1, *bin_contents.shape[1:])
+    assert res.shape == (1, *binned_pdf.shape[1:])
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphInterpolator1D_extended_grid_extradims(bins):
@@ -372,13 +372,13 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins):
 
     grid = np.array([[20], [40], [60], [80]])
     target = np.array([25])
-    bin_contents = np.array(
+    binned_pdf = np.array(
         [
             [
-                [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)],
+                [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 1, 0, bins)],
                 [
-                    binned_normal_pdf([x + 10, 0], bins),
-                    binned_normal_pdf([x + 2, 0], bins),
+                    binned_normal_pdf(x + 10, 0, bins),
+                    binned_normal_pdf(x + 2, 0, bins),
                 ],
             ]
             for x in grid
@@ -386,78 +386,78 @@ def test_MomentMorphInterpolator1D_extended_grid_extradims(bins):
     )
 
     interp = MomentMorphInterpolator(
-        grid_points=grid, bin_edges=bins, bin_contents=bin_contents
+        grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf
     )
 
     truth = np.array(
         [
             [
-                binned_normal_pdf([target, 0], bins),
-                binned_normal_pdf([target + 1, 0], bins),
+                binned_normal_pdf(target, 0, bins),
+                binned_normal_pdf(target + 1, 0, bins),
             ],
             [
-                binned_normal_pdf([target + 10, 0], bins),
-                binned_normal_pdf([target + 2, 0], bins),
+                binned_normal_pdf(target + 10, 0, bins),
+                binned_normal_pdf(target + 2, 0, bins),
             ],
         ]
     )
 
     res = interp(target)
 
-    assert np.allclose(np.sum(res, axis=-1), 1)
+    np.testing.assert_allclose(np.sum(res * np.diff(bins), axis=-1), 1)
     assert np.all(np.isfinite(res))
-    assert res.shape == (1, *bin_contents.shape[1:])
+    assert res.shape == (1, *binned_pdf.shape[1:])
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-4, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphInterpolator2D(bins, simple_2D_data):
     from pyirf.interpolation import MomentMorphInterpolator
 
-    grid, target, bin_contents, truth = simple_2D_data.values()
+    grid, target, binned_pdf, truth = simple_2D_data.values()
 
     interp = MomentMorphInterpolator(
-        grid_points=grid, bin_edges=bins, bin_contents=bin_contents
+        grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf
     )
 
     res = interp(target)
 
-    assert np.isclose(np.sum(res), 1)
+    np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1)
     assert np.all(np.isfinite(res))
     assert res.shape == (1, len(bins) - 1)
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphInterpolator2D_partially_empty(bins, simple_2D_data):
     from pyirf.interpolation import MomentMorphInterpolator
 
-    grid, target, bin_contents, _ = simple_2D_data.values()
+    grid, target, binned_pdf, _ = simple_2D_data.values()
 
-    bin_contents[0, :] = np.zeros(len(bins) - 1)
+    binned_pdf[0, :] = np.zeros(len(bins) - 1)
 
     interp = MomentMorphInterpolator(
-        grid_points=grid, bin_edges=bins, bin_contents=bin_contents
+        grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf
     )
 
     res = interp(target)
 
-    assert np.allclose(res, 0)
+    np.testing.assert_allclose(res, 0)
 
 
 def test_MomentMorphInterpolator2D_all_empty(bins, simple_2D_data):
     from pyirf.interpolation import MomentMorphInterpolator
 
     grid, target, _, _ = simple_2D_data.values()
-    bin_contents = np.array([np.zeros(len(bins) - 1) for _ in grid])
+    binned_pdf = np.array([np.zeros(len(bins) - 1) for _ in grid])
 
     interp = MomentMorphInterpolator(
-        grid_points=grid, bin_edges=bins, bin_contents=bin_contents
+        grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf
     )
 
     res = interp(target)
 
-    assert np.allclose(res, 0)
+    np.testing.assert_allclose(res, 0)
 
 
 def test_MomentMorphInterpolator2D_mixed(bins):
@@ -467,34 +467,31 @@ def test_MomentMorphInterpolator2D_mixed(bins):
     target = np.array([25, 25])
 
     # Create template histograms
-    bin_contents = np.array(
+    binned_pdf = np.array(
         [
             [
-                [binned_normal_pdf(x, bins), binned_normal_pdf([x[0] + 1, x[1]], bins)],
-                [
-                    binned_normal_pdf([x[0] + 10, x[1]], bins),
-                    binned_normal_pdf([x[0] + 2, x[1]], bins),
-                ],
+                [binned_normal_pdf(a, b, bins), binned_normal_pdf(a + 1, b, bins)],
+                [binned_normal_pdf(a + 10, b, bins), binned_normal_pdf(a + 2, b, bins)],
             ]
-            for x in grid
+            for a, b in grid
         ]
     )
 
     # Make template histograms at indizes [:, 1, 1, :] all zeroed
-    bin_contents[:, 1, 1, :] = np.zeros(len(bins) - 1)
+    binned_pdf[:, 1, 1, :] = np.zeros(len(bins) - 1)
 
     # Zero template histogram at index [1, 0, 0, :]
-    bin_contents[1, 0, 0, :] = np.zeros(len(bins) - 1)
+    binned_pdf[1, 0, 0, :] = np.zeros(len(bins) - 1)
 
     truth = np.array(
         [
             [
-                binned_normal_pdf(target, bins),
-                binned_normal_pdf([target[0] + 1, target[1]], bins),
+                binned_normal_pdf(*target, bins),
+                binned_normal_pdf(target[0] + 1, target[1], bins),
             ],
             [
-                binned_normal_pdf([target[0] + 10, target[1]], bins),
-                binned_normal_pdf([target[0] + 2, target[1]], bins),
+                binned_normal_pdf(target[0] + 10, target[1], bins),
+                binned_normal_pdf(target[0] + 2, target[1], bins),
             ],
         ]
     )
@@ -504,17 +501,17 @@ def test_MomentMorphInterpolator2D_mixed(bins):
     truth[1, 1, :] = np.zeros(len(bins) - 1)
 
     interp = MomentMorphInterpolator(
-        grid_points=grid, bin_edges=bins, bin_contents=bin_contents
+        grid_points=grid, bin_edges=bins, binned_pdf=binned_pdf
     )
 
     res = interp(target)
 
-    expected_norms = np.array([[0, 1], [1, 0]])
-    assert np.allclose(np.sum(res, axis=-1), expected_norms)
+    expected_norms = np.array([[[0, 1], [1, 0]]])
+    np.testing.assert_allclose(np.sum(res * np.diff(bins), axis=-1), expected_norms)
     assert np.all(np.isfinite(res))
-    assert res.shape == (1, *bin_contents.shape[1:])
+    assert res.shape == (1, *binned_pdf.shape[1:])
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphInterpolator1D_extended_grid(bins):
@@ -522,22 +519,22 @@ def test_MomentMorphInterpolator1D_extended_grid(bins):
 
     grid = np.array([[20], [40], [60], [80]])
     target = np.array([25])
-    bin_contents = np.array([binned_normal_pdf([x, 0], bins) for x in grid])
+    binned_pdf = np.array([binned_normal_pdf(x, 0, bins) for x in grid])
 
     interp = MomentMorphInterpolator(
         grid_points=grid,
         bin_edges=bins,
-        bin_contents=bin_contents,
+        binned_pdf=binned_pdf,
     )
 
     res = interp(target)
-    truth = binned_normal_pdf([target, 0], bins)
+    truth = binned_normal_pdf(target, 0, bins)
 
-    assert np.isclose(np.sum(res), 1)
+    np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1)
     assert np.all(np.isfinite(res))
     assert res.shape == (1, len(bins) - 1)
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphInterpolator2D_extended_grid(bins):
@@ -545,68 +542,70 @@ def test_MomentMorphInterpolator2D_extended_grid(bins):
 
     grid = np.array([[20, 20], [40, 20], [30, 40], [50, 20], [45, 40]])
     target = np.array([25, 25])
-    bin_contents = np.array([binned_normal_pdf(x, bins) for x in grid])
+    binned_pdf = np.array([binned_normal_pdf(a, b, bins) for a, b in grid])
 
     interp = MomentMorphInterpolator(
         grid_points=grid,
         bin_edges=bins,
-        bin_contents=bin_contents,
+        binned_pdf=binned_pdf,
     )
 
     res = interp(target)
-    truth = binned_normal_pdf(target, bins)
+    truth = binned_normal_pdf(*target, bins)
 
-    assert np.isclose(np.sum(res), 1)
+    np.testing.assert_almost_equal(np.sum(res * np.diff(bins)), 1)
     assert np.all(np.isfinite(res))
     assert res.shape == (1, len(bins) - 1)
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphInterpolator2D_extended_grid_extradims(bins):
     from pyirf.interpolation import MomentMorphInterpolator
 
     grid = np.array([[20, 20], [40, 20], [30, 40], [50, 20], [45, 40]])
-    target = np.array([25, 25])
-    bin_contents = np.array(
+    a = 25
+    b = 25
+    target = np.array([a, b])
+    binned_pdf = np.array(
         [
             [
-                [binned_normal_pdf(x, bins), binned_normal_pdf([x[0] + 1, x[1]], bins)],
+                [binned_normal_pdf(a, b, bins), binned_normal_pdf(a + 1, b, bins)],
                 [
-                    binned_normal_pdf([x[0] + 10, x[1]], bins),
-                    binned_normal_pdf([x[0] + 2, x[1]+5], bins),
+                    binned_normal_pdf(a + 10, b, bins),
+                    binned_normal_pdf(a + 2, b + 5, bins),
                 ],
             ]
-            for x in grid
+            for a, b in grid
         ]
     )
 
     interp = MomentMorphInterpolator(
         grid_points=grid,
         bin_edges=bins,
-        bin_contents=bin_contents,
+        binned_pdf=binned_pdf,
     )
 
     truth = np.array(
         [
             [
-                binned_normal_pdf(target, bins),
-                binned_normal_pdf([target[0] + 1, target[1]], bins),
+                binned_normal_pdf(a, b, bins),
+                binned_normal_pdf(a + 1, b, bins),
             ],
             [
-                binned_normal_pdf([target[0] + 10, target[1]], bins),
-                binned_normal_pdf([target[0] + 2, target[1]+5], bins),
+                binned_normal_pdf(a + 10, b, bins),
+                binned_normal_pdf(a + 2, b + 5, bins),
             ],
         ]
     )
 
     res = interp(target)
 
-    assert np.allclose(np.sum(res, axis=-1), 1)
+    np.testing.assert_allclose(np.sum(res * np.diff(bins), axis=-1), 1)
     assert np.all(np.isfinite(res))
-    assert res.shape == (1, *bin_contents.shape[1:])
+    assert res.shape == (1, *binned_pdf.shape[1:])
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res.squeeze(), truth, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphInterpolator3D():
@@ -616,7 +615,7 @@ def test_MomentMorphInterpolator3D():
 
     grid = np.array([[0, 0, 0], [0, 20, 0], [20, 0, 0], [20, 20, 0], [10, 10, 10]])
 
-    bin_contents = np.array([np.ones(len(bins) - 1) / (len(bins) - 1) for _ in grid])
+    binned_pdf = np.array([np.ones(len(bins) - 1) / (len(bins) - 1) for _ in grid])
 
     with pytest.raises(
         NotImplementedError,
@@ -625,5 +624,5 @@ def test_MomentMorphInterpolator3D():
         MomentMorphInterpolator(
             grid_points=grid,
             bin_edges=bins,
-            bin_contents=bin_contents,
+            binned_pdf=binned_pdf,
         )
diff --git a/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py b/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py
index 66c25f45d..6a3fc57bf 100644
--- a/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py
+++ b/pyirf/interpolation/tests/test_nearest_neighbor_searcher.py
@@ -13,7 +13,7 @@ def grid_2d():
 
 
 @pytest.fixture
-def contents(grid_1d):
+def binned_pdf(grid_1d):
     return np.array(
         [
             [[np.full(10, x), np.full(10, x / 2)], [np.full(10, 2 * x), np.zeros(10)]]
@@ -22,83 +22,83 @@ def contents(grid_1d):
     )
 
 
-def test_BaseNearestNeighborSearcher_1DGrid(grid_1d, contents):
+def test_BaseNearestNeighborSearcher_1DGrid(grid_1d, binned_pdf):
     from pyirf.interpolation import BaseNearestNeighborSearcher
 
-    searcher = BaseNearestNeighborSearcher(grid_1d, contents, norm_ord=2)
+    searcher = BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord=2)
 
     target = np.array([0])
-    assert np.array_equal(searcher(target), contents[0, :])
+    assert np.array_equal(searcher(target), binned_pdf[0, :])
 
     target = np.array([[1.9]])
-    assert np.array_equal(searcher(target), contents[1, :])
+    assert np.array_equal(searcher(target), binned_pdf[1, :])
 
 
-def test_BaseNearestNeighborSearcher_2DGrid(grid_2d, contents):
+def test_BaseNearestNeighborSearcher_2DGrid(grid_2d, binned_pdf):
     from pyirf.interpolation import BaseNearestNeighborSearcher
 
-    searcher = BaseNearestNeighborSearcher(grid_2d, contents, norm_ord=2)
+    searcher = BaseNearestNeighborSearcher(grid_2d, binned_pdf, norm_ord=2)
 
     target = np.array([[0, 1]])
-    assert np.array_equal(searcher(target), contents[0, :])
+    assert np.array_equal(searcher(target), binned_pdf[0, :])
 
     target = np.array([3, 3])
-    assert np.array_equal(searcher(target), contents[-1, :])
+    assert np.array_equal(searcher(target), binned_pdf[-1, :])
 
 
-def test_BaseNearestNeighborSearcher_manhatten_norm(grid_2d, contents):
+def test_BaseNearestNeighborSearcher_manhatten_norm(grid_2d, binned_pdf):
     from pyirf.interpolation import BaseNearestNeighborSearcher
 
-    searcher = BaseNearestNeighborSearcher(grid_2d, contents, norm_ord=1)
+    searcher = BaseNearestNeighborSearcher(grid_2d, binned_pdf, norm_ord=1)
 
     target = np.array([[0, 1]])
-    assert np.array_equal(searcher(target), contents[0, :])
+    assert np.array_equal(searcher(target), binned_pdf[0, :])
 
     target = np.array([[3, 3]])
-    assert np.array_equal(searcher(target), contents[-1, :])
+    assert np.array_equal(searcher(target), binned_pdf[-1, :])
 
 
-def test_BaseNearestNeighborSearcher_wrong_norm(grid_1d, contents):
+def test_BaseNearestNeighborSearcher_wrong_norm(grid_1d, binned_pdf):
     from pyirf.interpolation import BaseNearestNeighborSearcher
 
     with pytest.raises(ValueError, match="Only positiv integers allowed for norm_ord"):
-        BaseNearestNeighborSearcher(grid_1d, contents, norm_ord=-2)
+        BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord=-2)
 
     with pytest.raises(ValueError, match="Only positiv integers allowed for norm_ord"):
-        BaseNearestNeighborSearcher(grid_1d, contents, norm_ord=1.5)
+        BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord=1.5)
 
     with pytest.raises(ValueError, match="Only positiv integers allowed for norm_ord"):
-        BaseNearestNeighborSearcher(grid_1d, contents, norm_ord=np.inf)
+        BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord=np.inf)
 
     with pytest.raises(ValueError, match="Only positiv integers allowed for norm_ord"):
-        BaseNearestNeighborSearcher(grid_1d, contents, norm_ord="nuc")
+        BaseNearestNeighborSearcher(grid_1d, binned_pdf, norm_ord="nuc")
 
 
-def test_DiscretePDFNearestNeighborSearcher(grid_2d, contents):
+def test_DiscretePDFNearestNeighborSearcher(grid_2d, binned_pdf):
     from pyirf.interpolation import DiscretePDFNearestNeighborSearcher
 
-    bin_edges = np.linspace(0, 1, contents.shape[-1] + 1)
+    bin_edges = np.linspace(0, 1, binned_pdf.shape[-1] + 1)
 
     searcher = DiscretePDFNearestNeighborSearcher(
-        grid_points=grid_2d, bin_edges=bin_edges, bin_contents=contents, norm_ord=1
+        grid_points=grid_2d, bin_edges=bin_edges, binned_pdf=binned_pdf, norm_ord=1
     )
 
     target = np.array([[0, 1]])
-    assert np.array_equal(searcher(target), contents[0, :])
+    assert np.array_equal(searcher(target), binned_pdf[0, :])
 
     target = np.array([[3, 3]])
-    assert np.array_equal(searcher(target), contents[-1, :])
+    assert np.array_equal(searcher(target), binned_pdf[-1, :])
 
 
-def test_ParametrizedNearestNeighborSearcher(grid_2d, contents):
+def test_ParametrizedNearestNeighborSearcher(grid_2d, binned_pdf):
     from pyirf.interpolation import ParametrizedNearestNeighborSearcher
 
     searcher = ParametrizedNearestNeighborSearcher(
-        grid_points=grid_2d, params=contents, norm_ord=1
+        grid_points=grid_2d, params=binned_pdf, norm_ord=1
     )
 
     target = np.array([[0, 1]])
-    assert np.array_equal(searcher(target), contents[0, :])
+    assert np.array_equal(searcher(target), binned_pdf[0, :])
 
     target = np.array([[3, 3]])
-    assert np.array_equal(searcher(target), contents[-1, :])
+    assert np.array_equal(searcher(target), binned_pdf[-1, :])
diff --git a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py
index 2ffc5bd8c..313cdc727 100644
--- a/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py
+++ b/pyirf/interpolation/tests/test_nearest_simplex_extrapolator.py
@@ -11,13 +11,13 @@ def expected_std(a, b):
     return 1 + 0.05 * (a + b)
 
 
-def binned_normal_pdf(mean_std_args, bins):
-    pdf = np.diff(
-        norm(loc=expected_mean(*mean_std_args), scale=expected_std(*mean_std_args)).cdf(
-            bins
-        )
+def binned_normal_pdf(a, b, bins):
+    dist = norm(
+        loc=expected_mean(a, b),
+        scale=expected_std(a, b),
     )
-    return pdf / pdf.sum()
+    pdf = np.diff(dist.cdf(bins))
+    return pdf / np.diff(bins)
 
 
 @pytest.fixture
@@ -45,17 +45,17 @@ def test_ParametrizedNearestSimplexExtrapolator_1DGrid():
     target_point1 = np.array([3])
     interpolant1 = interpolator(target_point1)
 
-    dummy_data_target1 = 3 * slope + 1
+    dummy_data_target1 = (3 * slope + 1)[np.newaxis, :]
 
-    assert np.allclose(interpolant1, dummy_data_target1)
+    np.testing.assert_allclose(interpolant1, dummy_data_target1)
     assert interpolant1.shape == (1, *dummy_data.shape[1:])
 
     target_point2 = np.array([[-2.5]])
     interpolant2 = interpolator(target_point2)
 
-    dummy_data_target2 = -2.5 * slope + 1
+    dummy_data_target2 = (-2.5 * slope + 1)[np.newaxis, :]
 
-    assert np.allclose(interpolant2, dummy_data_target2)
+    np.testing.assert_allclose(interpolant2, dummy_data_target2)
     assert interpolant2.shape == (1, *dummy_data.shape[1:])
 
 
@@ -96,7 +96,7 @@ def test_ParametrizedNearestSimplexExtrapolator_2DGrid():
         ]
     ).squeeze()
 
-    assert np.allclose(interpolant1.squeeze(), dummy_data_target1)
+    np.testing.assert_allclose(interpolant1.squeeze(), dummy_data_target1)
     assert interpolant1.shape == (1, *dummy_data.shape[1:])
 
     target_point2 = np.array([[-2.5, -5.5]])
@@ -104,12 +104,14 @@ def test_ParametrizedNearestSimplexExtrapolator_2DGrid():
 
     dummy_data_target2 = np.array(
         [
-            np.dot((m.T * target_point2 + n), np.array([1, 1]))
-            for m, n in zip(slope, intercept)
+            [
+                np.dot((m.T * target_point2 + n), np.array([1, 1]))
+                for m, n in zip(slope, intercept)
+            ]
         ]
-    ).squeeze()
+    )
 
-    assert np.allclose(interpolant2, dummy_data_target2)
+    np.testing.assert_allclose(interpolant2, dummy_data_target2)
     assert interpolant2.shape == (1, *dummy_data.shape[1:])
 
 
@@ -137,32 +139,34 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins):
 
     grid_points = np.array([[20], [30], [40]])
 
+    n_bins = len(bins) - 1
+    bin_width = np.diff(bins)
+
     # Create template histograms
-    bin_contents = np.array(
+    binned_pdf = np.array(
         [
             [
-                [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 1, 0], bins)],
+                [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 1, 0, bins)],
                 [
-                    binned_normal_pdf([x + 10, 0], bins),
-                    np.zeros(len(bins) - 1),
+                    binned_normal_pdf(x + 10, 0, bins),
+                    np.zeros(n_bins),
                 ],
-                [np.zeros(len(bins) - 1), np.ones(len(bins) - 1) / (len(bins) - 1)],
+                [np.zeros(n_bins), np.ones(n_bins) / n_bins / bin_width],
             ]
             for x in grid_points
         ]
     )
 
     # Include dirac-delta moving with grid_point x at [:, 2, 0, x]
-    bin_contents[0, 2, 0, 20] = 1
-    bin_contents[1, 2, 0, 30] = 1
-    bin_contents[2, 2, 0, 40] = 1
+    for i, bin_idx in enumerate(grid_points.ravel()):
+        binned_pdf[i, 2, 0, bin_idx] = 1 / bin_width[bin_idx]
 
     # Zero template histogram at index [0, 0, 0, :], extrapolations "lower"
     # then this bin has consequently to be zero in these bins.
-    bin_contents[0, 0, 0, :] = np.zeros(len(bins) - 1)
+    binned_pdf[0, 0, 0, :] = 0.0
 
     extrap = MomentMorphNearestSimplexExtrapolator(
-        grid_points=grid_points, bin_contents=bin_contents, bin_edges=bins
+        grid_points=grid_points, binned_pdf=binned_pdf, bin_edges=bins
     )
 
     target1 = np.array([10])
@@ -171,84 +175,86 @@ def test_MomentMorphNearestSimplexExtrapolator_1DGrid(bins):
     truth1 = np.array(
         [
             [
-                np.zeros(len(bins) - 1),
-                binned_normal_pdf([target1 + 1, 0], bins),
+                np.zeros(n_bins),
+                binned_normal_pdf(target1 + 1, 0, bins),
             ],
             [
-                binned_normal_pdf([target1 + 10, 0], bins),
-                np.zeros(len(bins) - 1),
+                binned_normal_pdf(target1 + 10, 0, bins),
+                np.zeros_like(bin_width),
             ],
-            [np.zeros(len(bins) - 1), np.ones(len(bins) - 1) / (len(bins) - 1)],
+            [np.zeros(n_bins), np.ones(n_bins) / n_bins / bin_width],
         ]
     )
 
     # Inlcude dirac-delta
-    truth1[2, 0, 10] = 1
+    delta = np.zeros(len(bins) - 1)
+    delta[target1[0]] = 1
+    delta /= bin_width
+    truth1[2, 0] = delta
 
     res1 = extrap(target1)
 
-    expected_norms1 = np.array([[0, 1], [1, 0], [1, 1]])
-    assert np.allclose(np.sum(res1, axis=-1), expected_norms1)
+    expected_norms1 = np.array([[[0, 1], [1, 0], [1, 1]]])
+    np.testing.assert_allclose(np.sum(res1 * np.diff(bins), axis=-1), expected_norms1)
     assert np.all(np.isfinite(res1))
-    assert res1.shape == (1, *bin_contents.shape[1:])
+    assert res1.shape == (1, *binned_pdf.shape[1:])
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res1.squeeze(), truth1, atol=4e-3, rtol=1e-5)
 
     target2 = np.array([45])
 
     truth2 = np.array(
         [
             [
-                binned_normal_pdf([target2, 0], bins),
-                binned_normal_pdf([target2 + 1, 0], bins),
+                binned_normal_pdf(target2, 0, bins),
+                binned_normal_pdf(target2 + 1, 0, bins),
             ],
             [
-                binned_normal_pdf([target2 + 10, 0], bins),
-                np.zeros(len(bins) - 1),
+                binned_normal_pdf(target2 + 10, 0, bins),
+                np.zeros(n_bins),
             ],
-            [np.zeros(len(bins) - 1), np.ones(len(bins) - 1) / (len(bins) - 1)],
+            [np.zeros(n_bins), np.ones(n_bins) / n_bins / bin_width],
         ]
     )
 
     # Inlcude dirac-delta
-    truth2[2, 0, 45] = 1
+    truth2[2, 0, 45] = 1 / bin_width[45]
 
     res2 = extrap(target2)
 
-    expected_norms2 = np.array([[1, 1], [1, 0], [1, 1]])
-    assert np.allclose(np.sum(res2, axis=-1), expected_norms2)
+    expected_norms2 = np.array([[[1, 1], [1, 0], [1, 1]]])
+    np.testing.assert_allclose(np.sum(res2 * bin_width, axis=-1), expected_norms2)
     assert np.all(np.isfinite(res2))
-    assert res2.shape == (1, *bin_contents.shape[1:])
+    assert res2.shape == (1, *binned_pdf.shape[1:])
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins):
     """Test ParametrizedNearestSimplexExtrapolator on a 2D Grid"""
     from pyirf.interpolation import MomentMorphNearestSimplexExtrapolator
 
+    n_bins = len(bins) - 1
+    bin_width = np.diff(bins)
     grid_points = np.array([[20, 20], [30, 10], [40, 20], [50, 10]])
 
     # Create template histograms
-    bin_contents = np.array(
+    binned_pdf = np.array(
         [
             [
-                [binned_normal_pdf(x, bins), binned_normal_pdf([x[0] + 1, x[1]], bins)],
-                [
-                    binned_normal_pdf([x[0] + 10, x[1]], bins),
-                    np.zeros(len(bins) - 1),
-                ],
+                [binned_normal_pdf(a, b, bins), binned_normal_pdf(a + 1, b, bins)],
+                [binned_normal_pdf(a + 10, b, bins), np.zeros(n_bins)],
             ]
-            for x in grid_points
+            for a, b in grid_points
         ]
     )
 
     # Zero template histogram at index [0, 0, 0, :] (the [20, 20] point),
     # so that targets interpolated from left simplex should also be zero at [0, 0, :]
-    bin_contents[0, 0, 0, :] = np.zeros(len(bins) - 1)
+    binned_pdf[0, 0, 0, :] = np.zeros(n_bins)
 
     extrap = MomentMorphNearestSimplexExtrapolator(
-        grid_points=grid_points, bin_contents=bin_contents, bin_edges=bins
+        grid_points=grid_points, binned_pdf=binned_pdf, bin_edges=bins
     )
 
     target1 = np.array([25, 10])
@@ -257,48 +263,48 @@ def test_MomentMorphNearestSimplexExtrapolator_2DGrid(bins):
     truth1 = np.array(
         [
             [
-                np.zeros(len(bins) - 1),
-                binned_normal_pdf([target1[0] + 1, target1[1]], bins),
+                np.zeros(n_bins),
+                binned_normal_pdf(target1[0] + 1, target1[1], bins),
             ],
             [
-                binned_normal_pdf([target1[0] + 10, target1[1]], bins),
-                np.zeros(len(bins) - 1),
+                binned_normal_pdf(target1[0] + 10, target1[1], bins),
+                np.zeros(n_bins),
             ],
         ]
     )
 
     res1 = extrap(target1)
 
-    expected_norms1 = np.array([[0, 1], [1, 0]])
-    assert np.allclose(np.sum(res1, axis=-1), expected_norms1)
+    expected_norms1 = np.array([[[0, 1], [1, 0]]])
+    np.testing.assert_allclose(np.sum(res1 * bin_width, axis=-1), expected_norms1)
     assert np.all(np.isfinite(res1))
-    assert res1.shape == (1, *bin_contents.shape[1:])
+    assert res1.shape == (1, *binned_pdf.shape[1:])
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res1.squeeze(), truth1, atol=1e-3, rtol=1e-5)
 
     target2 = np.array([45, 20])
 
     truth2 = np.array(
         [
             [
-                binned_normal_pdf(target2, bins),
-                binned_normal_pdf([target2[0] + 1, target2[1]], bins),
+                binned_normal_pdf(*target2, bins),
+                binned_normal_pdf(target2[0] + 1, target2[1], bins),
             ],
             [
-                binned_normal_pdf([target2[0] + 10, target2[1]], bins),
-                np.zeros(len(bins) - 1),
+                binned_normal_pdf(target2[0] + 10, target2[1], bins),
+                np.zeros(n_bins),
             ],
         ]
     )
 
     res2 = extrap(target2)
 
-    expected_norms2 = np.array([[1, 1], [1, 0]])
-    assert np.allclose(np.sum(res2, axis=-1), expected_norms2)
+    expected_norms2 = np.array([[[1, 1], [1, 0]]])
+    np.testing.assert_allclose(np.sum(res2 * bin_width, axis=-1), expected_norms2)
     assert np.all(np.isfinite(res2))
-    assert res2.shape == (1, *bin_contents.shape[1:])
+    assert res2.shape == (1, *binned_pdf.shape[1:])
     # Assert truth and result matching within +- 0.1%, atol dominates comparison
-    assert np.allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5)
+    np.testing.assert_allclose(res2.squeeze(), truth2, atol=1e-3, rtol=1e-5)
 
 
 def test_MomentMorphNearestSimplexExtrapolator_3DGrid():
@@ -315,7 +321,7 @@ def test_MomentMorphNearestSimplexExtrapolator_3DGrid():
         match="Extrapolation in more then two dimension not impemented.",
     ):
         MomentMorphNearestSimplexExtrapolator(
-            grid_points=grid_points, bin_contents=dummy_data, bin_edges=bin_edges
+            grid_points=grid_points, binned_pdf=dummy_data, bin_edges=bin_edges
         )
 
 
@@ -330,15 +336,15 @@ def test_MomentMorphNearestSimplexExtrapolator_below_zero_warning(bins):
     grid_points = np.array([[30], [40]])
     bins = np.linspace(-10, 40, 51)
 
-    bin_contents = np.array(
+    binned_pdf = np.array(
         [
-            [binned_normal_pdf([x, 0], bins), binned_normal_pdf([x + 10, 0], bins)]
+            [binned_normal_pdf(x, 0, bins), binned_normal_pdf(x + 10, 0, bins)]
             for x in grid_points
         ]
     )
 
     extrap = MomentMorphNearestSimplexExtrapolator(
-        grid_points=grid_points, bin_contents=bin_contents, bin_edges=bins
+        grid_points=grid_points, binned_pdf=binned_pdf, bin_edges=bins
     )
 
     with pytest.warns(
@@ -347,5 +353,5 @@ def test_MomentMorphNearestSimplexExtrapolator_below_zero_warning(bins):
     ):
         res = extrap(np.array([0]))
 
-    assert np.allclose(np.sum(res, axis=-1), 1)
+    np.testing.assert_allclose(np.sum(res, axis=-1), 1)
     assert np.all(res >= 0)
diff --git a/pyirf/interpolation/tests/test_quantile_interpolator.py b/pyirf/interpolation/tests/test_quantile_interpolator.py
index 195d95fdc..2a1892cd0 100644
--- a/pyirf/interpolation/tests/test_quantile_interpolator.py
+++ b/pyirf/interpolation/tests/test_quantile_interpolator.py
@@ -3,6 +3,8 @@
 from pyirf.binning import bin_center
 from scipy.stats import norm
 
+from pyirf.interpolation.base_interpolators import PDFNormalization
+
 
 @pytest.fixture
 def data():
@@ -12,6 +14,8 @@ def data():
 
     # create binned pdfs by interpolation of bin content
     binned_pdfs = np.array([np.diff(dist.cdf(bin_edges)) for dist in distributions])
+    bin_width = np.diff(bin_edges)
+    binned_pdfs /= bin_width[np.newaxis, :]
 
     dataset = {
         "bin_edges": bin_edges,
@@ -28,11 +32,19 @@ def data():
 def test_cdf_values(data):
     from pyirf.interpolation.quantile_interpolator import cdf_values
 
-    cdf_est = cdf_values(data["binned_pdfs"][0])
-
     # Assert empty histograms result in cdf containing zeros
-    assert np.all(cdf_values(np.zeros(shape=5)) == 0)
+    np.testing.assert_array_equal(
+        cdf_values(
+            np.zeros_like(data["binned_pdfs"])[0],
+            data["bin_edges"],
+            PDFNormalization.AREA,
+        ),
+        0,
+    )
 
+    cdf_est = cdf_values(
+        data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA
+    )
     # Assert cdf is increasing or constant for actual pdfs
     assert np.all(np.diff(cdf_est) >= 0)
 
@@ -40,7 +52,8 @@ def test_cdf_values(data):
     assert np.max(cdf_est) == 1
 
     # Assert estimated and true cdf are matching
-    assert np.allclose(cdf_est, data["distributions"][0].cdf(data["bin_edges"][1:]))
+    true_cdf = data["distributions"][0].cdf(data["bin_edges"][1:])
+    np.testing.assert_allclose(cdf_est, true_cdf, atol=1e-12)
 
 
 def test_ppf_values(data):
@@ -54,11 +67,13 @@ def test_ppf_values(data):
     bin_mids = bin_center(data["bin_edges"])
 
     # Estimated ppf-values
-    cdf_est = cdf_values(data["binned_pdfs"][0])
+    cdf_est = cdf_values(
+        data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA
+    )
     ppf_est = ppf_values(bin_mids, cdf_est, quantiles)
 
     # Assert truth and estimation match allowing for +- bin_width deviation
-    assert np.allclose(ppf_true, ppf_est, atol=np.diff(data["bin_edges"])[0])
+    np.testing.assert_allclose(ppf_true, ppf_est, atol=np.diff(data["bin_edges"])[0])
 
 
 def test_pdf_from_ppf(data):
@@ -73,21 +88,38 @@ def test_pdf_from_ppf(data):
     bin_mids = bin_center(data["bin_edges"])
 
     # Estimate ppf-values
-    cdf_est = cdf_values(data["binned_pdfs"][0])
+    cdf_est = cdf_values(
+        data["binned_pdfs"][0], data["bin_edges"], PDFNormalization.AREA
+    )
     ppf_est = ppf_values(bin_mids, cdf_est, quantiles)
 
     # Compute pdf_values
     pdf_est = pdf_from_ppf(data["bin_edges"], ppf_est, quantiles)
 
     # Assert pdf-values matching true pdf within +-1%
-    assert np.allclose(pdf_est, data["binned_pdfs"][0], atol=1e-2)
+    np.testing.assert_allclose(pdf_est, data["binned_pdfs"][0], atol=1e-2)
 
 
 def test_norm_pdf(data):
     from pyirf.interpolation.quantile_interpolator import norm_pdf
 
-    assert np.allclose(norm_pdf(2 * data["binned_pdfs"][0]), data["binned_pdfs"][0])
-    assert np.allclose(norm_pdf(np.zeros(5)), 0)
+    result = norm_pdf(
+        2 * data["binned_pdfs"][0],
+        data["bin_edges"],
+        PDFNormalization.AREA,
+    )
+    np.testing.assert_allclose(np.sum(result * np.diff(data["bin_edges"])), 1)
+    np.testing.assert_allclose(
+        result,
+        data["binned_pdfs"][0],
+    )
+
+    result = norm_pdf(
+        np.zeros(len(data["bin_edges"]) - 1),
+        data["bin_edges"],
+        PDFNormalization.AREA,
+    )
+    np.testing.assert_allclose(result, 0)
 
 
 def test_interpolate_binned_pdf(data):
@@ -96,7 +128,7 @@ def test_interpolate_binned_pdf(data):
     interpolator = QuantileInterpolator(
         grid_points=data["grid_points"][[0, 2]],
         bin_edges=data["bin_edges"],
-        bin_contents=data["binned_pdfs"][[0, 2], :],
+        binned_pdf=data["binned_pdfs"][[0, 2], :],
         quantile_resolution=1e-3,
     )
 
diff --git a/pyirf/interpolation/tests/test_utils.py b/pyirf/interpolation/tests/test_utils.py
index 60fe5f7df..cba06d6cf 100644
--- a/pyirf/interpolation/tests/test_utils.py
+++ b/pyirf/interpolation/tests/test_utils.py
@@ -147,3 +147,15 @@ def test_find_nearest_simplex(non_rect_grid):
     assert find_nearest_simplex(non_rect_grid, np.array([-10, -10])) == 0
     assert find_nearest_simplex(non_rect_grid, np.array([10, 30])) == 1
     assert find_nearest_simplex(non_rect_grid, np.array([20.00000000001, -10])) == 2
+
+
+def test_get_bin_width():
+    from pyirf.interpolation.utils import get_bin_width
+    from pyirf.interpolation import PDFNormalization
+
+    bins = np.array([0, 1, 3])
+    np.testing.assert_allclose(get_bin_width(bins, PDFNormalization.AREA), [1, 2])
+
+    bins = np.array([0, np.pi / 3, np.pi / 2])
+    width = get_bin_width(bins, PDFNormalization.CONE_SOLID_ANGLE)
+    np.testing.assert_allclose(width, [np.pi, np.pi])
diff --git a/pyirf/interpolation/utils.py b/pyirf/interpolation/utils.py
index aa54ad100..5936ec63a 100644
--- a/pyirf/interpolation/utils.py
+++ b/pyirf/interpolation/utils.py
@@ -1,5 +1,19 @@
+import astropy.units as u
 import numpy as np
 
+from ..utils import cone_solid_angle
+from .base_interpolators import PDFNormalization
+
+
+def get_bin_width(bin_edges, normalization):
+    if normalization is PDFNormalization.AREA:
+        return np.diff(bin_edges)
+
+    if normalization is PDFNormalization.CONE_SOLID_ANGLE:
+        return np.diff(cone_solid_angle(bin_edges).to_value(u.sr))
+
+    raise ValueError(f"Invalid PDF normalization: {normalization}")
+
 
 def plumb_point_dist(line, target):
     """
diff --git a/pyirf/irf/energy_dispersion.py b/pyirf/irf/energy_dispersion.py
index 9bb203999..de631e426 100644
--- a/pyirf/irf/energy_dispersion.py
+++ b/pyirf/irf/energy_dispersion.py
@@ -10,9 +10,10 @@
 ]
 
 
-def _normalize_hist(hist):
+def _normalize_hist(hist, migration_bins):
     # make sure we do not mutate the input array
     hist = hist.copy()
+    bin_width = np.diff(migration_bins)
 
     # calculate number of events along the N_MIGRA axis to get events
     # per energy per fov
@@ -21,8 +22,8 @@ def _normalize_hist(hist):
     with np.errstate(invalid="ignore"):
         # hist shape is (N_E, N_MIGRA, N_FOV), norm shape is (N_E, N_FOV)
         # so we need to add a new axis in the middle to get (N_E, 1, N_FOV)
-        # for broadcasting
-        hist = hist / norm[:, np.newaxis, :]
+        # bin_width is 1d, so we need newaxis, use the values, newaxis
+        hist = hist / norm[:, np.newaxis, :] / bin_width[np.newaxis, :, np.newaxis]
 
     return np.nan_to_num(hist)
 
@@ -74,8 +75,7 @@ def energy_dispersion(
         ],
     )
 
-    n_events_per_energy = energy_dispersion.sum(axis=1)
-    energy_dispersion = _normalize_hist(energy_dispersion)
+    energy_dispersion = _normalize_hist(energy_dispersion, migration_bins)
 
     return energy_dispersion
 
diff --git a/pyirf/irf/tests/__init__.py b/pyirf/irf/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/pyirf/irf/tests/test_energy_dispersion.py b/pyirf/irf/tests/test_energy_dispersion.py
index b7dba84a5..036d106de 100644
--- a/pyirf/irf/tests/test_energy_dispersion.py
+++ b/pyirf/irf/tests/test_energy_dispersion.py
@@ -50,35 +50,31 @@ def test_energy_dispersion():
     )
 
     assert result.shape == (3, 1000, 2)
-    assert np.isclose(result.sum(), 6.0)
 
-    cumulative_sum = np.cumsum(result, axis=1)
+    bin_width = np.diff(migration_bins)
     bin_centers = 0.5 * (migration_bins[1:] + migration_bins[:-1])
+    # edisp shape is (N_E, N_MIGRA, N_FOV), we need to integrate over migration axis
+    integral = (result * bin_width[np.newaxis, :, np.newaxis]).sum(axis=1)
+    np.testing.assert_allclose(integral, 1.0)
+
+    cdf = np.cumsum(result * bin_width[np.newaxis, :, np.newaxis], axis=1)
+
+    def ppf(cdf, bins, value):
+        return np.interp(value, cdf, bins[1:])
+
     assert np.isclose(
         TRUE_SIGMA_1,
-        (
-            bin_centers[np.where(cumulative_sum[0, :, :] >= 0.84)[0][0]]
-            - bin_centers[np.where(cumulative_sum[0, :, :] >= 0.16)[0][0]]
-        )
-        / 2,
+        0.5 * (ppf(cdf[0, :, 0], migration_bins, 0.84) - ppf(cdf[0, :, 0], migration_bins, 0.16)),
         rtol=0.1,
     )
     assert np.isclose(
         TRUE_SIGMA_2,
-        (
-            bin_centers[np.where(cumulative_sum[1, :, :] >= 0.84)[0][0]]
-            - bin_centers[np.where(cumulative_sum[1, :, :] >= 0.16)[0][0]]
-        )
-        / 2,
+        0.5 * (ppf(cdf[1, :, 0], migration_bins, 0.84) - ppf(cdf[1, :, 0], migration_bins, 0.16)),
         rtol=0.1,
     )
     assert np.isclose(
         TRUE_SIGMA_3,
-        (
-            bin_centers[np.where(cumulative_sum[2, :, :] >= 0.84)[0][0]]
-            - bin_centers[np.where(cumulative_sum[2, :, :] >= 0.16)[0][0]]
-        )
-        / 2,
+        0.5 * (ppf(cdf[2, :, 0], migration_bins, 0.84) - ppf(cdf[2, :, 0], migration_bins, 0.16)),
         rtol=0.1,
     )
 
diff --git a/pyproject.toml b/pyproject.toml
index 51f3b9859..95c0fe5d2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -32,4 +32,11 @@ build-backend = "setuptools.build_meta"
 
     [[tool.towncrier.section]]
         name = ""
-        path = ""
\ No newline at end of file
+        path = ""
+
+
+[tool.pytest.ini_options]
+filterwarnings = [
+    "ignore:pkg_resources is deprecated:DeprecationWarning",
+    "ignore:Deprecated call to `pkg_resources:DeprecationWarning",
+]