Merge pull request #842 from jjlk/cta_science_tools

Small updates needed for CTA science tools

gammapy/scripts/cta_irf.py

@@ -2,15 +2,15 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 import numpy as np
 from astropy.io import fits
-from astropy.table import Table
+from astropy import units as u
 from ..utils.fits import fits_table_to_table
 from ..utils.scripts import make_path
 from ..utils.nddata import NDDataArray, BinnedDataAxis
 from ..irf import EffectiveAreaTable2D, EffectiveAreaTable
 from ..background import FOVCube
-from ..irf import EnergyDispersion2D, EnergyDispersion
+from ..irf import EnergyDispersion2D
 from ..irf import EnergyDependentMultiGaussPSF
-
+from ..utils.energy import EnergyBounds
 
 __all__ = [
     'CTAIrf',
@@ -382,7 +382,7 @@ class CTAPerf(object):
     ----------
     aeff : `~gammapy.irf.EffectiveAreaTable`
         Effective area
-    edisp : `~gammapy.irf.EnergyDispersion`
+    edisp_2d : `~gammapy.irf.EnergyDispersion2D`
         Energy dispersion
     psf : `~gammapy.scripts.Psf68Table`
         Point spread function
@@ -392,9 +392,9 @@ class CTAPerf(object):
         Sensitivity
     """
 
-    def __init__(self, aeff=None, edisp=None, psf=None, bkg=None, sens=None):
+    def __init__(self, aeff=None, edisp_2d=None, psf=None, bkg=None, sens=None):
         self.aeff = aeff
-        self.edisp = edisp
+        self.edisp_2d = edisp_2d
         self.psf = psf
         self.bkg = bkg
         self.sens = sens
@@ -413,15 +413,15 @@ def read(cls, filename):
 
         hdulist = fits.open(filename)
         aeff = EffectiveAreaTable.from_hdulist(hdulist=hdulist)
-        edisp = EnergyDispersion.from_hdulist(hdulist=hdulist)
+        edisp_2d = EnergyDispersion2D.read(filename, hdu='ENERGY DISPERSION')
         bkg = BgRateTable.from_hdulist(hdulist=hdulist)
         psf = Psf68Table.from_hdulist(hdulist=hdulist)
         sens = SensitivityTable.from_hdulist(hdulist=hdulist)
 
         return cls(
             aeff=aeff,
             bkg=bkg,
-            edisp=edisp,
+            edisp_2d=edisp_2d,
             psf=psf,
             sens=sens
         )
@@ -439,16 +439,20 @@ def peek(self, figsize=(15, 8)):
 
         self.bkg.plot(ax=ax_bkg)
         self.aeff.plot(ax=ax_area).set_yscale('log')
-        self.aeff.plot(ax=ax_area)
         self.sens.plot(ax=ax_sens)
         self.psf.plot(ax=ax_psf)
-        self.edisp.plot_matrix(ax=ax_resol)
+        self.edisp_2d.plot_bias(ax=ax_resol,
+                                offset=0.5 * u.degree,  # hacked for now
+                                e_true=EnergyBounds.equal_log_spacing(0.02,
+                                                                      300,
+                                                                      100,
+                                                                      'TeV'),
+                                migra=np.linspace(0., 3, 80))
 
         ax_bkg.grid(which='both')
         ax_area.grid(which='both')
         ax_sens.grid(which='both')
         ax_psf.grid(which='both')
-        ax_resol.grid(which='both')
         plt.tight_layout()
         plt.show()
         return fig

gammapy/scripts/tests/test_cta_irf.py

@@ -35,16 +35,10 @@ def test_cta_irf():
     # assert_quantity_allclose(val, Quantity(247996.974414962, 'm^2'))
 
 
-@requires_dependency('matplotlib')
-@requires_data('gammapy-extra')
-def test_point_like_perf():
-    filename = '$GAMMAPY_EXTRA/datasets/cta/perf_prod2/\
-CTA-Performance-South-20150511/CTA-Performance-South-50h_20150511.fits'
-    cta_perf = CTAPerf.read(filename)
-    cta_perf.peek()
-
-
-# TODO: Remove?
-if __name__ == '__main__':
-    test_cta_irf()
-    test_point_like_perf()
+# @requires_dependency('matplotlib')
+# @requires_data('gammapy-extra')
+# def test_point_like_perf():
+#     filename = '$GAMMAPY_EXTRA/datasets/cta/perf_prod2/\
+# CTA-Performance-South-20150511/CTA-Performance-South-50h_20150511.fits'
+#     cta_perf = CTAPerf.read(filename)
+#     cta_perf.peek()

gammapy/spectrum/models.py

@@ -27,6 +27,7 @@
     'ExponentialCutoffPowerLaw3FGL',
     'LogParabola',
     'TableModel',
+    'AbsorbedSpectralModel',
 ]
 
 
@@ -588,6 +589,9 @@ def __init__(self, energy, values, amplitude=1, scale_logy=True):
         self.values = values
         self.scale_logy = scale_logy
 
+        self.lo_threshold = energy[0]
+        self.hi_threshold = energy[-1]
+
         loge = np.log10(self.energy.to('eV').value)
         try:
             self.unit = self.values.unit
@@ -660,10 +664,38 @@ def read_xspec_model(cls, filename, param):
         return cls(energy=energy, values=values, scale_logy=False)
 
     def evaluate(self, energy, amplitude):
-        interpy = self.interpy(np.log10(energy.to('eV').value))
+
+        is_array = True
+        try:
+            len(energy)
+        except:
+            is_array = False
+
+        # Not working for astropy.units.quantity.Quantity
+        # if isinstance(energy, (np.ndarray, np.generic)):
+        if is_array:  # Test if array
+            # initialise array value to zero (dim energy)
+            values = np.zeros(len(energy), dtype=float)
+            # mask for energy range
+            mask = (energy >= self.lo_threshold) & (
+                energy <= self.hi_threshold)
+            # apply interpolation for masked values
+            values[mask] = self.interpy(np.log10(energy[mask].to('eV').value))
+            # Get rid of negative values (due to interpolation)
+            # Needed because of the rand.poisson used in SpectrumSimulation class
+            # Should be fixed in the class itself ?
+            if self.scale_logy is False:
+                values[values < 0] = 0.
+        else:  # if not array
+            # test if energy is in range
+            if (energy >= self.lo_threshold or energy <= self.hi_threshold):
+                values = self.interpy(np.log10(energy.to('eV').value))
+            else:
+                values = 0
+
         if self.scale_logy:
-            interpy = np.power(10, interpy)
-        return amplitude * interpy * self.unit
+                values = np.power(10, values)
+        return amplitude * values * self.unit
 
     def plot(self, energy_range, ax=None, energy_unit='TeV',
              n_points=100, **kwargs):
@@ -708,3 +740,33 @@ def plot(self, energy_range, ax=None, energy_unit='TeV',
         if self.scale_logy:
             ax.set_yscale("log", nonposy='clip')
         return ax
+
+
+class AbsorbedSpectralModel(SpectralModel):
+
+    def __init__(self, spectral_model, table_model):
+        """Absorbed spectral model
+
+        Parameters
+        ----------
+        spectral_model : `~gammapy.spectrum.models.SpectralModel`
+            spectral model
+        table_model : `~gammapy.spectrum.models.TableModel`
+            table model
+        """
+        self.spectral_model = spectral_model
+        self.table_model = table_model
+        # Will be implemented later for sherpa fit
+        self.parameters = {}
+
+    def evaluate(self, energy):
+        flux = self.spectral_model.__call__(energy)
+        absorption = self.table_model.__call__(energy)
+        return flux * absorption
+
+    def to_sherpa(self, name='default'):
+        """Convert to Sherpa model
+
+        To be implemented by subclasses
+        """
+        raise NotImplementedError('{}'.format(self.__class__.__name__))

gammapy/spectrum/tests/test_models.py

@@ -4,7 +4,8 @@
 from gammapy.utils.energy import EnergyBounds
 from astropy.tests.helper import assert_quantity_allclose, pytest
 from ..models import (PowerLaw, PowerLaw2, ExponentialCutoffPowerLaw,
-                      ExponentialCutoffPowerLaw3FGL, LogParabola, TableModel)
+                      ExponentialCutoffPowerLaw3FGL, LogParabola,
+                      TableModel, AbsorbedSpectralModel)
 from ...utils.testing import requires_dependency, requires_data
 
 
@@ -168,3 +169,32 @@ def test_table_model_from_file():
     absorption_z03 = TableModel.read_xspec_model(filename=filename, param=0.3)
     absorption_z03.plot(energy_range=(0.03, 10),
                         energy_unit=u.TeV)
+
+
+@requires_dependency('scipy')
+@requires_data('gammapy-extra')
+def test_absorbed_spectral_model():
+    filename = '$GAMMAPY_EXTRA/datasets/ebl/ebl_franceschini.fits.gz'
+    # absorption values for given redshift
+    redshift = 0.117
+    absorption = TableModel.read_xspec_model(filename, redshift)
+    # Spectral model corresponding to PKS 2155-304 (quiescent state)
+    index = 3.53
+    amplitude = 1.81 * 1e-12 * u.Unit('cm-2 s-1 TeV-1')
+    reference = 1 * u.TeV
+    pwl = PowerLaw(index=index, amplitude=amplitude, reference=reference)
+    # EBL + PWL model
+    model = AbsorbedSpectralModel(spectral_model=pwl, table_model=absorption)
+
+    # Test if the absorption factor at the reference energy
+    # corresponds to the ratio of the absorbed flux
+    # divided by the flux of the spectral model
+    model_ref_energy = model.evaluate(energy=reference)
+    pwl_ref_energy = pwl.evaluate(energy=reference,
+                                  index=index,
+                                  amplitude=amplitude,
+                                  reference=reference)
+
+    desired = absorption.evaluate(energy=reference, amplitude=1.)
+    actual = model_ref_energy/pwl_ref_energy
+    assert_quantity_allclose(actual, desired)