From 9a31a65a22fb5588688d89f05e129f3b3787d807 Mon Sep 17 00:00:00 2001
From: Terrier <rterrier@apc.in2p3.fr>
Date: Tue, 24 Sep 2019 15:11:44 +0200
Subject: [PATCH] Added MapDataset.to_spectrum_dataset

---
 gammapy/cube/fit.py            | 53 +++++++++++++++++++++++++++++++++-
 gammapy/cube/tests/test_fit.py | 14 +++++++++
 gammapy/spectrum/__init__.py   |  2 +-
 gammapy/spectrum/make.py       |  2 +-
 4 files changed, 68 insertions(+), 3 deletions(-)

diff --git a/gammapy/cube/fit.py b/gammapy/cube/fit.py
index 732d8b73e33..894d939a150 100644
--- a/gammapy/cube/fit.py
+++ b/gammapy/cube/fit.py
@@ -9,11 +9,12 @@
 from gammapy.cube.psf_kernel import PSFKernel
 from gammapy.cube.psf_map import PSFMap
 from gammapy.data import GTI
-from gammapy.irf import EnergyDispersion
+from gammapy.irf import EnergyDispersion, EffectiveAreaTable
 from gammapy.maps import Map, MapAxis
 from gammapy.modeling import Dataset, Parameters
 from gammapy.modeling.models import BackgroundModel, SkyModel, SkyModels, SkyPointSource
 from gammapy.stats import cash, cash_sum_cython, cstat, cstat_sum_cython
+from gammapy.spectrum import SpectrumDataset
 from gammapy.utils.random import get_random_state
 from gammapy.utils.scripts import make_path
 
@@ -683,6 +684,56 @@ def to_dict(self, filename=""):
         data["filename"] = filename
         return data
 
+    def to_spectrum_dataset(self, on_region):
+        """Return a ~gammapy.spectrum.SpectrumDataset from on_region.
+
+        Counts and background are summed in the on_region.
+
+        Effective area is taken from the average exposure divided by the livetime.
+        Here we assume it is the sum of the GTIs.
+
+        EnergyDispersion is obtained at the on_region center
+
+        Parameters
+        ----------
+        on_region : `~regions.SkyRegion`
+            the input ON region on which to extract the spectrum
+        Returns
+        -------
+        dataset : `~gammapy.spectrum.SpectrumDataset`
+            the resulting reduced dataset
+        """
+        if self.gti is not None:
+            livetime = self.gti.time_sum
+        else:
+            raise ValueError("No GTI in MapDataset. Impossible to compute livetime")
+
+        if self.counts is not None:
+            counts = self.counts.get_spectrum(on_region, np.sum)
+        else:
+            counts = None
+
+        #TODO : change to proper evaluation of BackgroundModel
+        if self.background_model is not None:
+            background = self.background_model.map.get_spectrum(on_region, np.sum)
+        else:
+            background = None
+
+        if self.exposure is not None:
+            exposure = self.exposure.get_spectrum(on_region, np.mean)
+            aeff = EffectiveAreaTable(
+                energy_lo=exposure.energy.edges[:-1], energy_hi=exposure.energy.edges[1:],
+                data=exposure.data/livetime
+            )
+        else:
+            aeff = None
+
+        # TODO: add obs_id once it is correctly defined in the MapDataset
+        return SpectrumDataset(counts=counts,
+                               background=background,
+                               aeff=aeff,
+                               livetime=livetime,
+                               gti=self.gti)
 
 class MapEvaluator:
     """Sky model evaluation on maps.
diff --git a/gammapy/cube/tests/test_fit.py b/gammapy/cube/tests/test_fit.py
index a6a54078e1b..b6e9e423966 100644
--- a/gammapy/cube/tests/test_fit.py
+++ b/gammapy/cube/tests/test_fit.py
@@ -5,6 +5,7 @@
 import astropy.units as u
 from astropy.coordinates import SkyCoord
 from regions import CircleSkyRegion
+from gammapy.data import GTI
 from gammapy.cube import MapDataset, PSFKernel, make_map_exposure_true_energy
 from gammapy.irf import (
     EffectiveAreaTable2D,
@@ -152,6 +153,19 @@ def test_different_exposure_unit(sky_model, geom):
 
     assert_allclose(npred.data[0, 50, 50], npred_ref.data[0, 50, 50])
 
+@requires_data()
+def test_to_spectrum_dataset(sky_model, geom):
+    dataset_ref = get_map_dataset(sky_model, geom, geom, edisp=True)
+    dataset_ref.counts = dataset_ref.background_model.map*0.
+    dataset_ref.counts.data[1,50,50]=1
+    dataset_ref.counts.data[1,60,50]=1
+
+    gti = GTI.create( [0*u.s], [10000*u.s], reference_time = "2010-01-01T00:00:00")
+    dataset_ref.gti = gti
+    on_region = CircleSkyRegion(center = geom.center_skydir, radius=0.1*u.deg)
+    spectrum_dataset = dataset_ref.to_spectrum_dataset(on_region)
+
+    assert np.sum(spectrum_dataset.counts.data) == 1
 
 @requires_data()
 def test_map_dataset_fits_io(tmpdir, sky_model, geom, geom_etrue):
diff --git a/gammapy/spectrum/__init__.py b/gammapy/spectrum/__init__.py
index 9812452fc13..cada65c62db 100644
--- a/gammapy/spectrum/__init__.py
+++ b/gammapy/spectrum/__init__.py
@@ -8,4 +8,4 @@
 from .phase import *
 from .reflected import *
 from .sensitivity import *
-from .make import *
\ No newline at end of file
+#from .make import *
\ No newline at end of file
diff --git a/gammapy/spectrum/make.py b/gammapy/spectrum/make.py
index 5ce72a7e094..a9cac4647f5 100644
--- a/gammapy/spectrum/make.py
+++ b/gammapy/spectrum/make.py
@@ -164,7 +164,7 @@ def extract_bkg(self):
 
         maps = self.mapmaker.run(["background"])
         mask = self.cutout.to_image().region_mask([self.on_region])
-        bkg_data = maps.background_model.quantity[..., mask].sum(axis=1)
+        bkg_data = maps.background_model.map.get_spectrum(self.on_region, np.sum) #quantity[..., mask].sum(axis=1)
 
         self.dataset.background = CountsSpectrum(
             energy_lo=self.e_reco[:-1], energy_hi=self.e_reco[1:], data=bkg_data