Merge pull request #2346 from registerrier/add_spectrumdataset_create

Add `SpectrumDatasetOnOff.create()`

gammapy/spectrum/dataset.py

@@ -364,6 +364,8 @@ def create(cls, e_reco, e_true=None, reference_time="2000-01-01"):
         Empty containers are created with the correct geometry.
         counts, background and aeff are zero and edisp is diagonal.
 
+        The safe_mask is set to False in every bin.
+
         Parameters
         ----------
         e_reco : `~astropy.units.Quantity`
@@ -551,6 +553,48 @@ def fake(self, background_model, random_state="random-seed"):
         npred_off.data = random_state.poisson(npred_off.data)
         self.counts_off = npred_off
 
+    @classmethod
+    def create(cls, e_reco, e_true=None, reference_time="2000-01-01"):
+        """Creates empty SpectrumDatasetOnOff
+
+        Empty containers are created with the correct geometry.
+        counts, counts_off and aeff are zero and edisp is diagonal.
+
+        The safe_mask is set to False in every bin.
+
+        Parameters
+        ----------
+        e_reco : `~astropy.units.Quantity`
+            edges of counts vector
+        e_true : `~astropy.units.Quantity`
+            edges of effective area table. If not set use reco energy values. Default : None
+        reference_time : `~astropy.time.Time`
+            reference time of the dataset, Default is "2000-01-01"
+        """
+        if e_true is None:
+            e_true = e_reco
+        counts = CountsSpectrum(e_reco[:-1], e_reco[1:])
+        counts_off = CountsSpectrum(e_reco[:-1], e_reco[1:])
+        aeff = EffectiveAreaTable(e_true[:-1],e_true[1:], np.zeros(e_true[:-1].shape)*u.m**2)
+        edisp = EnergyDispersion.from_diagonal_response(e_true, e_reco)
+        mask_safe = np.zeros_like(counts.data, 'bool')
+        gti = GTI.create(u.Quantity([],'s'), u.Quantity([],'s'), reference_time)
+        livetime = gti.time_sum
+        acceptance = np.ones_like(counts.data, int)
+        acceptance_off = np.ones_like(counts.data, int)
+
+        return SpectrumDatasetOnOff(counts=counts,
+                                    counts_off=counts_off,
+                                    aeff=aeff,
+                                    edisp=edisp,
+                                    mask_safe=mask_safe,
+                                    acceptance=acceptance,
+                                    acceptance_off=acceptance_off,
+                                    livetime=livetime,
+                                    gti=gti,
+                                   )
+
+
     @classmethod
     def read(cls, filename):
         """Read from file

gammapy/spectrum/tests/test_dataset.py

@@ -122,6 +122,20 @@ def test_set_model(self):
     def test_str(self):
         assert "SpectrumDataset" in str(self.dataset)
 
+    def test_spectrumdataset_create(self):
+        e_reco = u.Quantity([0.1,1,10.], 'TeV')
+        e_true = u.Quantity([0.05, 0.5, 5, 20.], 'TeV')
+        empty_dataset = SpectrumDataset.create(e_reco, e_true)
+
+        assert empty_dataset.counts.total_counts == 0
+        assert empty_dataset.data_shape[0] == 2
+        assert empty_dataset.background.total_counts == 0
+        assert empty_dataset.background.energy.nbin == 2
+        assert empty_dataset.aeff.data.axis("energy").nbin == 3
+        assert empty_dataset.edisp.data.axis("e_reco").nbin == 2
+        assert empty_dataset.livetime.value == 0
+        assert len(empty_dataset.gti.table) == 0
+
 
 class TestSpectrumOnOff:
     """ Test ON OFF SpectrumDataset"""
@@ -154,17 +168,20 @@ def setup(self):
             obs_id="test",
         )
 
-    def test_spectrumdataset_create(self):
-        e_reco = u.Quantity([0.1, 1, 10.0], "TeV")
-        e_true = u.Quantity([0.05, 0.5, 5, 20.0], "TeV")
-        empty_dataset = SpectrumDataset.create(e_reco, e_true)
+
+    def test_spectrumdatasetonoff_create(self):
+        e_reco = u.Quantity([0.1,1,10.], 'TeV')
+        e_true = u.Quantity([0.05, 0.5, 5, 20.], 'TeV')
+        empty_dataset = SpectrumDatasetOnOff.create(e_reco, e_true)
 
         assert empty_dataset.counts.total_counts == 0
         assert empty_dataset.data_shape[0] == 2
-        assert empty_dataset.background.total_counts == 0
-        assert empty_dataset.background.energy.nbin == 2
-        assert empty_dataset.aeff.data.axis("energy").nbin == 3
-        assert empty_dataset.edisp.data.axis("e_reco").nbin == 2
+        assert empty_dataset.counts_off.total_counts == 0
+        assert empty_dataset.counts_off.energy.nbin == 2
+        assert_allclose(empty_dataset.acceptance_off,1)
+        assert_allclose(empty_dataset.acceptance,1)
+        assert empty_dataset.acceptance.shape[0] == 2
+        assert empty_dataset.acceptance_off.shape[0] == 2
         assert empty_dataset.livetime.value == 0
         assert len(empty_dataset.gti.table) == 0