apply containment correction in true energy

gammapy/spectrum/extract.py

@@ -69,6 +69,12 @@ def __init__(self, obs_list, bkg_estimate, e_reco=None, e_true=None,
         self.use_recommended_erange = use_recommended_erange
         self.observations = SpectrumObservationList()
 
+        self.containment = None
+        self._on_vector = None
+        self._off_vector = None
+        self._aeff = None
+        self._edisp = None
+
     def run(self, outdir=None, use_sherpa=False):
         """Run all steps.
 
@@ -119,6 +125,8 @@ def process(self, obs, bkg):
 
         if self.containment_correction:
             self.apply_containment_correction(obs, bkg)
+        else:
+            self.containment = np.ones(self._aeff.energy.nbins)
 
         spectrum_observation = SpectrumObservation(
             on_vector=self._on_vector,
@@ -218,22 +226,20 @@ def apply_containment_correction(self, obs, bkg):
         else:
             psf = obs.psf.to_energy_dependent_table_psf(offset, angles)
 
-        center_energies = self._on_vector.energy.nodes
-        areascal = []
+        center_energies = self._aeff.energy.nodes
+        containment = []
         for index, energy in enumerate(center_energies):
             try:
-                correction = psf.integral(energy,
-                                          0. * u.deg,
-                                          bkg.on_region.radius)
+                cont_ = psf.integral(energy, 0. * u.deg, bkg.on_region.radius)
             except:
                 msg = 'Containment correction failed for bin {}, energy {}.'
                 log.warning(msg.format(index, energy))
-                correction = 1
+                cont_ = 1
             finally:
-                areascal.append(correction)
+                containment.append(cont_)
 
-        self._on_vector.areascal = areascal
-        self._off_vector.areascal = areascal
+        self.containment = np.array(containment)
+        self._aeff.data.data *= self.containment
 
     def compute_energy_threshold(self, **kwargs):
         """Compute and set the safe energy threshold for all observations.

gammapy/spectrum/tests/test_extract.py

@@ -43,9 +43,9 @@ class TestSpectrumExtraction:
                                                   )),
         (dict(containment_correction=True), dict(n_on=172,
                                                  sigma=24.98,
-                                                 aeff=549861.8 * u.m ** 2,
+                                                 aeff=412731.8 * u.m ** 2,
                                                  edisp=0.2595896944765074,
-                                                 containment=0.8525390663792395,
+                                                 containment=0.7645777148101338,
                                                  ))
     ])
     def test_extract(self, pars, results, obs_list, bkg_estimate):
@@ -63,10 +63,11 @@ def test_extract(self, pars, results, obs_list, bkg_estimate):
         assert_quantity_allclose(aeff_actual, results['aeff'], rtol=1e-3)
         assert_quantity_allclose(edisp_actual, results['edisp'], rtol=1e-3)
 
+        containment_actual = extraction.containment[60]
+
         # TODO: Introduce assert_stats_allclose
         n_on_actual = obs.total_stats.n_on
         sigma_actual = obs.total_stats.sigma
-        containment_actual = obs.on_vector.areascal[60]
 
         assert n_on_actual == results['n_on']
         assert_allclose(sigma_actual, results['sigma'], atol=1e-2)

gammapy/spectrum/tests/test_fit.py

@@ -224,16 +224,6 @@ def test_compound(self):
         assert_quantity_allclose(pars['amplitude'].quantity,
                                  1.0243449507421302e-7 * u.Unit('m-2 s-1 TeV-1'))
 
-    def test_areascal(self):
-        areascal = np.ones(self.fit.obs_list[0].e_reco.nbins)
-        areascal *= 0.5
-        self.fit.obs_list[0].on_vector.areascal = areascal
-        self.fit.fit()
-        pars = self.fit.result[0].model.parameters
-        assert_quantity_allclose(pars['index'].value, 2.2542312883476465)
-        assert_quantity_allclose(pars['amplitude'].quantity,
-                                 4.0165729155447114e-7 * u.Unit('m-2 s-1 TeV-1'))
-
     def test_npred(self):
         self.fit.fit()
         actual = self.fit.obs_list[0].predicted_counts(