Repository: gammapy/gammapy
Pull request #5382 (2024-07-19 07:38:02)

Fix and tweaks for the stat_scan in TSMapEstimator
Fix and tweaks for the stat_scan in TSMapEstimator.
Follow up of #5285.


Patch:
diff --git a/gammapy/estimators/map/ts.py b/gammapy/estimators/map/ts.py
index 4ee810a10..57f05ca02 100644
--- a/gammapy/estimators/map/ts.py
+++ b/gammapy/estimators/map/ts.py
@@ -850,22 +850,24 @@ class BrentqFluxEstimator(Estimator):
         spline = InterpolatedUnivariateSpline(
             sparse_norms, stat_scan_local, k=1, ext="raise", check_finite=True
         )
-        stat_scan = spline(self.norm.scan_values)
+
+        norms = np.unique(np.concatenate((sparse_norms, self.norm.scan_values)))
+        stat_scan = spline(norms)
 
         ts = -stat_scan.min()
         ind = stat_scan.argmin()
-        norm = self.norm.scan_values[ind]
+        norm = norms[ind]
 
-        maskp = self.norm.scan_values > norm
+        maskp = norms > 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
+        norm_errn = norm - norms[~maskp][ind]
+
+        ind = np.abs(stat_diff - self.n_sigma**2)[maskp].argmin()
+        norm_errp = norms[maskp][ind] - norm
 
-        ind = np.abs(stat_diff - self.n_sigma_ul**2)[~maskp].argmin()
-        norm_ul = self.norm.scan_values[~maskp][ind]
+        ind = np.abs(stat_diff - self.n_sigma_ul**2)[maskp].argmin()
+        norm_ul = norms[maskp][ind]
 
         norm_err = (norm_errn + norm_errp) / 2
 
diff --git a/gammapy/estimators/utils.py b/gammapy/estimators/utils.py
index f3df04d3b..b5d4c0cb2 100644
--- a/gammapy/estimators/utils.py
+++ b/gammapy/estimators/utils.py
@@ -849,7 +849,7 @@ def combine_flux_maps(maps, method="gaussian_errors", reference_model=None):
     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):
+def _generate_scan_values(power_min=-6, 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`.
 


Test patch:
diff --git a/gammapy/estimators/map/tests/test_ts.py b/gammapy/estimators/map/tests/test_ts.py
index af30cd396..7a531abb4 100644
--- a/gammapy/estimators/map/tests/test_ts.py
+++ b/gammapy/estimators/map/tests/test_ts.py
@@ -286,6 +286,7 @@ def test_ts_map_stat_scan(fake_dataset):
         model,
         kernel_width="0.3 deg",
         energy_edges=[200, 3500] * u.GeV,
+        selection_optional=["errn-errp", "ul"],
     )
 
     estimator = TSMapEstimator(
@@ -311,6 +312,17 @@ def test_ts_map_stat_scan(fake_dataset):
     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)
 
+    assert_allclose(maps.norm.data[success], maps_ref.norm.data[success], rtol=1e-4)
+    assert_allclose(
+        maps.norm_errn.data[success], maps_ref.norm_errn.data[success], rtol=5e-2
+    )
+    assert_allclose(
+        maps.norm_errp.data[success], maps_ref.norm_errp.data[success], rtol=5e-2
+    )
+    assert_allclose(
+        maps.norm_ul.data[success], maps_ref.norm_ul.data[success], rtol=5e-2
+    )
+
 
 def test_ts_map_with_model(fake_dataset):
     model = fake_dataset.models["source"]