Repository: gammapy/gammapy

Issue:
3D (Fermi) + 1D (LST/MAGIC PointSkyRegion-style) analysis. Global fit ignores the 1D. 
**Gammapy version**
0.20.1

Are you using the latest stable version of Gammapy? (see https://gammapy.org/news.html)
yes

**Bug description**
Issue discussed already in slack (and I should probably have opened this ticket back then, my bad). A 3D (e.g. Fermi) and a 1D ON/OFF datasets (PointSkyRegion type, e.g. LST or MAGIC) are analyzed together with N skymodels for the 3D (with spatial and spectral components, including those for the source of interest) and 1 skymodel (the same model as in the 3D for ONLY the source of interest, including spectral and spatial components). The result, after the likelihood fit, is that one of the datasets is completely ignored in the fit. However, no error or warning is raised, so I have no idea how to debug it. 

<img width="730" alt="Captura de Pantalla 2022-09-29 a las 11 35 47" src="https://user-images.githubusercontent.com/2607760/193009687-d5d1bee3-9a8e-4fa0-852e-d625918a63cd.png">

Partial workaround: 
1) define a skymodel for the 1D dataset consisting only on a spectral component, no spatial ones. 
2) The parameters of the spectral components are linked/fixed to the corresponding spectral component of the source model in the 3D dataset (lets call the source name in the 3D dataset XXX). 
3) The source name of the skymodel then needs to be different, e.g. renamed to XXX-2, otherwise gammapy returns an error saying that the names should be unique.
4) perform the fit, and then voila both datasets contribute to the result, the confidence band is much narrower.

<img width="738" alt="Captura de Pantalla 2022-09-29 a las 11 35 57" src="https://user-images.githubusercontent.com/2607760/193009783-7fc169bf-0e19-4e89-887f-4b8cd6335b7f.png">

However, this is not a perfect solution, because the FluxPointEstimator needs to be called for each dataset with their unique source names (even when it is the same source and both datasets should contribute to the energy ranges where both instruments are sensitive). So you have Fermi points on one hand and IACT points on the other. At the overlapping energy range (~20-300 GeV)  you may have 'duplicated' points with larger error bars (because only one dataset contributes at a time). 

**Expected behavior**

The first option, i.e. defining 'datasets' to include both the 3D and the 1D (PointSkyRegion) and then making the fit should work fine, with the only difference with respect to a 3D that for the 1D model the spatial component should not contribute to the likelihood. Or if there is any reason why it shouldn't, it should at least raise an error instead of silently ignoring one of the datasets.

Another possible solution would be to try to get the FluxPointEstimation to also use the two datasets to compute each flux point, but that would involve dealing with two datasets with a source with different names in each. If the global fit works this way, no reason why the flux point estimation proceedure, which basically re-optimizes the 'norm' in each bin, shouldn't. But it sounds more 'hacky' to me than just fixing the first option and allowing the minimizer to work with this 1D+3D datasets

**To Reproduce**
1) Build a 3D dataset (suggested from a Fermi analysis done with the Fermitools).
2) Build a 1D dataset (energy-dependent / PointSkyRegion style, I ignore if the standard CircleSkyRegion analysis reproduces the issue)
3) Assign all the spatial+spectral skymodels for each source in the field of the 3D dataset.
4) Assign the spatial+spectral for the source of interest to match the one in the 3D dataset. 
5) Build the datasets([1D,3D]) object.
6) Make the likelihood fit with the models

