Repository: gammapy/gammapy
Pull request #5385 (2024-07-21 14:27:09)

Fix reco_exposure computation with mask safe in ExcessMapEstimator
Fix reco_exposure computation with mask safe in ExcessMapEstimator


Patch:
diff --git a/gammapy/estimators/map/excess.py b/gammapy/estimators/map/excess.py
index c50d6e68a..1c2b8cbe4 100644
--- a/gammapy/estimators/map/excess.py
+++ b/gammapy/estimators/map/excess.py
@@ -228,6 +228,17 @@ class ExcessMapEstimator(Estimator):
         resampled_dataset = dataset.resample_energy_axis(
             energy_axis=axis, name=dataset.name
         )
+
+        if dataset.exposure:
+            reco_exposure = estimate_exposure_reco_energy(
+                dataset, self.spectral_model, normalize=False
+            )
+            reco_exposure = reco_exposure.resample_axis(
+                axis=axis, weights=dataset.mask_safe
+            )
+        else:
+            reco_exposure = None
+
         if isinstance(dataset, MapDatasetOnOff):
             resampled_dataset.models = dataset.models
         else:
@@ -236,7 +247,7 @@ class ExcessMapEstimator(Estimator):
             )
             resampled_dataset.models = None
 
-        result = self.estimate_excess_map(resampled_dataset)
+        result = self.estimate_excess_map(resampled_dataset, reco_exposure)
         return result
 
     def estimate_kernel(self, dataset):
@@ -283,7 +294,7 @@ class ExcessMapEstimator(Estimator):
             mask = Map.from_geom(dataset.counts.geom, data=True, dtype=bool)
         return mask
 
-    def estimate_exposure_reco_energy(self, dataset, kernel, mask):
+    def estimate_exposure_reco_energy(self, dataset, kernel, mask, reco_exposure):
         """Estimate exposure map in reconstructed energy for a single dataset
            assuming the given spectral_model shape.
 
@@ -302,18 +313,16 @@ class ExcessMapEstimator(Estimator):
             Reconstructed exposure map.
         """
         if dataset.exposure:
-            reco_exposure = estimate_exposure_reco_energy(
-                dataset, self.spectral_model, normalize=False
-            )
             with np.errstate(invalid="ignore", divide="ignore"):
                 reco_exposure = reco_exposure.convolve(kernel.data) / mask.convolve(
                     kernel.data
                 )
         else:
             reco_exposure = 1
+
         return reco_exposure
 
