Repository: gammapy/gammapy
Pull request #5397 (2024-07-25 16:24:33)

Backport PR #5385 on branch v1.0.x (Fix reco_exposure computation with mask safe in ExcessMapEstimator)
Backport PR #5385 on branch v1.0.x : 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 de61acaf3..8ff27e4e3 100644
--- a/gammapy/estimators/map/excess.py
+++ b/gammapy/estimators/map/excess.py
@@ -186,6 +186,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:
@@ -194,10 +205,10 @@ 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_excess_map(self, dataset):
+    def estimate_excess_map(self, dataset, reco_exposure):
         """Estimate excess and ts maps for single dataset.
 
         If exposure is defined, a flux map is also computed.
@@ -234,9 +245,6 @@ class ExcessMapEstimator(Estimator):
         maps["sqrt_ts"] = Map.from_geom(geom, data=counts_stat.sqrt_ts)
 
         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.array) / mask.convolve(
                     kernel.array


Test patch:
diff --git a/gammapy/estimators/map/tests/test_excess.py b/gammapy/estimators/map/tests/test_excess.py
index a704a2cd8..cc8208b13 100644
--- a/gammapy/estimators/map/tests/test_excess.py
+++ b/gammapy/estimators/map/tests/test_excess.py
@@ -41,9 +41,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
 
 
@@ -155,20 +158,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):