Repository: gammapy/gammapy
Pull request #5285 (2024-06-03 13:32:07)

Add support for stat_scan in TSMapEstimator
Add support for stat_scan in TSMapEstimator with 2 methods:
* "stat_scan_local": estimate likelihood profile around best norm
* "stat_scan": estimate likelihood profile on a fixed norm range


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/gammapy/gammapy/pull/5285?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) Report
Attention: Patch coverage is `88.23529%` with `8 lines` in your changes are missing coverage. Please review.
> Project coverage is 94.45%. Comparing base [(`f85fb7b`)](https://app.codecov.io/gh/gammapy/gammapy/commit/f85fb7b40717542f697d1707ce18d0aec211faaf?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) to head [(`05a73c2`)](https://app.codecov.io/gh/gammapy/gammapy/commit/05a73c2404783b1780a25c3afa3b5ee14b0c6efc?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).
> Report is 1 commits behind head on main.

> :exclamation: **Current head 05a73c2 differs from pull request most recent head af644ad**
> 
> Please [upload](https://docs.codecov.com/docs/codecov-uploader) reports for the commit af644ad to get more accurate results.

| [Files](https://app.codecov.io/gh/gammapy/gammapy/pull/5285?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) | Patch % | Lines |
|---|---|---|
| [gammapy/estimators/map/ts.py](https://app.codecov.io/gh/gammapy/gammapy/pull/5285?src=pr&el=tree&filepath=gammapy%2Festimators%2Fmap%2Fts.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9lc3RpbWF0b3JzL21hcC90cy5weQ==) | 88.23% | [8 Missing :warning: ](https://app.codecov.io/gh/gammapy/gammapy/pull/5285?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #5285      +/-   ##
==========================================
- Coverage   94.46%   94.45%   -0.02%     
==========================================
  Files         232      232              
  Lines       35636    35698      +62     
==========================================
+ Hits        33664    33718      +54     
- Misses       1972     1980       +8     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/gammapy/gammapy/pull/5285?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).

QRemy: > Why do you want to provide two options `stat_scan_local` and `stat_scan`? Why not a similar API as for the `FluxPointsEstimator` where we have a configurable `norm` parameter which, by default corresponds to the `stat_scan_local` case ?

The API is the same, there is  a `norm_local` and `norm` the difference it how they are applied.
The problem is to define a norm range and sampling that is acceptable for all the bin of the map so I have two approach :

- The `stat_scan_local` uses the  `norm_local` parameter  that is multiplied with the best fit norm map and has a small range and only few points because we start close to the solution. In that case the outputs have  are "stat_scan_local" and "norm_scan_values" a map that gives the different norm values for each (E, l,b) pixels.

- The `stat_scan` uses the `norm` parameter that apply to the reference_model and has a wider range of values and finer sampling because a priori we don't know how far the reference_model is from the best model. In that case the output is just "stat_scan" a map where the norm values are given by the norm axis (same for all l,b pixels). This is quite convenient because if two flux maps were computed with the same reference_model and norm parameter their "stat_scan" maps can be summed to compute a joint flux maps as done in #5289.
QRemy: > Can you please describe again the motivation for the second "local" stat scan? I read your comment above, however I have still not understood why it is needed. Stat scans can have arbitrary sampling so if we need a finer sampling close to the minimum we could still do it with a single array...

The main motivations are better performances and easier configurability.

The best-fit norm can wary largely accros the map  so if we want to use a single norm values array with a finer binning and a wide range as implemented in "stat_scan" the computation time can be very high. Moreover we have to run a first time the estimator without the stat_scan to get the optimal norm range and then update estimator norm parameter and run the stat_scan. 
 
As "stat_scan_local" is defined with respect to the best-fit norm of each (l,b) pixel even if we scan a very small range of norms and it should be well defined for all pixels.
adonath: Thanks @QRemy, that makes sense. I guess for flux points this was much less of a problem than for a flux map. However one main motivation for the TS profiles is the possibility to combine them afterwards, which requires them to be on the same grid (at least it is much simpler). Alternatively one could allow for a spatially varying reference model, but I guess this is currently not possible with `FluxMaps`.

Naively I would assume a log grid with like ~11 bins should be sufficient. Do you maybe have a distribution of norm values for a typical map? 

In any case it will be good to check, what "gttscube" does and who this problem is solved for the "TS Cube" format I linked above.
QRemy: > Thanks a lot @QRemy, this is really nice! Please also take a look at https://gamma-astro-data-formats.readthedocs.io/en/v0.3/spectra/binned_likelihoods/index.html (Section "Likelihood SED Cube") to see whether the data could be serialized in the documented format.

the required columns for the scan are:
``
dloglike_scan – ndim: 3, Dimension: nrows x nebins x nnorms
norm_scan – ndim: 3, Dimension: nrows x nebins x nnorms
```
This is similar to the output of "stat_scan_local" except that the spatial dimensions are flatten.
For "stat_scan" I return only the stat map because the norm values are simply defined by the geom.axes["norm"].center

It should be easy to convert the gammapy output to this format if the reshaping from 2d to 1d is done in the same way. 


Patch:
diff --git a/gammapy/estimators/flux.py b/gammapy/estimators/flux.py
index 76127d342..e4c7a2662 100644
--- a/gammapy/estimators/flux.py
+++ b/gammapy/estimators/flux.py
@@ -4,8 +4,8 @@ import numpy as np
 from gammapy.datasets import Datasets
 from gammapy.datasets.actors import DatasetsActor
 from gammapy.estimators.parameter import ParameterEstimator
+from gammapy.estimators.utils import _get_default_norm
 from gammapy.maps import Map, MapAxis
-from gammapy.modeling import Parameter
 from gammapy.modeling.models import ScaleSpectralModel
 from gammapy.utils.deprecation import deprecated_attribute, deprecated_renamed_argument
 
@@ -97,27 +97,15 @@ class FluxEstimator(ParameterEstimator):
     ):
 
         self.source = source
-        if norm is None or isinstance(norm, dict):
-            norm_kwargs = dict(
-                name="norm",
-                value=1,
-                unit="",
-                interp="log",
-                frozen=False,
-                scan_min=norm_min,
-                scan_max=norm_max,
-                scan_n_values=norm_n_values,
-                scan_values=norm_values,
-            )
-            if isinstance(norm, dict):
-                norm_kwargs.update(norm)
-            try:
-                norm = Parameter(**norm_kwargs)
-            except TypeError as error:
-                raise TypeError(f"Invalid dict key for norm init : {error}")
-        if norm.name != "norm":
-            raise ValueError("norm.name is not 'norm'")
-        self.norm = norm
+
+        self.norm = _get_default_norm(
+            norm,
+            scan_min=norm_min,
+            scan_max=norm_max,
+            scan_n_values=norm_n_values,
+            scan_values=norm_values,
+            interp="log",
+        )
 
         super().__init__(
             null_value=0,
diff --git a/gammapy/estimators/map/core.py b/gammapy/estimators/map/core.py
index cee19ce84..77445a643 100644
--- a/gammapy/estimators/map/core.py
+++ b/gammapy/estimators/map/core.py
@@ -78,6 +78,8 @@ VALID_QUANTITIES = [
     "stat",
     "stat_scan",
     "stat_null",
+    "stat_scan",
+    "dnde_scan_values",
     "niter",
     "is_ul",
     "counts",
@@ -93,6 +95,8 @@ OPTIONAL_QUANTITIES_COMMON = [
     "npred_excess",
     "stat",
     "stat_null",
+    "stat_scan",
+    "dnde_scan_values",
     "niter",
     "is_ul",
     "counts",
@@ -488,6 +492,12 @@ class FluxMaps:
         self._check_quantity("stat_scan")
         return self._data["stat_scan"]
 
+    @property
+    def dnde_scan_values(self):
+        """Fit statistic norm scan values."""
+        self._check_quantity("dnde_scan_values")
+        return self._data["dnde_scan_values"]
+
     @property
     def stat(self):
         """Fit statistic value."""
diff --git a/gammapy/estimators/map/ts.py b/gammapy/estimators/map/ts.py
index 900c6af97..55b602b17 100644
--- a/gammapy/estimators/map/ts.py
+++ b/gammapy/estimators/map/ts.py
@@ -4,19 +4,25 @@ import warnings
 from itertools import repeat
 import numpy as np
 import scipy.optimize
+from scipy.interpolate import InterpolatedUnivariateSpline
 from astropy.coordinates import Angle
 from astropy.utils import lazyproperty
 import gammapy.utils.parallel as parallel
 from gammapy.datasets.map import MapEvaluator
 from gammapy.datasets.utils import get_nearest_valid_exposure_position
-from gammapy.maps import Map, Maps
+from gammapy.maps import Map, MapAxis, Maps
 from gammapy.modeling.models import PointSpatialModel, PowerLawSpectralModel, SkyModel
-from gammapy.stats import cash_sum_cython, f_cash_root_cython, norm_bounds_cython
+from gammapy.stats import cash, cash_sum_cython, f_cash_root_cython, norm_bounds_cython
 from gammapy.utils.array import shape_2N, symmetric_crop_pad_width
 from gammapy.utils.pbar import progress_bar
 from gammapy.utils.roots import find_roots
 from ..core import Estimator
-from ..utils import estimate_exposure_reco_energy
+from ..utils import (
+    _generate_scan_values,
+    _get_default_norm,
+    _get_norm_scan_values,
+    estimate_exposure_reco_energy,
+)
 from .core import FluxMaps
 
 __all__ = ["TSMapEstimator"]
@@ -84,6 +90,7 @@ class TSMapEstimator(Estimator, parallel.ParallelMixin):
             * "all": all the optional steps are executed
             * "errn-errp": estimate asymmetric error on flux.
             * "ul": estimate upper limits on flux.
+            * "stat_scan": estimate likelihood profile
 
         Default is None so the optional steps are not executed.
     energy_edges : list of `~astropy.units.Quantity`, optional
@@ -94,6 +101,14 @@ class TSMapEstimator(Estimator, parallel.ParallelMixin):
     sum_over_energy_groups : bool
         Whether to sum over the energy groups or fit the norm on the full energy
         cube.
+    norm : `~gammapy.modeling.Parameter` or dict
+        Norm parameter used for the likelihood profile computation on a fixed norm range.
+        Only used for "stat_scan" in `selection_optional`.
+        Default is None and a new parameter is created automatically,
+        with value=1, name="norm", scan_min=-100, scan_max=100,
+        and values sampled such as we can probe a 0.1% relative error on the norm.
+        If a dict is given the entries should be a subset of
+        `~gammapy.modeling.Parameter` arguments.
     n_jobs : int
         Number of processes used in parallel for the computation. Default is one,
         unless `~gammapy.utils.parallel.N_JOBS_DEFAULT` was modified. The number
@@ -150,7 +165,7 @@ class TSMapEstimator(Estimator, parallel.ParallelMixin):
     """
 
     tag = "TSMapEstimator"
-    _available_selection_optional = ["errn-errp", "ul"]
+    _available_selection_optional = ["errn-errp", "ul", "stat_scan"]
 
     def __init__(
         self,
@@ -166,12 +181,15 @@ class TSMapEstimator(Estimator, parallel.ParallelMixin):
         sum_over_energy_groups=True,
         n_jobs=None,
         parallel_backend=None,
+        norm=None,
     ):
         if kernel_width is not None:
             kernel_width = Angle(kernel_width)
 
         self.kernel_width = kernel_width
 
+        self.norm = _get_default_norm(norm, scan_values=_generate_scan_values())
+
         if model is None:
             model = SkyModel(
                 spectral_model=PowerLawSpectralModel(),
@@ -197,6 +215,7 @@ class TSMapEstimator(Estimator, parallel.ParallelMixin):
             n_sigma_ul=self.n_sigma_ul,
             selection_optional=selection_optional,
             ts_threshold=threshold,
+            norm=self.norm,
         )
 
     @property
@@ -214,11 +233,19 @@ class TSMapEstimator(Estimator, parallel.ParallelMixin):
             "success",
         ]
 
-        if "errn-errp" in self.selection_optional:
-            selection += ["norm_errp", "norm_errn"]
-
-        if "ul" in self.selection_optional:
-            selection += ["norm_ul"]
+        if "stat_scan" in self.selection_optional:
+            selection += [
+                "dnde_scan_values",
+                "stat_scan",
+                "norm_errp",
+                "norm_errn",
+                "norm_ul",
+            ]
+        else:
+            if "errn-errp" in self.selection_optional:
+                selection += ["norm_errp", "norm_errn"]
+            if "ul" in self.selection_optional:
+                selection += ["norm_ul"]
 
         return selection
 
@@ -443,10 +470,34 @@ class TSMapEstimator(Estimator, parallel.ParallelMixin):
         j, i = zip(*positions)
 
         geom = maps["counts"].geom.squash(axis_name="energy")
+        energy_axis = maps["counts"].geom.axes["energy"]
+        dnde_ref = self.model.spectral_model(energy_axis.center)
 
         for name in self.selection_all:
-            m = Map.from_geom(geom=geom, data=np.nan, unit="")
-            m.data[0, j, i] = [_[name] for _ in results]
+            if name in ["dnde_scan_values", "stat_scan"]:
+                norm_bin_axis = MapAxis(
+                    range(len(results[0]["dnde_scan_values"])),
+                    interp="lin",
+                    node_type="center",
+                    name="dnde_bin",
+                )
+
+                axes = geom.axes + [norm_bin_axis]
+                geom_scan = geom.to_image().to_cube(axes)
+
+                if name == "dnde_scan_values":
+                    unit = dnde_ref.unit
+                    factor = dnde_ref.value
+                else:
+                    unit = ""
+                    factor = 1
+
+                m = Map.from_geom(geom_scan, data=np.nan, unit=unit)
+                m.data[:, 0, j, i] = np.array([_[name] for _ in results]).T * factor
+
+            else:
+                m = Map.from_geom(geom=geom, data=np.nan, unit="")
+                m.data[0, j, i] = [_[name] for _ in results]
             result[name] = m
 
         return result
@@ -602,7 +653,7 @@ class SimpleMapDataset:
 class BrentqFluxEstimator(Estimator):
     """Single parameter flux estimator."""
 
-    _available_selection_optional = ["errn-errp", "ul"]
+    _available_selection_optional = ["errn-errp", "ul", "stat_scan"]
     tag = "BrentqFluxEstimator"
 
     def __init__(
@@ -613,6 +664,7 @@ class BrentqFluxEstimator(Estimator):
         selection_optional=None,
         max_niter=20,
         ts_threshold=None,
+        norm=None,
     ):
         self.rtol = rtol
         self.n_sigma = n_sigma
@@ -620,6 +672,7 @@ class BrentqFluxEstimator(Estimator):
         self.selection_optional = selection_optional
         self.max_niter = max_niter
         self.ts_threshold = ts_threshold
+        self.norm = norm
 
     def estimate_best_fit(self, dataset):
         """Estimate best fit norm parameter.
@@ -745,6 +798,63 @@ class BrentqFluxEstimator(Estimator):
         )
         return {"norm_errn": flux_errn, "norm_errp": flux_errp}
 
+    def estimate_scan(self, dataset, result):
+        """Compute likelihood profile.
+
+        Parameters
+        ----------
+        dataset : `SimpleMapDataset`
+            Simple map dataset.
+
+        Returns
+        -------
+        result : dict
+            Result dictionary including 'stat_scan'.
+        """
+
+        sparse_norms = _get_norm_scan_values(self.norm, result)
+
+        scale = sparse_norms[None, :]
+        model = dataset.model.ravel()[:, None]
+        background = dataset.background.ravel()[:, None]
+        counts = dataset.counts.ravel()[:, None]
+        stat_scan = cash(counts, model * scale + background)
+
+        stat_scan_local = stat_scan.sum(axis=0) - result["stat_null"]
+
+        spline = InterpolatedUnivariateSpline(
+            sparse_norms, stat_scan_local, k=1, ext="raise", check_finite=True
+        )
+        stat_scan = spline(self.norm.scan_values)
+
+        ts = -stat_scan.min()
+        ind = stat_scan.argmin()
+        norm = self.norm.scan_values[ind]
+
+        maskp = self.norm.scan_values > norm
+        stat_diff = stat_scan - stat_scan.min()
+        ind = np.abs(stat_diff - self.n_sigma**2)[maskp].argmin()
+        norm_errn = norm - self.norm.scan_values[maskp][ind]
+
+        ind = np.abs(stat_diff - self.n_sigma**2)[~maskp].argmin()
+        norm_errp = self.norm.scan_values[~maskp][ind] - norm
+
+        ind = np.abs(stat_diff - self.n_sigma_ul**2)[~maskp].argmin()
+        norm_ul = self.norm.scan_values[~maskp][ind]
+
+        norm_err = (norm_errn + norm_errp) / 2
+
+        return dict(
+            ts=ts,
+            norm=norm,
+            norm_err=norm_err,
+            norm_errn=norm_errn,
+            norm_errp=norm_errp,
+            norm_ul=norm_ul,
+            stat_scan=stat_scan_local,
+            dnde_scan_values=sparse_norms,
+        )
+
     def estimate_default(self, dataset):
         """Estimate default norm.
 
@@ -796,16 +906,20 @@ class BrentqFluxEstimator(Estimator):
         else:
             result = self.estimate_best_fit(dataset)
 
-        norm = result["norm"]
-        result["npred"] = dataset.npred(norm=norm).sum()
-        result["npred_excess"] = result["npred"] - dataset.npred(norm=0).sum()
-
         if "ul" in self.selection_optional:
             result.update(self.estimate_ul(dataset, result))
 
         if "errn-errp" in self.selection_optional:
             result.update(self.estimate_errn_errp(dataset, result))
 
+        if "stat_scan" in self.selection_optional:
+            result.update(self.estimate_scan(dataset, result))
+
+        norm = result["norm"]
+        result["npred"] = dataset.npred(norm=norm).sum()
+        result["npred_excess"] = result["npred"] - dataset.npred(norm=0).sum()
+        result["stat"] = dataset.stat_sum(norm=norm)
+
         return result
 
 
diff --git a/gammapy/estimators/utils.py b/gammapy/estimators/utils.py
index bb2cd2fe6..f3df04d3b 100644
--- a/gammapy/estimators/utils.py
+++ b/gammapy/estimators/utils.py
@@ -8,6 +8,7 @@ from astropy.table import Table
 from gammapy.datasets import SpectrumDataset, SpectrumDatasetOnOff
 from gammapy.datasets.map import MapEvaluator
 from gammapy.maps import Map, MapAxis, TimeMapAxis, WcsNDMap
+from gammapy.modeling import Parameter
 from gammapy.modeling.models import (
     ConstantFluxSpatialModel,
     PowerLawSpectralModel,
@@ -846,3 +847,75 @@ def combine_flux_maps(maps, method="gaussian_errors", reference_model=None):
 
     kwargs = dict(sed_type="dnde", gti=gti, reference_model=reference_model, meta=meta)
     return FluxMaps.from_maps(dict(dnde=mean, dnde_err=sigma, ts=ts), **kwargs)
+
+
+def _generate_scan_values(power_min=-4, power_max=2, relative_error=1e-2):
+    """Values sampled such as we can probe a given `relative_error` on the norm
+    between 10**`power_min` and 10**`power_max`.
+
+    """
+    arrays = []
+    for power in range(power_min, power_max):
+        vmin = 10**power
+        vmax = 10 ** (power + 1)
+        bin_per_decade = int((vmax - vmin) / (vmin * relative_error))
+        arrays.append(np.linspace(vmin, vmax, bin_per_decade + 1, dtype=np.float32))
+    scan_1side = np.unique(np.concatenate(arrays))
+    return np.concatenate((-scan_1side[::-1], [0], scan_1side))
+
+
+def _get_default_norm(
+    norm,
+    scan_min=None,
+    scan_max=None,
+    scan_n_values=None,
+    scan_values=None,
+    interp="lin",
+):
+    """create default norm parameter"""
+    if norm is None or isinstance(norm, dict):
+        norm_kwargs = dict(
+            name="norm",
+            value=1,
+            unit="",
+            interp=interp,
+            frozen=False,
+            scan_min=scan_min,
+            scan_max=scan_max,
+            scan_n_values=scan_n_values,
+            scan_values=scan_values,
+        )
+        if isinstance(norm, dict):
+            norm_kwargs.update(norm)
+        try:
+            norm = Parameter(**norm_kwargs)
+        except TypeError as error:
+            raise TypeError(f"Invalid dict key for norm init : {error}")
+    if norm.name != "norm":
+        raise ValueError("norm.name is not 'norm'")
+    return norm
+
+
+def _get_norm_scan_values(norm, result):
+    """Compute norms based on the fit result to sample the stat profile at different scales."""
+
+    norm_err = result["norm_err"]
+    norm_value = result["norm"]
+    if ~np.isfinite(norm_err) or norm_err == 0:
+        norm_err = 0.1
+    if ~np.isfinite(norm_value) or norm_value == 0:
+        norm_value = 1.0
+    sparse_norms = np.concatenate(
+        (
+            norm_value + np.linspace(-2.5, 2.5, 51) * norm_err,
+            norm_value + np.linspace(-10, 10, 21) * norm_err,
+            np.abs(norm_value) * np.linspace(-10, 10, 21),
+            np.linspace(-10, 10, 21),
+            np.linspace(norm.scan_values[0], norm.scan_values[-1], 2),
+        )
+    )
+    sparse_norms = np.unique(sparse_norms)
+    if len(sparse_norms) != 109:
+        rand_norms = 20 * np.random.rand(109 - len(sparse_norms)) - 10
+        sparse_norms = np.concatenate((sparse_norms, rand_norms))
+    return np.sort(sparse_norms)


Test patch:
diff --git a/gammapy/estimators/map/tests/test_ts.py b/gammapy/estimators/map/tests/test_ts.py
index 40913e913..a6b0e3950 100644
--- a/gammapy/estimators/map/tests/test_ts.py
+++ b/gammapy/estimators/map/tests/test_ts.py
@@ -199,7 +199,7 @@ def test_compute_ts_map_psf(fermi_dataset):
     model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)
 
     estimator = TSMapEstimator(
-        model=model, kernel_width="1 deg", selection_optional="all"
+        model=model, kernel_width="1 deg", selection_optional=["ul", "errn-errp"]
     )
     result = estimator.run(fermi_dataset)
 
@@ -277,6 +277,41 @@ def test_compute_ts_map_downsampled(input_dataset):
     assert np.isnan(result["ts"].data[0, 30, 40])
 
 
+def test_ts_map_stat_scan(fake_dataset):
+    model = fake_dataset.models["source"]
+
+    dataset = fake_dataset.downsample(25)
+
+    estimator_ref = TSMapEstimator(
+        model,
+        kernel_width="0.3 deg",
+        energy_edges=[200, 3500] * u.GeV,
+    )
+
+    estimator = TSMapEstimator(
+        model,
+        kernel_width="0.3 deg",
+        selection_optional=["stat_scan"],
+        energy_edges=[200, 3500] * u.GeV,
+    )
+
+    maps_ref = estimator_ref.run(dataset)
+    maps = estimator.run(dataset)
+    success = maps.success.data
+
+    assert maps.stat_scan.geom.data_shape == (1, 109, 2, 2)
+    ts = np.abs(maps["stat_scan"].data.min(axis=1))
+    assert_allclose(ts[success], maps_ref.ts.data[success], rtol=1e-3)
+
+    ind_best = maps.stat_scan.data.argmin(axis=1)
+    ij, ik, il = np.indices(ind_best.shape)
+
+    dnde_ref = maps.dnde_ref.squeeze()
+    assert maps.dnde_scan_values.unit == dnde_ref.unit
+    norm = maps.dnde_scan_values.data[ij, ind_best, ik, il] / dnde_ref.value
+    assert_allclose(norm[success], maps_ref.norm.data[success], rtol=1e-5)
+
+
 def test_ts_map_with_model(fake_dataset):
     model = fake_dataset.models["source"]
 
@@ -285,7 +320,7 @@ def test_ts_map_with_model(fake_dataset):
     estimator = TSMapEstimator(
         model,
         kernel_width="0.3 deg",
-        selection_optional=["all"],
+        selection_optional=["ul", "errn-errp"],
         energy_edges=[200, 3500] * u.GeV,
     )
     maps = estimator.run(fake_dataset)