Repository: gammapy/gammapy

Issue:
Load global RAD_MAX header keyword into RAD_MAX hdu
* The PR #3515 added support for the `RAD_MAX_2D` IRF, stored as the RAD_MAX table in GADF.

* PR #3684  uses this to for energy dependent on region calculation relying on this RAD_MAX


It should be possible to use the mechanisms introduced in #3684 also when using a fixed RAD_MAX that is stored in the `RAD_MAX` header keyword in GADF.

Possibly the best way would be to read this value if present into a `RAD_MAX_2D` object with a single energy / offset bin covering the same min/max values as the corresponding IRF.


Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/gammapy/gammapy/pull/3731?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) Report
> Merging [#3731](https://codecov.io/gh/gammapy/gammapy/pull/3731?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) (c2d5913) into [master](https://codecov.io/gh/gammapy/gammapy/commit/a46ee2df7cbc844e8f0aa28a0bb8877f8e62718d?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) (a46ee2d) will **decrease** coverage by `0.00%`.
> The diff coverage is `93.75%`.

[![Impacted file tree graph](https://codecov.io/gh/gammapy/gammapy/pull/3731/graphs/tree.svg?width=650&height=150&src=pr&token=OUjhqZUZzq&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy)](https://codecov.io/gh/gammapy/gammapy/pull/3731?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy)

```diff
@@            Coverage Diff             @@
##           master    #3731      +/-   ##
==========================================
- Coverage   93.77%   93.77%   -0.01%     
==========================================
  Files         162      162              
  Lines       19566    19582      +16     
==========================================
+ Hits        18348    18363      +15     
- Misses       1218     1219       +1     
```


| [Impacted Files](https://codecov.io/gh/gammapy/gammapy/pull/3731?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) | Coverage Δ | |
|---|---|---|
| [gammapy/irf/io.py](https://codecov.io/gh/gammapy/gammapy/pull/3731/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9pcmYvaW8ucHk=) | `96.36% <93.75%> (-1.08%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/gammapy/gammapy/pull/3731?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/gammapy/gammapy/pull/3731?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy). Last update [a46ee2d...c2d5913](https://codecov.io/gh/gammapy/gammapy/pull/3731?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).

maxnoe: @adonath To know whether we had any point-like IRFs in the datastore, we needed to actually load at least one IRF in the `DataStore.obs` method, this partly removing the lazy loading.

Do you know if there is a better solution for this? Neither obs index nor hdu index contain the information about point-like vs. full-enclosure.
maxnoe: @adonath @registerrier We are not really sure this is the correct approach. Another one would simply be adding a `rad_max` attribute or property to the IRF classes that can either be a `RadMax2D` or a simple value.
maxnoe: @adonath I implemented the property, seems to work nicely.
maxnoe: I updated the PR with more tests to make the coverage pass (it was complaining justly that I didn't have tests for the error conditions)
maxnoe: @adonath @registerrier Could you review?
bkhelifi: @maxnoe Today, there is no metadata into the obs_table or in the hdu_table specifying which type of IRFs is used for a production (or an observation)? If not, it might be a subject of thinking... as it facilities a lot the 'findability' of data (the tricks that Gammapy does is great, not not ideal in reality)
maxnoe: > @maxnoe Today, there is no metadata into the obs_table or in the hdu_table specifying which type of IRFs is used for a production (or an observation)? I

No, the information of whether IRFs are point-like or not are not stored in the HDU index if only a global RAD_MAX was used.
In the newest version, there is the option to add the `RAD_MAX_2D` HDU to the HDU index.

Since there is no `HDU` in case of a global `RAD_MAX`, I don't think the HDU index table is the correct place. Also, the information is there, it only needs to be read. 

So I don't think we need a change in the standard here. In general, you don't need `RAD_MAX` before actually reading the IRF, so reading it from the IRF is fine.
maxnoe: > He does not want to have any type of production on this machine... To be discussed else where.

In my mind, this will work like this:
* A user goes to the observatory website where they can download data
* They select what source / dates / event types / kind of analysis they want
* The data organization system creates the corresponding HDU and OBS index on the fly and prepares a download for the user

I think using a giant obs index or hdu index as the index for a whole archive is impractical. That sort of information should be stored in a proper database, from which the queries can be answered and the index files can be created.

Patch:
diff --git a/gammapy/data/observations.py b/gammapy/data/observations.py
index cfa126665..add34e57c 100644
--- a/gammapy/data/observations.py
+++ b/gammapy/data/observations.py
@@ -35,9 +35,9 @@ class Observation:
         Point spread function
     bkg : `~gammapy.irf.Background3D`
         Background rate model
-    rad_max: `~gammapy.irf.RadMax2D` or `~astropy.units.Quantity`
+    rad_max: `~gammapy.irf.RadMax2D`
         Only for point-like IRFs: RAD_MAX table (energy dependent RAD_MAX)
-        or a single angle (global RAD_MAX)
+        For a fixed RAD_MAX, create a RadMax2D with a single bin.
     gti : `~gammapy.data.GTI`
         Table with GTI start and stop time
     events : `~gammapy.data.EventList`
@@ -50,7 +50,7 @@ class Observation:
     edisp = LazyFitsData(cache=False)
     psf = LazyFitsData(cache=False)
     bkg = LazyFitsData(cache=False)
-    rad_max = LazyFitsData(cache=False)
+    _rad_max = LazyFitsData(cache=False)
     _events = LazyFitsData(cache=False)
     _gti = LazyFitsData(cache=False)
 
@@ -73,11 +73,32 @@ class Observation:
         self.edisp = edisp
         self.psf = psf
         self.bkg = bkg
-        self.rad_max = rad_max
+        self._rad_max = rad_max
         self._gti = gti
         self._events = events
         self.obs_filter = obs_filter or ObservationFilter()
 
+
+    @property
+    def rad_max(self):
+        # prevent circular import
+        from gammapy.irf import RadMax2D
+
+        if self._rad_max is not None:
+            return self._rad_max
+
+        # load once to avoid trigger lazy loading it three times
+        aeff = self.aeff
+        if aeff is not None and aeff.is_pointlike:
+            self._rad_max = RadMax2D.from_irf(aeff)
+            return self._rad_max
+
+        edisp = self.edisp
+        if edisp is not None and edisp.is_pointlike:
+            self._rad_max = RadMax2D.from_irf(self.edisp)
+
+        return self._rad_max
+
     @property
     def available_irfs(self):
         """Which irfs are available"""
diff --git a/gammapy/irf/io.py b/gammapy/irf/io.py
index 7c3817be4..8396de9d0 100644
--- a/gammapy/irf/io.py
+++ b/gammapy/irf/io.py
@@ -205,18 +205,24 @@ def load_irf_dict_from_file(filename):
         dictionary with instances of the Gammapy objects corresponding
         to the IRF components
     """
+    from .rad_max import RadMax2D
+
     filename = make_path(filename)
 
     hdulist = fits.open(make_path(filename))
 
     irf_dict = {}
 
+    is_pointlike = False
+
     for hdu in hdulist:
         hdu_class = hdu.header.get("HDUCLAS1", "").lower()
 
         if hdu_class == "response":
             hdu_class = hdu.header.get("HDUCLAS4", "").lower()
 
+            is_pointlike |= hdu.header["HDUCLAS3"] == "POINT-LIKE"
+
             loc = HDULocation(
                 hdu_class=hdu_class,
                 hdu_name=hdu.name,
@@ -237,4 +243,8 @@ def load_irf_dict_from_file(filename):
                     irf_dict[name] = data
         else:  # not an IRF component
             continue
+
+    if is_pointlike and "rad_max" not in irf_dict:
+        irf_dict["rad_max"] = RadMax2D.from_irf(irf_dict['aeff'])
+
     return irf_dict
diff --git a/gammapy/irf/rad_max.py b/gammapy/irf/rad_max.py
index 2d142fc2e..e5db98990 100644
--- a/gammapy/irf/rad_max.py
+++ b/gammapy/irf/rad_max.py
@@ -1,3 +1,5 @@
+import astropy.units as u
+
 from .core import IRF
 
 __all__ = [
@@ -28,3 +30,48 @@ class RadMax2D(IRF):
 
     tag = "rad_max_2d"
     required_axes = ["energy", "offset"]
+
+
+    @classmethod
+    def from_irf(cls, irf):
+        '''
+        Create a RadMax2D instance from another IRF component.
+
+        This reads the RAD_MAX metadata keyword from the irf and creates
+        a RadMax2D with a single bin in energy and offset using the
+        ranges from the input irf.
+
+        Parameters
+        ----------
+        irf: `~gammapy.irf.EffectiveAreaTable2D` or `~gammapy.irf.EnergyDispersion2D`
+            IRF instance from which to read the RAD_MAX and limit information
+
+        Returns
+        -------
+        rad_max: `RadMax2D`
+            `RadMax2D` object with a single bin corresponding to the fixed
+            RAD_MAX cut.
+
+        Notes
+        -----
+        This assumes the true energy axis limits are also valid for the
+        reco energy limits.
+        '''
+        if not irf.is_pointlike:
+            raise ValueError('RadMax2D.from_irf is only valid for point-like irfs')
+
+        if 'RAD_MAX' not in irf.meta:
+            raise ValueError('irf does not contain RAD_MAX keyword')
+
+        rad_max_value = irf.meta["RAD_MAX"]
+        if not isinstance(rad_max_value, float):
+            raise ValueError('RAD_MAX must be a float')
+
+        energy_axis = irf.axes["energy_true"].squash()
+        energy_axis.name = 'energy'
+        offset_axis = irf.axes["offset"].squash()
+
+        return cls(
+            data=u.Quantity([[rad_max_value]], u.deg),
+            axes=[energy_axis, offset_axis],
+        )
diff --git a/gammapy/makers/background/fov.py b/gammapy/makers/background/fov.py
index d7ffde894..15630b23a 100644
--- a/gammapy/makers/background/fov.py
+++ b/gammapy/makers/background/fov.py
@@ -20,8 +20,8 @@ class FoVBackgroundMaker(Maker):
 
     The normalization is performed outside the exclusion mask that is passed on init.
 
-    If a SkyModel is set on the input dataset its parameters
-    are frozen during the fov re-normalization.
+    If a SkyModel is set on the input dataset and method is 'fit', it' parameters
+    are frozen during the fov normalization fit.
 
     If the requirement (greater than) of either min_counts or min_npred_background is not satisfied,
     the background will not be normalised
@@ -36,7 +36,7 @@ class FoVBackgroundMaker(Maker):
         Reference norm spectral model to use for the `FoVBackgroundModel`, if none is defined
         on the dataset. By default, use pl-norm.
     min_counts : int
-        Minimum number of counts, or residuals counts if a SkyModel is set, required outside the exclusion region
+        Minimum number of counts required outside the exclusion region
     min_npred_background : float
        Minimum number of predicted background counts required outside the exclusion region
     """
@@ -132,50 +132,24 @@ class FoVBackgroundMaker(Maker):
             mask = Map.from_geom(geom=geom, data=1, dtype=bool)
         return mask
 
-    def _make_masked_summed_counts(self, dataset):
-        """"Compute the sums of the counts, npred, and bacground maps within the mask"""
-
-        npred = dataset.npred()
-        mask = dataset.mask & ~np.isnan(npred)
-        count_tot = dataset.counts.data[mask].sum()
-        npred_tot = npred.data[mask].sum()
-        bkg_tot = dataset.npred_background().data[mask].sum()
-        return {
-            "counts": count_tot,
-            "npred": npred_tot,
-            "bkg": bkg_tot,
-        }
-
     def _verify_requirements(self, dataset):
         """"Verify that the requirements of min_counts
         and min_npred_background are satisfied"""
 
-        total = self._make_masked_summed_counts(dataset)
-        not_bkg_tot = total["npred"] - total["bkg"]
+        mask = dataset.mask
+        count_tot = dataset.counts.data[mask].sum()
+        bkg_tot = dataset.npred_background().data[mask].sum()
 
-        value = (total["counts"] - not_bkg_tot) / total["bkg"]
-        if not np.isfinite(value):
+        if count_tot <= self.min_counts:
             log.warning(
-                f"FoVBackgroundMaker failed. Non-finite normalisation value for {dataset.name}. "
-                "Setting mask to False."
+                f"FoVBackgroundMaker failed. Only {int(count_tot)} counts outside exclusion mask for {dataset.name}. "
+                f"Setting mask to False."
             )
             return False
-        elif total["bkg"] <= self.min_npred_background:
+        elif bkg_tot <= self.min_npred_background:
             log.warning(
-                f"FoVBackgroundMaker failed. Only {int(total['bkg'])} background counts outside exclusion mask for {dataset.name}. "
-                "Setting mask to False."
-            )
-            return False
-        elif total["counts"] <= self.min_counts:
-            log.warning(
-                f"FoVBackgroundMaker failed. Only {int(total['counts'])} counts outside exclusion mask for {dataset.name}. "
-                "Setting mask to False."
-            )
-            return False
-        elif total["counts"] - not_bkg_tot <= 0:
-            log.warning(
-                f"FoVBackgroundMaker failed. Negative residuals counts for {dataset.name}. "
-                "Setting mask to False."
+                f"FoVBackgroundMaker failed. Only {int(bkg_tot)} background counts outside exclusion mask for {dataset.name}. "
+                f"Setting mask to False."
             )
             return False
         else:
@@ -235,7 +209,7 @@ class FoVBackgroundMaker(Maker):
             if not fit_result.success:
                 log.warning(
                     f"FoVBackgroundMaker failed. Fit did not converge for {dataset.name}. "
-                    "Setting mask to False."
+                    f"Setting mask to False."
                 )
                 dataset.mask_safe.data[...] = False
 
@@ -255,12 +229,12 @@ class FoVBackgroundMaker(Maker):
             Map dataset with scaled background model
 
         """
-        total = self._make_masked_summed_counts(dataset)
-        not_bkg_tot = total["npred"] - total["bkg"]
+        mask = dataset.mask
+        count_tot = dataset.counts.data[mask].sum()
+        bkg_tot = dataset.npred_background().data[mask].sum()
 
-        value = (total["counts"] - not_bkg_tot) / total["bkg"]
-        error = np.sqrt(total["counts"] - not_bkg_tot) / total["bkg"]
+        value = count_tot / bkg_tot
+        err = np.sqrt(count_tot) / bkg_tot
         dataset.models[f"{dataset.name}-bkg"].spectral_model.norm.value = value
-        dataset.models[f"{dataset.name}-bkg"].spectral_model.norm.error = error
-
+        dataset.models[f"{dataset.name}-bkg"].spectral_model.norm.error = err
         return dataset


Test patch:
diff --git a/gammapy/data/tests/test_data_store.py b/gammapy/data/tests/test_data_store.py
index 4de710089..6af5f4ffe 100644
--- a/gammapy/data/tests/test_data_store.py
+++ b/gammapy/data/tests/test_data_store.py
@@ -4,6 +4,7 @@ from pathlib import Path
 import pytest
 import numpy as np
 from astropy.io import fits
+import astropy.units as u
 from gammapy.data import DataStore
 from gammapy.utils.scripts import make_path
 from gammapy.utils.testing import requires_data
@@ -187,3 +188,29 @@ class TestDataStoreMaker:
         assert obs.bkg.__class__.__name__ == "Background3D"
         assert obs.edisp.__class__.__name__ == "EnergyDispersion2D"
         assert obs.psf.__class__.__name__ == "EnergyDependentMultiGaussPSF"
+
+
+@requires_data('gammapy-data')
+def test_datastore_fixed_rad_max():
+    data_store = DataStore.from_dir("$GAMMAPY_DATA/joint-crab/dl3/magic")
+    observations = data_store.get_observations([5029748], required_irf=['aeff', 'edisp'])
+
+    assert len(observations) == 1
+    obs = observations[0]
+
+    assert obs.rad_max is not None
+    assert obs.rad_max.quantity.shape == (1, 1)
+    assert u.allclose(obs.rad_max.quantity, np.sqrt(0.02) * u.deg)
+
+    # test it also works with edisp (removing aeff)
+    obs = data_store.get_observations([5029748], required_irf=['aeff', 'edisp'])[0]
+    obs.aeff = None
+    assert obs.rad_max is not None
+    assert obs.rad_max.quantity.shape == (1, 1)
+    assert u.allclose(obs.rad_max.quantity, 0.1414213 * u.deg)
+
+    # removing the last irf means we have no rad_max info
+    obs = data_store.get_observations([5029748], required_irf=['aeff', 'edisp'])[0]
+    obs.aeff = None
+    obs.edisp = None
+    assert obs.rad_max is None
diff --git a/gammapy/data/tests/test_observations.py b/gammapy/data/tests/test_observations.py
index 43009afa5..d948b5d81 100644
--- a/gammapy/data/tests/test_observations.py
+++ b/gammapy/data/tests/test_observations.py
@@ -281,6 +281,19 @@ def test_observation_read_single_file():
     assert val.unit == "m2"
 
 
+@requires_data()
+def test_observation_read_single_file_fixed_rad_max():
+    """check that for a point-like observation without the RAD_MAX_2D table
+    a RadMax2D object is generated from the RAD_MAX keyword"""
+    obs = Observation.read(
+        "$GAMMAPY_DATA/joint-crab/dl3/magic/run_05029748_DL3.fits"
+    )
+
+    assert obs.rad_max is not None
+    assert obs.rad_max.quantity.shape == (1, 1)
+    assert u.allclose(obs.rad_max.quantity, 0.1414213 * u.deg)
+
+
 @requires_data()
 class TestObservationChecker:
     def setup(self):
diff --git a/gammapy/irf/tests/test_io.py b/gammapy/irf/tests/test_io.py
index 45aa75849..4379c50b3 100644
--- a/gammapy/irf/tests/test_io.py
+++ b/gammapy/irf/tests/test_io.py
@@ -5,6 +5,7 @@ import astropy.units as u
 from astropy.io import fits
 from astropy.units import Quantity
 from gammapy.irf import (
+    RadMax2D,
     Background3D,
     EffectiveAreaTable2D,
     EnergyDispersion2D,
@@ -150,6 +151,24 @@ def test_irf_dict_from_file_duplicate_irfs(caplog, tmp_path):
     assert "loaded the PSF HDU in the dictionary" in caplog.text
 
 
+@requires_data()
+def test_irf_dict_from_file_fixed_rad_max():
+    """test that for point-like IRF without RAD_MAX_2D HDU a RadMax2D with a
+    single value is generated from the RAD_MAX header keyword"""
+    irf = load_irf_dict_from_file(
+        "$GAMMAPY_DATA/joint-crab/dl3/magic/run_05029748_DL3.fits"
+    )
+
+    assert "RAD_MAX" in irf["aeff"].meta
+    assert "rad_max" in irf
+    assert isinstance(irf["rad_max"], RadMax2D)
+
+    # check that has a single-bin in energy and offset
+    assert irf["rad_max"].axes["energy"].nbin == 1
+    assert irf["rad_max"].axes["offset"].nbin == 1
+    assert irf["rad_max"].quantity.to_value("deg") == irf["aeff"].meta["RAD_MAX"]
+
+
 class TestIRFWrite:
     def setup(self):
         self.energy_lo = np.logspace(0, 1, 10)[:-1] * u.TeV
diff --git a/gammapy/irf/tests/test_rad_max.py b/gammapy/irf/tests/test_rad_max.py
index 6071f2715..7e2c51e0c 100644
--- a/gammapy/irf/tests/test_rad_max.py
+++ b/gammapy/irf/tests/test_rad_max.py
@@ -1,6 +1,7 @@
 import numpy as np
 import astropy.units as u
 from gammapy.maps import MapAxis
+import pytest
 
 
 def test_rad_max_roundtrip(tmp_path):
@@ -31,3 +32,45 @@ def test_rad_max_roundtrip(tmp_path):
 
     assert np.all(rad_max_read.data.data == rad_max)
     assert np.all(rad_max_read.data.data == rad_max_read.data.data)
+
+
+def test_rad_max_from_irf():
+    from gammapy.irf import RadMax2D, EffectiveAreaTable2D
+
+    e_bins = 3
+    o_bins = 2
+    energy_axis = MapAxis.from_energy_bounds(1 * u.TeV, 10 * u.TeV, nbin=e_bins, name='energy_true')
+    offset_axis = MapAxis.from_bounds(0 * u.deg, 3 * u.deg, nbin=o_bins, name='offset')
+    aeff = EffectiveAreaTable2D(
+        data=u.Quantity(np.ones((e_bins, o_bins)), u.m**2, copy=False),
+        axes=[energy_axis, offset_axis],
+    )
+
+    with pytest.raises(ValueError):
+        # not a point-like IRF
+        RadMax2D.from_irf(aeff)
+
+
+    aeff.meta['is_pointlike'] = True
+
+    with pytest.raises(ValueError):
+        # missing rad_max
+        RadMax2D.from_irf(aeff)
+
+
+    aeff.meta['RAD_MAX'] = '0.2 deg'
+    with pytest.raises(ValueError):
+        # invalid format
+        RadMax2D.from_irf(aeff)
+
+    aeff.meta['RAD_MAX'] = 0.2
+    rad_max = RadMax2D.from_irf(aeff)
+
+    assert rad_max.axes['energy'].nbin == 1
+    assert rad_max.axes['offset'].nbin ==  1
+    assert rad_max.axes['energy'].edges[0] == aeff.axes['energy_true'].edges[0]
+    assert rad_max.axes['energy'].edges[1] == aeff.axes['energy_true'].edges[-1]
+    assert rad_max.axes['offset'].edges[0] == aeff.axes['offset'].edges[0]
+    assert rad_max.axes['offset'].edges[1] == aeff.axes['offset'].edges[-1]
+    assert rad_max.quantity.shape == (1, 1)
+    assert rad_max.quantity[0, 0] == 0.2 * u.deg
diff --git a/gammapy/makers/background/tests/test_fov.py b/gammapy/makers/background/tests/test_fov.py
index fd69d2410..8382e7afc 100644
--- a/gammapy/makers/background/tests/test_fov.py
+++ b/gammapy/makers/background/tests/test_fov.py
@@ -7,10 +7,9 @@ from gammapy.data import DataStore
 from gammapy.datasets import MapDataset
 from gammapy.makers import FoVBackgroundMaker, MapDatasetMaker, SafeMaskMaker
 from gammapy.maps import MapAxis, WcsGeom
-from gammapy.modeling import Fit
 from gammapy.modeling.models import (
     FoVBackgroundModel,
-    PointSpatialModel,
+    GaussianSpatialModel,
     PowerLawNormSpectralModel,
     PowerLawSpectralModel,
     SkyModel,
@@ -100,6 +99,7 @@ def test_fov_bkg_maker_scale_nocounts(obs_dataset, exclusion_mask, caplog):
     message1 = "FoVBackgroundMaker failed. Only 0 counts outside exclusion mask for test-fov. Setting mask to False."
     assert message1 in [_.message for _ in caplog.records]
 
+
 @requires_data()
 @requires_dependency("iminuit")
 def test_fov_bkg_maker_fit(obs_dataset, exclusion_mask):
@@ -161,16 +161,15 @@ def test_fov_bkg_maker_fit_nocounts(obs_dataset, exclusion_mask, caplog):
 
 @requires_data()
 @requires_dependency("iminuit")
-def test_fov_bkg_maker_with_source_model(obs_dataset, exclusion_mask, caplog):
+def test_fov_bkg_maker_fit_with_source_model(obs_dataset, exclusion_mask):
+    fov_bkg_maker = FoVBackgroundMaker(method="fit", exclusion_mask=exclusion_mask)
 
     test_dataset = obs_dataset.copy(name="test-fov")
-
-    # crab model
-    spatial_model = PointSpatialModel(
-        lon_0="83.619deg", lat_0="22.024deg", frame="icrs"
+    spatial_model = GaussianSpatialModel(
+        lon_0="0.2 deg", lat_0="0.1 deg", sigma="0.2 deg", frame="galactic"
     )
     spectral_model = PowerLawSpectralModel(
-        index=2.6, amplitude="4.5906e-11 cm-2 s-1 TeV-1", reference="1 TeV"
+        index=3, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
     )
     model = SkyModel(
         spatial_model=spatial_model, spectral_model=spectral_model, name="test-source"
@@ -179,61 +178,17 @@ def test_fov_bkg_maker_with_source_model(obs_dataset, exclusion_mask, caplog):
     bkg_model = FoVBackgroundModel(dataset_name="test-fov")
     test_dataset.models = [model, bkg_model]
 
-    # pre-fit both source and background to get reference model
-    Fit().run(test_dataset)
-    bkg_model_spec = test_dataset.models[f"{test_dataset.name}-bkg"].spectral_model
-    norm_ref = 0.897
-    assert not bkg_model_spec.norm.frozen
-    assert_allclose(bkg_model_spec.norm.value, norm_ref, rtol=1e-4)
-    assert_allclose(bkg_model_spec.tilt.value, 0.0, rtol=1e-4)
-
-    # apply scale method with pre-fitted source model and no exclusion_mask
-    bkg_model_spec.norm.value = 1
-    fov_bkg_maker = FoVBackgroundMaker(method="scale", exclusion_mask=None)
-    dataset = fov_bkg_maker.run(test_dataset)
-
-    bkg_model_spec = test_dataset.models[f"{dataset.name}-bkg"].spectral_model
-    assert_allclose(bkg_model_spec.norm.value, norm_ref, rtol=1e-4)
-    assert_allclose(bkg_model_spec.tilt.value, 0.0, rtol=1e-4)
-
-    # apply fit method with pre-fitted source model and no exlusion mask
-    bkg_model_spec.norm.value = 1
-    fov_bkg_maker = FoVBackgroundMaker(method="fit", exclusion_mask=None)
     dataset = fov_bkg_maker.run(test_dataset)
 
-    bkg_model_spec = test_dataset.models[f"{dataset.name}-bkg"].spectral_model
-    assert_allclose(bkg_model_spec.norm.value, norm_ref, rtol=1e-4)
-    assert_allclose(bkg_model_spec.tilt.value, 0.0, rtol=1e-4)
-
-    # apply scale method with pre-fitted source model and exclusion_mask
-    bkg_model_spec.norm.value = 1
-    fov_bkg_maker = FoVBackgroundMaker(method="scale", exclusion_mask=exclusion_mask)
-    dataset = fov_bkg_maker.run(test_dataset)
-
-    bkg_model_spec = test_dataset.models[f"{dataset.name}-bkg"].spectral_model
-    assert_allclose(bkg_model_spec.norm.value, 0.830779, rtol=1e-4)
-    assert_allclose(bkg_model_spec.tilt.value, 0.0, rtol=1e-4)
-
-    # apply fit method with pre-fitted source model and exlusion mask
-    bkg_model_spec.norm.value = 1
-    fov_bkg_maker = FoVBackgroundMaker(method="fit", exclusion_mask=exclusion_mask)
-    dataset = fov_bkg_maker.run(test_dataset)
-
-    bkg_model_spec = test_dataset.models[f"{dataset.name}-bkg"].spectral_model
-    assert_allclose(bkg_model_spec.norm.value, 0.830779, rtol=1e-4)
-    assert_allclose(bkg_model_spec.tilt.value, 0.0, rtol=1e-4)
-
     # Here we check that source parameters are correctly thawed after fit.
     assert not dataset.models.parameters["index"].frozen
     assert not dataset.models.parameters["lon_0"].frozen
 
-    #test  
-    model.spectral_model.amplitude.value *= 1e5
-    fov_bkg_maker = FoVBackgroundMaker(method="scale")
-    dataset = fov_bkg_maker.run(test_dataset)
-    assert "WARNING" in [_.levelname for _ in caplog.records]
-    message1 = 'FoVBackgroundMaker failed. Negative residuals counts for test-fov. Setting mask to False.'
-    assert message1 in [_.message for _ in caplog.records]
+    model = dataset.models[f"{dataset.name}-bkg"].spectral_model
+    assert not model.norm.frozen
+    assert_allclose(model.norm.value, 0.830789, rtol=1e-4)
+    assert_allclose(model.tilt.value, 0.0, rtol=1e-4)
+
 
 @requires_data()
 @requires_dependency("iminuit")
@@ -266,7 +221,7 @@ def test_fov_bkg_maker_fit_fail(obs_dataset, exclusion_mask, caplog):
     model = dataset.models[f"{dataset.name}-bkg"].spectral_model
     assert_allclose(model.norm.value, 1, rtol=1e-4)
     assert "WARNING" in [_.levelname for _ in caplog.records]
-    message1 = "FoVBackgroundMaker failed. Non-finite normalisation value for test-fov. Setting mask to False."
+    message1 = "FoVBackgroundMaker failed. Only 0 background counts outside exclusion mask for test-fov. Setting mask to False."
     assert message1 in [_.message for _ in caplog.records]