Fix SensitivityEstimator after IRF API change

gammapy/spectrum/sensitivity.py

@@ -16,11 +16,15 @@ class SensitivityEstimator(object):
 
     Parameters
     ----------
-    irf : `~gammapy.scripts.CTAPerf`
-        IRF object
+    arf : `~gammapy.irf.EffectiveAreaTable`
+        1D effective area
+    rmf : `~gammapy.irf.EnergyDispersion`
+        energy dispersion table
+    bkg : `~gammapy.spectrum.CountsSpectrum`
+        the background array
     livetime : `~astropy.units.Quantity`
         Livetime (object with the units of time), e.g. 5*u.h
-    slope : float, optional
+    index : float, optional
         Index of the spectral shape (Power-law), should be positive (>0)
     alpha : float, optional
         On/OFF normalisation
@@ -31,42 +35,32 @@ class SensitivityEstimator(object):
     bkg_sys : float, optional
         Fraction of Background systematics relative to the number of ON counts
 
-    Examples
-    --------
-    Compute and plot a sensitivity curve for CTA::
-
-        from gammapy.irf import CTAPerf
-        from gammapy.spectrum import SensitivityEstimator
-
-        filename = '$GAMMAPY_DATA/cta/perf_prod2/point_like_non_smoothed/South_5h.fits.gz'
-        irf = CTAPerf.read(filename)
-        sensitivity_estimator = SensitivityEstimator(irf=irf, livetime='5h')
-        sensitivity_estimator.run()
-        print(sensitivity_estimator.results_table)
-
-    Further examples in :gp-extra-notebook:`cta_sensitivity` .
+    An example can be found in :gp-extra-notebook:`cta_sensitivity` .
 
     Notes
     -----
