Merge pull request #1073 from joleroi/fix_no_edisp

Fix spectrum fit for case of no EDISP

gammapy/spectrum/fit.py

@@ -276,7 +276,10 @@ def _calc_statval_helper(self, obs, prediction):
             # Store the result of the profile likelihood as bkg prediction
             mu_bkg = stats.get_wstat_mu_bkg(**kwargs)
             prediction[1] = mu_bkg * obs.alpha
-            on_stat = stats.wstat(**kwargs)
+            on_stat_ = stats.wstat(**kwargs)
+            # The on_stat sometime contains nan values
+            # TODO: Handle properly
+            on_stat = np.nan_to_num(on_stat_)
             off_stat = np.zeros_like(on_stat)
         else:
             raise NotImplementedError('{}'.format(self.stat))
@@ -305,7 +308,7 @@ def _check_valid_fit(self):
         """Helper function to give useful error messages."""
         # Assume that settings are the same for all observations
         test_obs = self.obs_list[0]
-        irfs_exist = test_obs.aeff is not None and test_obs.edisp is not None
+        irfs_exist = test_obs.aeff is not None or test_obs.edisp is not None
         if self.forward_folded and not irfs_exist:
             raise ValueError('IRFs required for forward folded fit')
         if self.stat == 'wstat' and self.obs_list[0].off_vector is None:

gammapy/spectrum/observation.py

@@ -696,7 +696,7 @@ def stack_counts_spectrum(counts_spectrum_list):
         stacked_quality = np.ones(energy.nbins)
 
         for spec in counts_spectrum_list:
-            stacked_data += spec.counts_in_safe_range
+            stacked_data += spec.counts_in_safe_range.value
             temp = np.logical_and(stacked_quality,
                                   spec.quality)
             stacked_quality = np.array(temp, dtype=int)
@@ -718,7 +718,7 @@ def stack_backscal(self):
             bkscal_on_data = obs.on_vector._backscal_array.copy()
 
             bkscal_off_data = obs.off_vector._backscal_array.copy()
-            bkscal_off += (bkscal_on_data / bkscal_off_data) * obs.off_vector.counts_in_safe_range
+            bkscal_off += (bkscal_on_data / bkscal_off_data) * obs.off_vector.counts_in_safe_range.value
 
         stacked_bkscal_off = self.stacked_off_vector.data.data.value / bkscal_off
 

gammapy/spectrum/tests/test_fit.py

@@ -10,6 +10,7 @@
     requires_data,
 )
 from ...utils.random import get_random_state
+from ...irf import EffectiveAreaTable
 from ...datasets import gammapy_extra
 from ...spectrum import (
     PHACountsSpectrum,
@@ -257,6 +258,19 @@ def test_fit_method(self):
         assert_quantity_allclose(result.model.parametersr['index'].value,
                                  2.2395184727047788)
 
+    def test_no_edisp(self):
+        obs = self.obs_list[0]
+        # Bring aeff in RECO space
+        data = obs.aeff.data.evaluate(energy=obs.on_vector.energy.nodes)
+        obs.aeff = EffectiveAreaTable(data=data,
+                                      energy_lo=obs.on_vector.energy.lo,
+                                      energy_hi=obs.on_vector.energy.hi)
+        obs.edisp = None
+        fit = SpectrumFit(obs_list=obs, model=self.pwl)
+        fit.fit()
+        assert_quantity_allclose(fit.result[0].model.parameters['index'].value,
+                                 2.2960518556630887, atol=0.02)
+
     def test_ecpl_fit(self):
         fit = SpectrumFit(self.obs_list[0], self.ecpl)
         fit.fit()