Stack EDISP for a set of observations

gammapy/data/data_store.py

@@ -12,10 +12,11 @@
 from astropy.time import Time
 from astropy.coordinates import SkyCoord
 from ..utils.scripts import make_path
+from ..utils.energy import Energy
 from .obs_table import ObservationTable
 from .hdu_index_table import HDUIndexTable
 from .utils import _earth_location_from_dict
-from ..irf import EnergyDependentTablePSF
+from ..irf import EnergyDependentTablePSF, IRFStacker
 
 __all__ = [
     'DataStore',
@@ -668,6 +669,7 @@ class ObservationList(UserList):
 
     Could be extended to hold a more generic class of observations
     """
+
     def __str__(self):
         s = self.__class__.__name__ + '\n'
         s += 'Number of observations: {}\n'.format(len(self))
@@ -711,3 +713,59 @@ def make_psf(self, position, energy=None, theta=None):
         psf_tot = EnergyDependentTablePSF(energy=energy, offset=theta, exposure=exposure,
                                           psf_value=psf_value.T)
         return psf_tot
+
+    def make_mean_edisp(self, position, e_true, e_reco, low_reco_threshold=None, high_reco_threshold=None):
+        r"""Make mean rmf for a given position and a set of observations.
+
+        Compute the mean rmf of a set of observations j at a given position
+
+        The stacking of :math:`j` observations is implemented as follows.  :math:`k`
+        and :math:`l` denote a bin in reconstructed and true energy, respectively.
+
+        .. math::
+
+            \epsilon_{jk} =\left\{\begin{array}{cl} 1, & \mbox{if
+                bin k is inside the energy thresholds}\\ 0, & \mbox{otherwise} \end{array}\right.
+
+            \overline{\mathrm{edisp}}_{kl} = \frac{\sum_{j} \mathrm{edisp}_{jkl}
+                \cdot \mathrm{aeff}_{jl} \cdot t_j \cdot \epsilon_{jk}}{\sum_{j} \mathrm{aeff}_{jl}
+                \cdot t_j}
+
+        Parameters
+        ----------
+        position : `~astropy.coordinates.SkyCoord`
+            Position at which to compute the mean EDISP
+        e_true : `~gammapy.utils.energy.EnergyBounds`
+            True energy axis
+        e_reco : `~gammapy.utils.energy.EnergyBounds`
+            Reconstructed energy axis
+        low_reco_threshold : `~gammapy.utils.energy.Energy`
+            low energy threshold in reco energy, default 0.002 TeV
+        high_reco_threshold : `~gammapy.utils.energy.Energy`
+            high energy threshold in reco energy , default 150 TeV
+        Returns
+        -------
+        stacked_rmf: `~gammapy.irf.EnergyDispersion`
+            Stacked RMF for a set of observation
+        """
+
+        list_arf = list()
+        list_edisp = list()
+        list_livetime = list()
+        if not low_reco_threshold:
+            low_reco_threshold = Energy(0.002, "TeV")
+        if not high_reco_threshold:
+            high_reco_threshold = Energy(150, "TeV")
+        list_low_threshold = [low_reco_threshold] * len(self)
+        list_high_threshold = [high_reco_threshold] * len(self)
+        for obs in self:
+            offset = position.separation(obs.pointing_radec)
+            list_arf.append(obs.aeff.to_effective_area_table(offset, energy=e_true))
+            list_edisp.append(obs.edisp.to_energy_dispersion(offset, e_reco=e_reco, e_true=e_true))
+            list_livetime.append(obs.observation_live_time_duration)
+
+        irf_stack = IRFStacker(list_arf=list_arf, list_edisp=list_edisp, list_livetime=list_livetime,
+                               list_low_threshold=list_low_threshold, list_high_threshold=list_high_threshold)
+        irf_stack.mean_edisp()
+
+        return irf_stack.stacked_edisp

gammapy/data/tests/test_data_store.py

@@ -1,13 +1,14 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import absolute_import, division, print_function, unicode_literals
 import numpy as np
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_equal
 from astropy.tests.helper import pytest, assert_quantity_allclose
 from astropy.coordinates import Angle, SkyCoord
 from astropy.units import Quantity
 import astropy.units as u
-from ...data import DataStore, DataManager
+from ...data import DataStore, DataManager, ObservationList
 from ...utils.testing import data_manager, requires_data, requires_dependency
+from ...utils.energy import Energy
 from ...datasets import gammapy_extra
 from ...utils.energy import EnergyBounds
 
@@ -150,3 +151,38 @@ def test_make_psf(pars, result):
     assert_quantity_allclose(psf.energy[10], result["psf_energy"])
     assert_quantity_allclose(psf.exposure[10], result["psf_exposure"])
     assert_quantity_allclose(psf.psf_value[10, 50], result["psf_value"])
+
+
+def test_make_meanedisp(tmpdir):
+    position = SkyCoord(83.63, 22.01, unit='deg')
+    store = gammapy_extra.filename("datasets/hess-crab4-hd-hap-prod2")
+    data_store = DataStore.from_dir(store)
+
+    obs1 = data_store.obs(23523)
+    obs2 = data_store.obs(23592)
+    obslist = ObservationList([obs1, obs2])
+
+    e_true = EnergyBounds.equal_log_spacing(0.01, 150, 80, "TeV")
+    e_reco = EnergyBounds.equal_log_spacing(0.5, 100, 15, "TeV")
+    rmf = obslist.make_mean_edisp(position=position, e_true=e_true, e_reco=e_reco)
+
+    assert len(rmf.e_true.nodes) == 80
+    assert len(rmf.e_reco.nodes) == 15
+    assert_quantity_allclose(rmf.data[53, 8], 0.0559785805550798)
+
+    rmf2 = obslist.make_mean_edisp(position=position, e_true=e_true, e_reco=e_reco,
+                                 low_reco_threshold=Energy(1, "TeV"), high_reco_threshold=Energy(60, "TeV"))
+    # Test that the energy dispersion is equal to zero for the reco energy bin inferior to low_reco_threshold and
+    #  superior to high_reco_threshold
+    i2 = np.where(rmf2.data[:, 13] != 0)[0]
+    assert len(i2) == 0
+    i2 = np.where(rmf2.data[:, 1] != 0)[0]
+    assert len(i2) == 0
+    # Test that for the reco energy bin inside the thresholds, the edisp defined without threshold in the same than
+    # the edisp defined with the low and high reco energy thresholds
+    i = np.where(rmf.data[:, 12] != 0)[0]
+    i2 = np.where(rmf2.data[:, 12] != 0)[0]
+    assert_equal(i, i2)
+    i = np.where(rmf.data[:, 2] != 0)[0]
+    i2 = np.where(rmf2.data[:, 2] != 0)[0]
+    assert_equal(i, i2)

gammapy/irf/irf_stack.py

@@ -43,7 +43,7 @@ class IRFStacker(object):
         list of `~gammapy.irf.EffectiveAreaTable`
     list_livetime: list
         list of `~astropy.units.Quantity` (livetime)
-    list_rmf: list
+    list_edisp: list
         list of `~gammapy.irf.EnergyDispersion`
     list_low_threshold: list
         list of low energy threshold, optional for effective area mean computation
@@ -53,11 +53,11 @@ class IRFStacker(object):
 
     """
 
-    def __init__(self, list_arf, list_livetime, list_rmf=None,
+    def __init__(self, list_arf, list_livetime, list_edisp=None,
                  list_low_threshold=None, list_high_threshold=None):
         self.list_arf = list_arf
         self.list_livetime = Quantity(list_livetime)
-        self.list_rmf = list_rmf
+        self.list_edisp = list_edisp
         self.list_low_threshold = list_low_threshold
         self.list_high_threshold = list_high_threshold
         self.stacked_aeff = None
@@ -80,27 +80,27 @@ def stack_aeff(self):
         self.stacked_aeff = EffectiveAreaTable(energy=self.list_arf[0].energy,
                                                data=stacked_data.to('cm2'))
 
-    def mean_rmf(self):
+    def mean_edisp(self):
         """
         Compute mean energy dispersion
         """
-        reco_bins = self.list_rmf[0].e_reco.nbins
-        true_bins = self.list_rmf[0].e_true.nbins
+        reco_bins = self.list_edisp[0].e_reco.nbins
+        true_bins = self.list_edisp[0].e_true.nbins
 
         aefft = Quantity(np.zeros(true_bins), 'cm2 s')
         temp = np.zeros(shape=(reco_bins, true_bins))
         aefftedisp = Quantity(temp, 'cm2 s')
 
-        for i, rmf in enumerate(self.list_rmf):
+        for i, edisp in enumerate(self.list_edisp):
             aeff_data = self.list_arf[i].evaluate(fill_nan=True)
             aefft_current = aeff_data * self.list_livetime[i]
             aefft += aefft_current
-            edisp_data = rmf.pdf_in_safe_range(self.list_low_threshold[i],
-                                               self.list_high_threshold[i])
+            edisp_data = edisp.pdf_in_safe_range(self.list_low_threshold[i],
+                                                 self.list_high_threshold[i])
 
             aefftedisp += edisp_data.transpose() * aefft_current
 
         stacked_edisp = np.nan_to_num(aefftedisp / aefft)
-        self.stacked_edisp = EnergyDispersion(e_true=self.list_rmf[0].e_true,
-                                              e_reco=self.list_rmf[0].e_reco,
+        self.stacked_edisp = EnergyDispersion(e_true=self.list_edisp[0].e_true,
+                                              e_reco=self.list_edisp[0].e_reco,
                                               data=stacked_edisp.transpose())

gammapy/spectrum/observation.py

@@ -644,22 +644,22 @@ def stack_edisp(self):
         calls :func:`~gammapy.irf.IRFStacker.stack_edisp`
         """
         list_arf = list()
-        list_rmf = list()
+        list_edisp = list()
         list_livetime = list()
         list_elo_threshold = list()
         list_ehi_threshold = list()
         for o in self.obs_list:
             list_arf.append(o.aeff)
             list_livetime.append(o.livetime)
-            list_rmf.append(o.edisp)
+            list_edisp.append(o.edisp)
             list_elo_threshold.append(o.lo_threshold)
             list_ehi_threshold.append(o.hi_threshold)
         irf_stack = IRFStacker(list_arf=list_arf,
                                list_livetime=list_livetime,
-                               list_rmf=list_rmf,
+                               list_edisp=list_edisp,
                                list_low_threshold=list_elo_threshold,
                                list_high_threshold=list_ehi_threshold)
-        irf_stack.mean_rmf()
+        irf_stack.mean_edisp()
         self.stacked_edisp = irf_stack.stacked_edisp
 
     def stack_obs(self):