-    def estimate_excess_map(self, dataset):
+    def estimate_excess_map(self, dataset, reco_exposure):
         """Estimate excess and test statistic maps for a single dataset.
 
         If exposure is defined, a flux map is also computed.
@@ -342,7 +351,9 @@ class ExcessMapEstimator(Estimator):
         maps["ts"] = Map.from_geom(geom, data=counts_stat.ts)
         maps["sqrt_ts"] = Map.from_geom(geom, data=counts_stat.sqrt_ts)
 
-        reco_exposure = self.estimate_exposure_reco_energy(dataset, kernel, mask)
+        reco_exposure = self.estimate_exposure_reco_energy(
+            dataset, kernel, mask, reco_exposure
+        )
 
         with np.errstate(invalid="ignore", divide="ignore"):
             maps["norm"] = maps["npred_excess"] / reco_exposure


Test patch:
diff --git a/gammapy/estimators/map/tests/test_excess.py b/gammapy/estimators/map/tests/test_excess.py
index abb190e07..d03e03a27 100644
--- a/gammapy/estimators/map/tests/test_excess.py
+++ b/gammapy/estimators/map/tests/test_excess.py
@@ -44,9 +44,12 @@ def simple_dataset_mask_safe():
     geom = WcsGeom.create(npix=20, binsz=0.02, axes=[axis])
     dataset = MapDataset.create(geom)
     dataset.mask_safe += np.ones(dataset.data_shape, dtype=bool)
-    dataset.mask_safe.data[0, :, :] = False
+    dataset.mask_safe.data[0, :10, :] = False
     dataset.counts += 2
     dataset.background += 1
+    dataset.exposure.data += 1
+    dataset.exposure.data[0, :, :] = 2
+    dataset.exposure.data[2, :10, :] = 4
     return dataset
 
 
@@ -170,20 +173,44 @@ def test_significance_map_estimator_map_dataset(simple_dataset):
 
 
 def test_significance_map_estimator_map_dataset_mask_safe(simple_dataset_mask_safe):
-    simple_dataset_mask_safe.exposure = None
+
     estimator = ExcessMapEstimator(0.1 * u.deg, selection_optional=["all"])
 
     result = estimator.run(simple_dataset_mask_safe)
 
-    assert_allclose(result["npred"].data[0, 10, 10], 324)
-    assert_allclose(result["npred_excess"].data[0, 10, 10], 162)
-    assert_allclose(result["npred_background"].data[0, 10, 10], 162)
-    assert_allclose(result["sqrt_ts"].data[0, 10, 10], 11.187468, atol=1e-5)
+    assert_allclose(result["npred"].data[0, 10, 10], 416)
+    assert_allclose(result["npred_excess"].data[0, 10, 10], 208)
+    assert_allclose(result["npred_background"].data[0, 10, 10], 208)
+    assert_allclose(result["sqrt_ts"].data[0, 10, 10], 12.676681, atol=1e-5)
+
+    assert_allclose(result["npred_excess_err"].data[0, 10, 10], 20.396078, atol=1e-3)
+    assert_allclose(result["npred_excess_errp"].data[0, 10, 10], 20.730114, atol=1e-3)
+    assert_allclose(result["npred_excess_errn"].data[0, 10, 10], 20.063506, atol=1e-3)
+    assert_allclose(result["npred_excess_ul"].data[0, 10, 10], 250.136257, atol=1e-3)
+
+    assert_allclose(result["flux"].data[0, 0, 0], 0.016359, atol=1e-3)
+    assert_allclose(result["flux"].data[0, 10, 10], 0.018004, atol=1e-3)
+
+    reco_exposure = result["npred_excess"] / result["norm"]
+
+    assert_allclose(np.unique(reco_exposure).min(), 3.1469110e-08, rtol=1e-5)
+    assert_allclose(np.unique(reco_exposure).max(), 1.7745566e-07, rtol=1e-5)
+
+    simple_dataset_mask_safe.exposure = None
+
+    result = estimator.run(simple_dataset_mask_safe)
+
+    assert_allclose(result["npred_excess"].data[0, 10, 10], 208)
+    assert_allclose(result["npred_background"].data[0, 10, 10], 208)
+    assert_allclose(result["sqrt_ts"].data[0, 10, 10], 12.676681, atol=1e-5)
+
+    assert_allclose(result["npred_excess_err"].data[0, 10, 10], 20.396078, atol=1e-3)
+    assert_allclose(result["npred_excess_errp"].data[0, 10, 10], 20.730114, atol=1e-3)
+    assert_allclose(result["npred_excess_errn"].data[0, 10, 10], 20.063506, atol=1e-3)
+    assert_allclose(result["npred_excess_ul"].data[0, 10, 10], 250.136257, atol=1e-3)
 
-    assert_allclose(result["npred_excess_err"].data[0, 10, 10], 18.0, atol=1e-3)
-    assert_allclose(result["npred_excess_errp"].data[0, 10, 10], 18.334, atol=1e-3)
-    assert_allclose(result["npred_excess_errn"].data[0, 10, 10], 17.668, atol=1e-3)
-    assert_allclose(result["npred_excess_ul"].data[0, 10, 10], 199.345, atol=1e-3)
+    assert_allclose(result["flux"].data[0, 0, 0], 5.148e-10, atol=1e-3)
+    assert_allclose(result["flux"].data[0, 10, 10], 2.0592e-09, atol=1e-3)
 
 
 def test_significance_map_estimator_map_dataset_exposure(simple_dataset):