Merge pull request #268 from cta-observatory/fix_bias_resolution_from…

…_edisp

Fix energy_bias_resolution_from_energy_dispersion

docs/changes/268.bugfix.rst

@@ -0,0 +1,4 @@
+Fix ``pyirf.benchmarks.energy_bias_resolution_from_energy_dispersion``.
+This function was not adapted to the now correct normalization of the
+energy dispersion matrix, resulting in wrong results on the now correct
+matrices.

pyirf/benchmarks/energy_bias_resolution.py

@@ -131,15 +131,18 @@ def energy_bias_resolution_from_energy_dispersion(
         Bin edges for the relative energy migration (``reco_energy / true_energy``)
     """
 
-    cdf = np.cumsum(energy_dispersion, axis=1)
+    bin_width = np.diff(migration_bins)
+    cdf = np.cumsum(energy_dispersion * bin_width[np.newaxis, :, np.newaxis], axis=1)
 
     n_energy_bins, _, n_fov_bins = energy_dispersion.shape
 
-    bias = np.zeros((n_energy_bins, n_fov_bins))
-    resolution = np.zeros((n_energy_bins, n_fov_bins))
+    bias = np.full((n_energy_bins, n_fov_bins), np.nan)
+    resolution = np.full((n_energy_bins, n_fov_bins), np.nan)
 
     for energy_bin in range(n_energy_bins):
         for fov_bin in range(n_fov_bins):
+            if np.count_nonzero(cdf[energy_bin, :, fov_bin]) == 0:
+                continue
 
             low, median, high = np.interp(
                 [NORM_LOWER_SIGMA, MEDIAN, NORM_UPPER_SIGMA],

pyirf/benchmarks/tests/test_bias_resolution.py

@@ -71,7 +71,6 @@ def test_energy_bias_resolution():
     assert u.isclose(resolution[1], TRUE_RES_2, rtol=0.05)
 
 
-
 def test_energy_bias_resolution():
     from pyirf.benchmarks import energy_bias_resolution_from_energy_dispersion
     from pyirf.binning import bin_center
@@ -105,6 +104,8 @@ def test_energy_bias_resolution():
         cdf[energy_bin, :, fov_bin] = norm.cdf(reco_energy, mu, sigma)
 
     edisp = cdf[:, 1:, :] - cdf[:, :-1, :]
+    bin_width = np.diff(migra_bins)
+    edisp /= bin_width[np.newaxis, :, np.newaxis]
 
     bias, resolution = energy_bias_resolution_from_energy_dispersion(
         edisp,
@@ -113,3 +114,16 @@ def test_energy_bias_resolution():
 
     assert np.allclose(bias, true_bias, atol=0.01)
     assert np.allclose(resolution, true_resolution, atol=0.01)
+
+    with_empty = np.zeros((n_energy_bins + 1, n_migra_bins, n_fov_bins))
+    with_empty[1:, :, :] = edisp
+
+    bias, resolution = energy_bias_resolution_from_energy_dispersion(
+        with_empty,
+        migra_bins,
+    )
+
+    assert np.all(np.isnan(bias[0]))
+    assert np.all(np.isnan(resolution[0]))
+    assert np.allclose(bias[1:], true_bias, atol=0.01)
+    assert np.allclose(resolution[1:], true_resolution, atol=0.01)