**Other information**
Any other information you think will be useful for us to fix the issue can go here.


Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/gammapy/gammapy/pull/4133?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) Report
> Merging [#4133](https://codecov.io/gh/gammapy/gammapy/pull/4133?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) (682daba) into [main](https://codecov.io/gh/gammapy/gammapy/commit/fa09ca357460d651a8348e72df447ec80649c286?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) (fa09ca3) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

```diff
@@           Coverage Diff           @@
##             main    #4133   +/-   ##
=======================================
  Coverage   93.48%   93.48%           
=======================================
  Files         164      164           
  Lines       21916    21919    +3     
=======================================
+ Hits        20488    20491    +3     
  Misses       1428     1428           
```


| [Impacted Files](https://codecov.io/gh/gammapy/gammapy/pull/4133?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) | Coverage Δ | |
|---|---|---|
| [gammapy/maps/region/geom.py](https://codecov.io/gh/gammapy/gammapy/pull/4133/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9tYXBzL3JlZ2lvbi9nZW9tLnB5) | `92.93% <100.00%> (+0.04%)` | :arrow_up: |
| [gammapy/modeling/models/cube.py](https://codecov.io/gh/gammapy/gammapy/pull/4133/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9tb2RlbGluZy9tb2RlbHMvY3ViZS5weQ==) | `90.49% <100.00%> (+0.02%)` | :arrow_up: |

:mega: We’re building smart automated test selection to slash your CI/CD build times. [Learn more](https://about.codecov.io/iterative-testing/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy)


Patch:
diff --git a/gammapy/datasets/map.py b/gammapy/datasets/map.py
index 624f81e81..fd94e17d8 100644
--- a/gammapy/datasets/map.py
+++ b/gammapy/datasets/map.py
@@ -1357,11 +1357,12 @@ class MapDataset(Dataset):
 
         counts = 0
         background, excess, sqrt_ts = np.nan, np.nan, np.nan
+
         if self.counts:
-            summed_stat = self._counts_statistic[mask].sum()
-            counts = summed_stat.n_on
+            counts = self.counts.data[mask].sum()
 
             if self.background:
+                summed_stat = self._counts_statistic[mask].sum()
                 background = summed_stat.n_bkg
                 excess = summed_stat.n_sig
                 sqrt_ts = summed_stat.sqrt_ts
diff --git a/gammapy/maps/region/geom.py b/gammapy/maps/region/geom.py
index d5172ee36..cd43c873e 100644
--- a/gammapy/maps/region/geom.py
+++ b/gammapy/maps/region/geom.py
@@ -17,7 +17,6 @@ from regions import (
     PixCoord,
     PointSkyRegion,
     RectanglePixelRegion,
-    RectangleSkyRegion,
     Regions,
     SkyRegion,
 )
@@ -208,6 +207,8 @@ class RegionGeom(Geom):
         """Check if a given map coordinate is contained in the region.
         Requires the `.region` attribute to be set.
 
+        For `PointSkyRegion` the method always returns true.
+
         Parameters
         ----------
         coords : tuple, dict, `MapCoord` or `~astropy.coordinates.SkyCoord`
@@ -225,21 +226,16 @@ class RegionGeom(Geom):
 
         coords = MapCoord.create(coords, frame=self.frame, axis_names=self.axes.names)
 
-        if isinstance(self.region, PointSkyRegion):
-            # the point region is approximated by 1x1 pixel here
-            region = RectangleSkyRegion(
-                center=self.center_skydir,
-                width=self.binsz_wcs[0],
-                height=self.binsz_wcs[1],
-            )
-        else:
-            region = self.region
+        if self.is_all_point_sky_regions:
+            return np.ones(coords.skycoord.shape, dtype=bool)
 
-        return region.contains(coords.skycoord, self.wcs)
+        return self.region.contains(coords.skycoord, self.wcs)
 
     def contains_wcs_pix(self, pix):
         """Check if a given wcs pixel coordinate is contained in the region.
 
+        For `PointSkyRegion` the method always returns true.
+
         Parameters
         ----------
         pix : tuple
@@ -250,6 +246,9 @@ class RegionGeom(Geom):
         containment : `~numpy.ndarray`
             Bool array.
         """
+        if self.is_all_point_sky_regions:
+            return np.ones(pix[0].shape, dtype=bool)
+
         region_pix = self.region.to_pixel(self.wcs)
         return region_pix.contains(PixCoord(pix[0], pix[1]))
 
@@ -540,14 +539,7 @@ class RegionGeom(Geom):
         if self.region is None:
             pix = (0, 0)
         else:
-            # TODO: remove once fix is available in regions
-            if isinstance(self.region, PointSkyRegion):
-                point_region = self.region.to_pixel(self.wcs)
-                point_region.meta["include"] = False
-                pix_coord = PixCoord.from_sky(coords.skycoord, self.wcs)
-                in_region = point_region.contains(pix_coord)
-            else:
-                in_region = self.region.contains(coords.skycoord, wcs=self.wcs)
+            in_region = self.contains(coords.skycoord)
 
             x = np.zeros(coords.skycoord.shape)
             x[~in_region] = np.nan
diff --git a/gammapy/modeling/models/cube.py b/gammapy/modeling/models/cube.py
index d50060fb9..9ce56f9d1 100644
--- a/gammapy/modeling/models/cube.py
+++ b/gammapy/modeling/models/cube.py
@@ -271,6 +271,9 @@ class SkyModel(ModelBase):
             mask = mask.reduce_over_axes(func=np.logical_or)
 
         if mask.geom.is_region and mask.geom.region is not None:
+            if mask.geom.is_all_point_sky_regions:
+                return True
+
             geom = mask.geom.to_wcs_geom()
             mask = geom.region_mask([mask.geom.region])
 


Test patch:
diff --git a/gammapy/makers/tests/test_reduce.py b/gammapy/makers/tests/test_reduce.py
index 7cff459ff..0ae6adf76 100644
--- a/gammapy/makers/tests/test_reduce.py
+++ b/gammapy/makers/tests/test_reduce.py
@@ -307,7 +307,7 @@ def test_dataset_maker_spectrum_global_rad_max():
     assert counts.unit == ""
     assert counts_off.unit == ""
     assert_allclose(counts.data.sum(), 438, rtol=1e-5)
-    assert_allclose(counts_off.data.sum(), 273, rtol=1e-5)
+    assert_allclose(counts_off.data.sum(), 276, rtol=1e-5)
 
 
 @requires_data()
diff --git a/gammapy/maps/region/tests/test_ndmap.py b/gammapy/maps/region/tests/test_ndmap.py
index 1cf11d8af..e0d332084 100644
--- a/gammapy/maps/region/tests/test_ndmap.py
+++ b/gammapy/maps/region/tests/test_ndmap.py
@@ -32,6 +32,19 @@ def region_map():
     return m
 
 
+@pytest.fixture
+def region_map_no_region():
+    axis = MapAxis.from_energy_bounds("1 TeV", "10 TeV", nbin=6, name="energy")
+    m = Map.create(
+        region=None,
+        map_type="region",
+        axes=[axis],
+        unit="1/TeV",
+    )
+    m.data = np.arange(m.data.size, dtype=float).reshape(m.geom.data_shape)
+    return m
+
+
 @pytest.fixture
 def point_region_map():
     axis = MapAxis.from_energy_bounds("1 TeV", "10 TeV", nbin=6, name="energy")
@@ -211,10 +224,13 @@ def test_region_nd_map_get(region_map):
     values = region_map.get_by_coord((83.63, 21.51, energies[[0, -1]]))
     assert_allclose(values.squeeze(), [0, 5])
 
-    values = region_map.get_by_coord((energies[[0, -1]],))
+
+def test_region_nd_map_get_no_region(region_map_no_region):
+    energies = region_map_no_region.geom.axes[0].center
+    values = region_map_no_region.get_by_coord((energies[[0, -1]],))
     assert_allclose(values.squeeze(), [0, 5])
 
-    values = region_map.get_by_coord({"energy": energies[[0, -1]]})
+    values = region_map_no_region.get_by_coord({"energy": energies[[0, -1]]})
     assert_allclose(values.squeeze(), [0, 5])
 
 
diff --git a/gammapy/modeling/models/tests/test_cube.py b/gammapy/modeling/models/tests/test_cube.py
index 66e486a9a..279b38b3e 100644
--- a/gammapy/modeling/models/tests/test_cube.py
+++ b/gammapy/modeling/models/tests/test_cube.py
@@ -11,7 +11,7 @@ from regions import CircleSkyRegion
 from gammapy.data.gti import GTI
 from gammapy.datasets.map import MapEvaluator
 from gammapy.irf import EDispKernel, PSFKernel
-from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
+from gammapy.maps import Map, MapAxis, RegionGeom, RegionNDMap, WcsGeom
 from gammapy.modeling import Parameter
 from gammapy.modeling.models import (
     CompoundSpectralModel,
@@ -712,3 +712,11 @@ def test_compound_spectral_model(caplog):
         name="source-1",
     )
     assert_allclose(m.spectral_model(5 * u.TeV).value, 2.87e-12, rtol=1e-2)
+
+
+def test_sky_model_contributes_point_region():
+    model = SkyModel.create("pl", "point")
+
+    geom = RegionGeom.create("icrs;point(0.05, 0.05)", binsz_wcs="0.01 deg")
+    mask = RegionNDMap.from_geom(geom)
+    assert np.any(model.contributes(mask))