-    For the moment, only the differential point-like sensitivity is computed at a fixed offset.
     This class allows to determine for each reconstructed energy bin the flux associated to the number of gamma-ray
     events for which the significance is ``sigma``, and being larger than ``gamma_min`` and ``bkg_sys`` percent larger than the
     number of background events in the ON region.
     """
 
     def __init__(
         self,
-        irf,
+        arf,
+        rmf,
+        bkg,
         livetime,
-        slope=2.0,
+        index=2.0,
         alpha=0.2,
         sigma=5.0,
         gamma_min=10.0,
         bkg_sys=0.05,
     ):
-        self.irf = irf
+        self.arf = arf
+        self.rmf = rmf
+        self.bkg = bkg
         self.livetime = u.Quantity(livetime).to("s")
-        self.slope = slope
+        self.index = index
         self.alpha = alpha
         self.sigma = sigma
         self.gamma_min = gamma_min
@@ -84,9 +78,9 @@ def run(self):
 
         # TODO: let the user decide on energy binning
         # then integrate bkg model and gamma over those energy bins.
-        energy = self.irf.bkg.energy.log_center()
+        energy = self.rmf.e_reco.log_center()
 
-        bkg_counts = (self.irf.bkg.data.data * self.livetime).value
+        bkg_counts = (self.bkg.data.data.to('1/s') * self.livetime).value
 
         excess_counts = excess_matching_significance_on_off(
             n_off=bkg_counts / self.alpha, alpha=self.alpha, significance=self.sigma
@@ -95,12 +89,12 @@ def run(self):
         excess_counts[is_gamma_limited] = self.gamma_min
 
         model = PowerLaw(
-            index=self.slope, amplitude="1 cm-2 s-1 TeV-1", reference="1 TeV"
+            index=self.index, amplitude="1 cm-2 s-1 TeV-1", reference="1 TeV"
         )
 
         # TODO: simplify the following computation
         predictor = CountsPredictor(
-            model, aeff=self.irf.aeff, edisp=self.irf.rmf, livetime=self.livetime
+            model, aeff=self.arf, edisp=self.rmf, livetime=self.livetime
         )
         predictor.run()
         counts = predictor.npred.data.data.value
@@ -153,3 +147,4 @@ def run(self):
             ]
         )
         self._results_table = table
+        return table

gammapy/spectrum/tests/test_sensitivity.py

@@ -2,47 +2,55 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 import pytest
 from numpy.testing import assert_allclose
+import numpy as np
 import astropy.units as u
 from ...utils.testing import requires_data
+from ...irf import EffectiveAreaTable, EnergyDispersion
+from ..core import CountsSpectrum
 from ..sensitivity import SensitivityEstimator
 
 
 @pytest.fixture()
 def sens():
-    filename = "$GAMMAPY_EXTRA/datasets/cta/perf_prod2/point_like_non_smoothed/South_5h.fits.gz"
-    irf = CTAPerf.read(filename)
-    sens = SensitivityEstimator(irf=irf, livetime=5.0 * u.h)
+    etrue = np.logspace(0, 1, 21) * u.TeV
+    elo = etrue[:-1]
+    ehi = etrue[1:]
+    area = np.zeros(20) + 1e6 * u.m**2
+
+    arf = EffectiveAreaTable(energy_lo = elo, energy_hi = ehi, data = area)
+
+    ereco = np.logspace(0,1,5) * u.TeV
+    rmf = EnergyDispersion.from_diagonal_response(etrue, ereco)
+
+    bkg_array = np.ones(4)/u.s
+    bkg_array[-1] = 1e-3/u.s
+    bkg = CountsSpectrum(energy_lo=ereco[:-1], energy_hi=ereco[1:], data=bkg_array)
+
+    sens = SensitivityEstimator(arf=arf, rmf=rmf, bkg=bkg, livetime=1 * u.h, index=2, gamma_min=20, alpha=0.2)
     sens.run()
     return sens
 
 
-@pytest.mark.xfail
 @requires_data("gammapy-extra")
 def test_cta_sensitivity_estimator(sens):
     table = sens.results_table
 
-    assert len(table) == 21
+    assert len(table) == 4
     assert table.colnames == ["energy", "e2dnde", "excess", "background", "criterion"]
     assert table["energy"].unit == "TeV"
     assert table["e2dnde"].unit == "erg / (cm2 s)"
 
     row = table[0]
-    assert_allclose(row["energy"], 0.015848, rtol=1e-3)
-    assert_allclose(row["e2dnde"], 1.2656e-10, rtol=1e-3)
-    assert_allclose(row["excess"], 339.143, rtol=1e-3)
-    assert_allclose(row["background"], 3703.48, rtol=1e-3)
+    assert_allclose(row["energy"], 1.33352, rtol=1e-3)
+    assert_allclose(row["e2dnde"], 3.40101e-11, rtol=1e-3)
+    assert_allclose(row["excess"], 334.454, rtol=1e-3)
+    assert_allclose(row["background"], 3600, rtol=1e-3)
     assert row["criterion"] == "significance"
 
-    row = table[9]
-    assert_allclose(row["energy"], 1, rtol=1e-3)
-    assert_allclose(row["e2dnde"], 4.28759e-13, rtol=1e-3)
-    assert_allclose(row["excess"], 18.1072, rtol=1e-3)
-    assert_allclose(row["background"], 5.11857, rtol=1e-3)
-    assert row["criterion"] == "significance"
-
-    row = table[20]
-    assert_allclose(row["energy"], 158.489, rtol=1e-3)
-    assert_allclose(row["e2dnde"], 9.0483e-12, rtol=1e-3)
-    assert_allclose(row["excess"], 10, rtol=1e-3)
-    assert_allclose(row["background"], 0.00566093, rtol=1e-3)
+    row = table[3]
+    assert_allclose(row["energy"], 7.49894, rtol=1e-3)
+    assert_allclose(row["e2dnde"], 1.14367e-11, rtol=1e-3)
+    assert_allclose(row["excess"], 20, rtol=1e-3)
+    assert_allclose(row["background"], 3.6, rtol=1e-3)
     assert row["criterion"] == "gamma"
+