Implement additional methods for SafeMaskMaker

gammapy/cube/make.py

@@ -282,7 +282,7 @@ class SafeMaskMaker:
 
     Parameters
     ----------
-    methods : {"aeff-default", "aeff-max", "edisp-bias", "offset-max"}
+    methods : {"aeff-default", "aeff-max", "edisp-bias", "offset-max", "bkg-peak"}
         Method to use for the safe energy range. Can be a
         list with a combination of those. Resulting masks
         are combined with logical `and`. "aeff-default"
@@ -294,17 +294,27 @@ class SafeMaskMaker:
     bias_percent : float
         Percentage of the energy bias to be used as lower
         energy threshold for method "edisp-bias"
+    position : `~astropy.coordinates.SkyCoord`
+        Position at which the `aeff_percent` or `bias_percent` are computed. By default,
+        it uses the position of the center of the map.
     offset_max : str or `~astropy.units.Quantity`
         Maximum offset cut.
     """
 
-    available_methods = {"aeff-default", "aeff-max", "edisp-bias", "offset-max"}
+    available_methods = {
+        "aeff-default",
+        "aeff-max",
+        "edisp-bias",
+        "offset-max",
+        "bkg-peak",
+    }
 
     def __init__(
         self,
         methods=("aeff-default",),
         aeff_percent=10,
         bias_percent=10,
+        position=None,
         offset_max="3 deg",
     ):
         methods = set(methods)
@@ -316,6 +326,7 @@ def __init__(
         self.methods = methods
         self.aeff_percent = aeff_percent
         self.bias_percent = bias_percent
+        self.position = position
         self.offset_max = Angle(offset_max)
 
     def make_mask_offset_max(self, dataset, observation):
@@ -395,23 +406,60 @@ def make_mask_energy_edisp_bias(self, dataset):
 
         Parameters
         ----------
-        dataset : `SpectrumDataset` or `SpectrumDatasetOnOff`
+        dataset : `MapDataset`, `MapDatasetOnOff`, `SpectrumDataset` or `SpectrumDatasetOnOff`
             Dataset to compute mask for.
 
         Returns
         -------
         mask_safe : `~numpy.ndarray`
             Safe data range mask.
         """
+        edisp = dataset.edisp
+
         if isinstance(dataset, (MapDataset, MapDatasetOnOff)):
-            raise NotImplementedError(
-                "'edisp-bias' method currently only supported for spectral datasets"
-            )
+            position = self.position
+            if position is None:
+                position = dataset.counts.geom.center_skydir
+            e_reco = dataset.counts.geom.get_axis_by_name("energy").edges
+            edisp = edisp.get_energy_dispersion(position, e_reco)
+            counts = dataset.counts.geom
+        else:
+            counts = dataset.counts
 
-        e_min = dataset.edisp.get_bias_energy(self.bias_percent / 100)
-        return dataset.counts.energy_mask(emin=e_min)
+        e_min = edisp.get_bias_energy(self.bias_percent / 100)
+        return counts.energy_mask(emin=e_min)
+
+    @staticmethod
+    def make_mask_energy_bkg_peak(dataset):
+        """Make safe energy mask based on the binned background.
+
+        The energy threshold is defined as the upper edge of the energy
+        bin with the highest predicted background rate. This method is motivated
+        by its use in the HESS DL3 validation paper: https://arxiv.org/pdf/1910.08088.pdf
+
+        Parameters
+        ----------
+        dataset : `MapDataset`, `MapDatasetOnOff`, `SpectrumDataset` or `SpectrumDatasetOnOff`
+            Dataset to compute mask for.
+
+        Returns
+        -------
+        mask_safe : `~numpy.ndarray`
+            Safe data range mask.
+        """
+
+        if isinstance(dataset, (MapDataset, MapDatasetOnOff)):
+            background_spectrum = dataset.background_model.map.get_spectrum()
+            counts = dataset.counts.geom
+        else:
+            background_spectrum = dataset.background
+            counts = dataset.counts
+
+        idx = np.argmax(background_spectrum.data)
+        e_min = background_spectrum.energy.edges[idx + 1]
+        return counts.energy_mask(emin=e_min)
 
-    def run(self, dataset, observation):
+    def run(self, dataset, observation=None):
         """Make safe data range mask.
 
         Parameters
@@ -440,6 +488,9 @@ def run(self, dataset, observation):
         if "edisp-bias" in self.methods:
             mask_safe &= self.make_mask_energy_edisp_bias(dataset)
 
+        if "bkg-peak" in self.methods:
+            mask_safe &= self.make_mask_energy_bkg_peak(dataset)
+
         if isinstance(dataset, (MapDataset, MapDatasetOnOff)):
             mask_safe = Map.from_geom(dataset.counts.geom, data=mask_safe)
 

gammapy/cube/tests/test_make.py

@@ -203,12 +203,16 @@ def test_map_maker_obs(observations):
 def test_safe_mask_maker(observations):
     obs = observations[0]
 
-    axis = MapAxis.from_edges([0.1, 1, 10], name="energy", interp="log", unit="TeV")
+    axis = MapAxis.from_bounds(
+        0.1, 10, nbin=16, unit="TeV", name="energy", interp="log"
+    )
     geom = WcsGeom.create(npix=(11, 11), axes=[axis], skydir=obs.pointing_radec)
 
     empty_dataset = MapDataset.create(geom=geom)
     dataset_maker = MapDatasetMaker(offset_max="3 deg")
-    safe_mask_maker = SafeMaskMaker(offset_max="3 deg")
+    safe_mask_maker = SafeMaskMaker(
+        offset_max="3 deg", bias_percent=0.02, position=obs.pointing_radec
+    )
     dataset = dataset_maker.run(empty_dataset, obs)
 
     mask_offset = safe_mask_maker.make_mask_offset_max(dataset=dataset, observation=obs)
@@ -217,14 +221,15 @@ def test_safe_mask_maker(observations):
     mask_energy_aeff_default = safe_mask_maker.make_mask_energy_aeff_default(
         dataset=dataset, observation=obs
     )
-    assert_allclose(mask_energy_aeff_default.sum(), 242)
+    assert_allclose(mask_energy_aeff_default.sum(), 1936)
 
-    with pytest.raises(NotImplementedError) as excinfo:
-        safe_mask_maker.make_mask_energy_edisp_bias(dataset)
+    mask_edisp_bias = safe_mask_maker.make_mask_energy_edisp_bias(dataset)
+    assert_allclose(mask_edisp_bias.sum(), 1815)
 
-    assert "only supported" in str(excinfo.value)
+    mask_bkg_peak = safe_mask_maker.make_mask_energy_bkg_peak(dataset)
+    assert_allclose(mask_bkg_peak.sum(), 1815)
 
     with pytest.raises(NotImplementedError) as excinfo:
-        safe_mask_maker.make_mask_energy_edisp_bias(dataset)
+        safe_mask_maker.make_mask_energy_aeff_max(dataset)
 
     assert "only supported" in str(excinfo.value)

gammapy/spectrum/tests/test_make.py

@@ -128,6 +128,9 @@ def test_safe_mask_maker_dl3(spectrum_dataset_maker_crab, observations_hess_dl3)
     mask_safe = safe_mask_maker.make_mask_energy_edisp_bias(dataset)
     assert mask_safe.sum() == 3
 
+    mask_safe = safe_mask_maker.make_mask_energy_bkg_peak(dataset)
+    assert mask_safe.sum() == 3
+
 
 @requires_data()
 def test_safe_mask_maker_dc1(spectrum_dataset_maker_gc, observations_cta_dc1):