Source code for gammapy.spectrum.models

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Spectral models for Gammapy."""
from __future__ import absolute_import, division, print_function, unicode_literals
import operator
import numpy as np
from scipy.optimize import brentq
import astropy.units as u
from astropy.table import Table
from ..utils.energy import EnergyBounds
from ..utils.nddata import NDDataArray, BinnedDataAxis
from ..utils.scripts import make_path
from ..utils.fitting import Parameter, Parameters, Model
from ..utils.interpolation import ScaledRegularGridInterpolator
from .utils import integrate_spectrum

__all__ = [
    "SpectralModel",
    "ConstantModel",
    "CompoundSpectralModel",
    "PowerLaw",
    "PowerLaw2",
    "ExponentialCutoffPowerLaw",
    "ExponentialCutoffPowerLaw3FGL",
    "PLSuperExpCutoff3FGL",
    "LogParabola",
    "TableModel",
    "AbsorbedSpectralModel",
    "Absorption",
]


[docs]class SpectralModel(Model): """Spectral model base class. Derived classes should store their parameters as `~gammapy.utils.modeling.Parameters` See for example return pardict of `~gammapy.spectrum.models.PowerLaw`. """
[docs] def __call__(self, energy): """Call evaluate method of derived classes""" kwargs = dict() for par in self.parameters.parameters: kwargs[par.name] = par.quantity return self.evaluate(energy, **kwargs)
def __mul__(self, model): if not isinstance(model, SpectralModel): model = ConstantModel(const=model) return CompoundSpectralModel(self, model, operator.mul) def __rmul__(self, model): # This is needed to support e.g. 5 * model return self.__mul__(model) def __add__(self, model): if not isinstance(model, SpectralModel): model = ConstantModel(const=model) return CompoundSpectralModel(self, model, operator.add) def __radd__(self, model): return self.__add__(model) def __sub__(self, model): if not isinstance(model, SpectralModel): model = ConstantModel(const=model) return CompoundSpectralModel(self, model, operator.sub) def __rsub__(self, model): return self.__sub__(model) def __truediv__(self, model): if not isinstance(model, SpectralModel): model = ConstantModel(const=model) return CompoundSpectralModel(self, model, operator.truediv) def __rtruediv__(self, model): return self.__div__(model) def _parse_uarray(self, uarray): from uncertainties import unumpy values = unumpy.nominal_values(uarray) errors = unumpy.std_devs(uarray) return values, errors def _convert_energy(self, energy): if "reference" in self.parameters.names: return energy.to(self.parameters["reference"].unit) elif "emin" in self.parameters.names: return energy.to(self.parameters["emin"].unit) else: return energy
[docs] def evaluate_error(self, energy): """Evaluate spectral model with error propagation. Parameters ---------- energy : `~astropy.units.Quantity` Energy at which to evaluate Returns ------- flux, flux_error : tuple of `~astropy.units.Quantity` Tuple of flux and flux error. """ energy = self._convert_energy(energy) unit = self(energy).unit upars = self.parameters._ufloats uarray = self.evaluate(energy.value, **upars) return self._parse_uarray(uarray) * unit
[docs] def integral(self, emin, emax, **kwargs): r"""Integrate spectral model numerically. .. math:: F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} \phi(E) dE If array input for ``emin`` and ``emax`` is given you have to set ``intervals=True`` if you want the integral in each energy bin. Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. **kwargs : dict Keyword arguments passed to :func:`~gammapy.spectrum.integrate_spectrum` """ return integrate_spectrum(self, emin, emax, **kwargs)
[docs] def integral_error(self, emin, emax, **kwargs): """Integrate spectral model numerically with error propagation. Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower adn upper bound of integration range. **kwargs : dict Keyword arguments passed to func:`~gammapy.spectrum.integrate_spectrum` Returns ------- integral, integral_error : tuple of `~astropy.units.Quantity` Tuple of integral flux and integral flux error. """ emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.integral(emin, emax, **kwargs).unit upars = self.parameters._ufloats def f(x): return self.evaluate(x, **upars) uarray = integrate_spectrum(f, emin.value, emax.value, **kwargs) return self._parse_uarray(uarray) * unit
[docs] def energy_flux(self, emin, emax, **kwargs): r"""Compute energy flux in given energy range. .. math:: G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} E \phi(E) dE Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. **kwargs : dict Keyword arguments passed to func:`~gammapy.spectrum.integrate_spectrum` """ def f(x): return x * self(x) return integrate_spectrum(f, emin, emax, **kwargs)
[docs] def energy_flux_error(self, emin, emax, **kwargs): r"""Compute energy flux in given energy range with error propagation. .. math:: G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} E \phi(E) dE Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower bound of integration range. **kwargs : dict Keyword arguments passed to :func:`~gammapy.spectrum.integrate_spectrum` Returns ------- energy_flux, energy_flux_error : tuple of `~astropy.units.Quantity` Tuple of energy flux and energy flux error. """ emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.energy_flux(emin, emax, **kwargs).unit upars = self.parameters._ufloats def f(x): return x * self.evaluate(x, **upars) uarray = integrate_spectrum(f, emin.value, emax.value, **kwargs) return self._parse_uarray(uarray) * unit
[docs] def to_dict(self): """Convert to dict.""" retval = self.parameters.to_dict() retval["name"] = self.__class__.__name__ return retval
[docs] @classmethod def from_dict(cls, val): """Create from dict.""" classname = val.pop("name") parameters = Parameters.from_dict(val) model = globals()[classname]() model.parameters = parameters model.parameters.covariance = parameters.covariance return model
[docs] def plot( self, energy_range, ax=None, energy_unit="TeV", flux_unit="cm-2 s-1 TeV-1", energy_power=0, n_points=100, **kwargs ): """Plot spectral model curve. kwargs are forwarded to `matplotlib.pyplot.plot` By default a log-log scaling of the axes is used, if you want to change the y axis scaling to linear you can use:: from gammapy.spectrum.models import ExponentialCutoffPowerLaw from astropy import units as u pwl = ExponentialCutoffPowerLaw() ax = pwl.plot(energy_range=(0.1, 100) * u.TeV) ax.set_yscale('linear') Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Axis energy_range : `~astropy.units.Quantity` Plot range energy_unit : str, `~astropy.units.Unit`, optional Unit of the energy axis flux_unit : str, `~astropy.units.Unit`, optional Unit of the flux axis energy_power : int, optional Power of energy to multiply flux axis with n_points : int, optional Number of evaluation nodes Returns ------- ax : `~matplotlib.axes.Axes`, optional Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax emin, emax = energy_range energy = EnergyBounds.equal_log_spacing(emin, emax, n_points, energy_unit) # evaluate model flux = self(energy).to(flux_unit) y = self._plot_scale_flux(energy, flux, energy_power) ax.plot(energy.value, y.value, **kwargs) self._plot_format_ax(ax, energy, y, energy_power) return ax
[docs] def plot_error( self, energy_range, ax=None, energy_unit="TeV", flux_unit="cm-2 s-1 TeV-1", energy_power=0, n_points=100, **kwargs ): """Plot spectral model error band. .. note:: This method calls ``ax.set_yscale("log", nonposy='clip')`` and ``ax.set_xscale("log", nonposx='clip')`` to create a log-log representation. The additional argument ``nonposx='clip'`` avoids artefacts in the plot, when the error band extends to negative values (see also https://github.com/matplotlib/matplotlib/issues/8623). When you call ``plt.loglog()`` or ``plt.semilogy()`` explicitely in your plotting code and the error band extends to negative values, it is not shown correctly. To circumvent this issue also use ``plt.loglog(nonposx='clip', nonposy='clip')`` or ``plt.semilogy(nonposy='clip')``. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Axis energy_range : `~astropy.units.Quantity` Plot range energy_unit : str, `~astropy.units.Unit`, optional Unit of the energy axis flux_unit : str, `~astropy.units.Unit`, optional Unit of the flux axis energy_power : int, optional Power of energy to multiply flux axis with n_points : int, optional Number of evaluation nodes **kwargs : dict Keyword arguments forwarded to `matplotlib.pyplot.fill_between` Returns ------- ax : `~matplotlib.axes.Axes`, optional Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax kwargs.setdefault("facecolor", "black") kwargs.setdefault("alpha", 0.2) kwargs.setdefault("linewidth", 0) emin, emax = energy_range energy = EnergyBounds.equal_log_spacing(emin, emax, n_points, energy_unit) flux, flux_err = self.evaluate_error(energy).to(flux_unit) y_lo = self._plot_scale_flux(energy, flux - flux_err, energy_power) y_hi = self._plot_scale_flux(energy, flux + flux_err, energy_power) where = (energy >= energy_range[0]) & (energy <= energy_range[1]) ax.fill_between(energy.value, y_lo.value, y_hi.value, where=where, **kwargs) self._plot_format_ax(ax, energy, y_lo, energy_power) return ax
def _plot_format_ax(self, ax, energy, y, energy_power): ax.set_xlabel("Energy [{}]".format(energy.unit)) if energy_power > 0: ax.set_ylabel("E{} * Flux [{}]".format(energy_power, y.unit)) else: ax.set_ylabel("Flux [{}]".format(y.unit)) ax.set_xscale("log", nonposx="clip") ax.set_yscale("log", nonposy="clip") def _plot_scale_flux(self, energy, flux, energy_power): try: eunit = [_ for _ in flux.unit.bases if _.physical_type == "energy"][0] except IndexError: eunit = energy.unit y = flux * np.power(energy, energy_power) return y.to(flux.unit * eunit ** energy_power)
[docs] def spectral_index(self, energy, epsilon=1e-5): """Compute spectral index at given energy. Parameters ---------- energy : `~astropy.units.Quantity` Energy at which to estimate the index epsilon : float Fractional energy increment to use for determining the spectral index. Returns ------- index : float Estimated spectral index. """ f1 = self(energy) f2 = self(energy * (1 + epsilon)) return np.log(f1 / f2) / np.log(1 + epsilon)
[docs] def inverse(self, value, emin=0.1 * u.TeV, emax=100 * u.TeV): """Return energy for a given function value of the spectral model. Calls the `scipy.optimize.brentq` numerical root finding method. Parameters ---------- value : `~astropy.units.Quantity` Function value of the spectral model. emin : `~astropy.units.Quantity` Lower bracket value in case solution is not unique. emax : `~astropy.units.Quantity` Upper bracket value in case solution is not unique. Returns ------- energy : `~astropy.units.Quantity` Energies at which the model has the given ``value``. """ eunit = "TeV" energies = [] for val in np.atleast_1d(value): def f(x): # scale by 1e12 to achieve better precision energy = u.Quantity(x, eunit, copy=False) y = self(energy).to_value(value.unit) return 1e12 * (y - val.value) energy = brentq(f, emin.to_value(eunit), emax.to_value(eunit)) energies.append(energy) return u.Quantity(energies, eunit, copy=False)
[docs]class ConstantModel(SpectralModel): r"""Constant model .. math:: \phi(E) = k Parameters ---------- const : `~astropy.units.Quantity` :math:`k` """ def __init__(self, const): self.parameters = Parameters([Parameter("const", const)])
[docs] @staticmethod def evaluate(energy, const): """Evaluate the model (static function).""" return np.ones(np.atleast_1d(energy).shape) * const
[docs]class CompoundSpectralModel(SpectralModel): """Represents the algebraic combination of two `~gammapy.spectrum.models.SpectralModel` """ def __init__(self, model1, model2, operator): self.model1 = model1 self.model2 = model2 self.operator = operator # TODO: Think about how to deal with covariance matrix @property def parameters(self): val = self.model1.parameters.parameters + self.model2.parameters.parameters return Parameters(val) @parameters.setter def parameters(self, parameters): idx = len(self.model1.parameters.parameters) self.model1.parameters.parameters = parameters.parameters[:idx] self.model2.parameters.parameters = parameters.parameters[idx:] def __str__(self): ss = self.__class__.__name__ ss += "\n Component 1 : {}".format(self.model1) ss += "\n Component 2 : {}".format(self.model2) ss += "\n Operator : {}".format(self.operator) return ss
[docs] def __call__(self, energy): val1 = self.model1(energy) val2 = self.model2(energy) return self.operator(val1, val2)
[docs] def to_dict(self): retval = dict() retval["model1"] = self.model1.to_dict() retval["model2"] = self.model2.to_dict() retval["operator"] = self.operator
[docs]class PowerLaw(SpectralModel): r"""Spectral power-law model. .. math:: \phi(E) = \phi_0 \cdot \left( \frac{E}{E_0} \right)^{-\Gamma} Parameters ---------- index : `~astropy.units.Quantity` :math:`\Gamma` amplitude : `~astropy.units.Quantity` :math:`\phi_0` reference : `~astropy.units.Quantity` :math:`E_0` Examples -------- This is how to plot the default `PowerLaw` model:: from astropy import units as u from gammapy.spectrum.models import PowerLaw pwl = PowerLaw() pwl.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__(self, index=2.0, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV"): self.parameters = Parameters( [ Parameter("index", index), Parameter("amplitude", amplitude), Parameter("reference", reference, frozen=True), ] )
[docs] @staticmethod def evaluate(energy, index, amplitude, reference): """Evaluate the model (static function).""" return amplitude * np.power((energy / reference), -index)
[docs] def integral(self, emin, emax, **kwargs): r"""Integrate power law analytically. .. math:: F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}\phi(E)dE = \left. \phi_0 \frac{E_0}{-\Gamma + 1} \left( \frac{E}{E_0} \right)^{-\Gamma + 1} \right \vert _{E_{min}}^{E_{max}} Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range """ # kwargs are passed to this function but not used # this is to get a consistent API with SpectralModel.integral() pars = self.parameters if np.isclose(pars["index"].value, 1): e_unit = emin.unit prefactor = pars["amplitude"].quantity * pars["reference"].quantity.to( e_unit ) upper = np.log(emax.to_value(e_unit)) lower = np.log(emin.to_value(e_unit)) else: val = -1 * pars["index"].value + 1 prefactor = pars["amplitude"].quantity * pars["reference"].quantity / val upper = np.power((emax / pars["reference"].quantity), val) lower = np.power((emin / pars["reference"].quantity), val) return prefactor * (upper - lower)
[docs] def integral_error(self, emin, emax, **kwargs): r"""Integrate power law analytically with error propagation. Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. Returns ------- integral, integral_error : tuple of `~astropy.units.Quantity` Tuple of integral flux and integral flux error. """ # kwargs are passed to this function but not used # this is to get a consistent API with SpectralModel.integral() emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.integral(emin, emax, **kwargs).unit upars = self.parameters._ufloats if np.isclose(upars["index"].nominal_value, 1): prefactor = upars["amplitude"] * upars["reference"] upper = np.log(emax.value) lower = np.log(emin.value) else: val = -1 * upars["index"] + 1 prefactor = upars["amplitude"] * upars["reference"] / val upper = np.power((emax.value / upars["reference"]), val) lower = np.power((emin.value / upars["reference"]), val) uarray = prefactor * (upper - lower) return self._parse_uarray(uarray) * unit
[docs] def energy_flux(self, emin, emax): r"""Compute energy flux in given energy range analytically. .. math:: G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}E \phi(E)dE = \left. \phi_0 \frac{E_0^2}{-\Gamma + 2} \left( \frac{E}{E_0} \right)^{-\Gamma + 2} \right \vert _{E_{min}}^{E_{max}} Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. """ pars = self.parameters val = -1 * pars["index"].value + 2 if np.isclose(val, 0): # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2) # for reference temp = pars["amplitude"].quantity * pars["reference"].quantity ** 2 return temp * np.log(emax / emin) else: prefactor = ( pars["amplitude"].quantity * pars["reference"].quantity ** 2 / val ) upper = (emax / pars["reference"].quantity) ** val lower = (emin / pars["reference"].quantity) ** val return prefactor * (upper - lower)
[docs] def energy_flux_error(self, emin, emax, **kwargs): r"""Compute energy flux in given energy range analytically with error propagation. Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. Returns ------- energy_flux, energy_flux_error : tuple of `~astropy.units.Quantity` Tuple of energy flux and energy flux error. """ emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.energy_flux(emin, emax, **kwargs).unit upars = self.parameters._ufloats val = -1 * upars["index"] + 2 if np.isclose(val.nominal_value, 0): # see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2) # for reference temp = upars["amplitude"] * upars["reference"] ** 2 uarray = temp * np.log(emax.value / emin.value) else: prefactor = upars["amplitude"] * upars["reference"] ** 2 / val upper = (emax.value / upars["reference"]) ** val lower = (emin.value / upars["reference"]) ** val uarray = prefactor * (upper - lower) return self._parse_uarray(uarray) * unit
[docs] def inverse(self, value): """Return energy for a given function value of the spectral model. Parameters ---------- value : `~astropy.units.Quantity` Function value of the spectral model. """ p = self.parameters base = value / p["amplitude"].quantity return p["reference"].quantity * np.power(base, -1.0 / p["index"].value)
[docs]class PowerLaw2(SpectralModel): r"""Spectral power-law model with integral as amplitude parameter. See also: https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html .. math:: \phi(E) = F_0 \cdot \frac{\Gamma + 1}{E_{0, max}^{-\Gamma + 1} - E_{0, min}^{-\Gamma + 1}} \cdot E^{-\Gamma} Parameters ---------- index : `~astropy.units.Quantity` Spectral index :math:`\Gamma` amplitude : `~astropy.units.Quantity` Integral flux :math:`F_0`. emin : `~astropy.units.Quantity` Lower energy limit :math:`E_{0, min}`. emax : `~astropy.units.Quantity` Upper energy limit :math:`E_{0, max}`. Examples -------- This is how to plot the default `PowerLaw2` model:: from astropy import units as u from gammapy.spectrum.models import PowerLaw2 pwl2 = PowerLaw2() pwl2.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__( self, amplitude="1e-12 cm-2 s-1", index=2, emin="0.1 TeV", emax="100 TeV" ): self.parameters = Parameters( [ Parameter("amplitude", amplitude), Parameter("index", index), Parameter("emin", emin, frozen=True), Parameter("emax", emax, frozen=True), ] )
[docs] @staticmethod def evaluate(energy, amplitude, index, emin, emax): """Evaluate the model (static function).""" top = -index + 1 # to get the energies dimensionless we use a modified formula bottom = emax - emin * (emin / emax) ** (-index) return amplitude * (top / bottom) * np.power(energy / emax, -index)
[docs] def integral(self, emin, emax, **kwargs): r"""Integrate power law analytically. .. math:: F(E_{min}, E_{max}) = F_0 \cdot \frac{E_{max}^{\Gamma + 1} \ - E_{min}^{\Gamma + 1}}{E_{0, max}^{\Gamma + 1} \ - E_{0, min}^{\Gamma + 1}} Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. """ pars = self.parameters temp1 = np.power(emax, -pars["index"].value + 1) temp2 = np.power(emin, -pars["index"].value + 1) top = temp1 - temp2 temp1 = np.power(pars["emax"].quantity, -pars["index"].value + 1) temp2 = np.power(pars["emin"].quantity, -pars["index"].value + 1) bottom = temp1 - temp2 return pars["amplitude"].quantity * top / bottom
[docs] def integral_error(self, emin, emax, **kwargs): r"""Integrate power law analytically with error propagation. Parameters ---------- emin, emax : `~astropy.units.Quantity` Lower and upper bound of integration range. Returns ------- integral, integral_error : tuple of `~astropy.units.Quantity` Tuple of integral flux and integral flux error. """ emin = self._convert_energy(emin) emax = self._convert_energy(emax) unit = self.integral(emin, emax, **kwargs).unit upars = self.parameters._ufloats temp1 = np.power(emax.value, -upars["index"] + 1) temp2 = np.power(emin.value, -upars["index"] + 1) top = temp1 - temp2 temp1 = np.power(upars["emax"], -upars["index"] + 1) temp2 = np.power(upars["emin"], -upars["index"] + 1) bottom = temp1 - temp2 uarray = upars["amplitude"] * top / bottom return self._parse_uarray(uarray) * unit
[docs] def inverse(self, value): """Return energy for a given function value of the spectral model. Parameters ---------- value : `~astropy.units.Quantity` Function value of the spectral model. """ p = self.parameters amplitude, index, emin, emax = ( p["amplitude"].quantity, p["index"].value, p["emin"].quantity, p["emax"].quantity, ) # to get the energies dimensionless we use a modified formula top = -index + 1 bottom = emax - emin * (emin / emax) ** (-index) term = (bottom / top) * (value / amplitude) return np.power(term.to_value(""), -1.0 / index) * emax
[docs]class ExponentialCutoffPowerLaw(SpectralModel): r"""Spectral exponential cutoff power-law model. .. math:: \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma} \exp(-\lambda E) Parameters ---------- index : `~astropy.units.Quantity` :math:`\Gamma` amplitude : `~astropy.units.Quantity` :math:`\phi_0` reference : `~astropy.units.Quantity` :math:`E_0` lambda_ : `~astropy.units.Quantity` :math:`\lambda` Examples -------- This is how to plot the default `ExponentialCutoffPowerLaw` model:: from astropy import units as u from gammapy.spectrum.models import ExponentialCutoffPowerLaw ecpl = ExponentialCutoffPowerLaw() ecpl.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__( self, index=1.5, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV", lambda_="0.1 TeV-1", ): self.parameters = Parameters( [ Parameter("index", index), Parameter("amplitude", amplitude), Parameter("reference", reference, frozen=True), Parameter("lambda_", lambda_), ] )
[docs] @staticmethod def evaluate(energy, index, amplitude, reference, lambda_): """Evaluate the model (static function).""" pwl = amplitude * (energy / reference) ** (-index) try: cutoff = np.exp(-energy * lambda_) except AttributeError: from uncertainties.unumpy import exp cutoff = exp(-energy * lambda_) return pwl * cutoff
@property def e_peak(self): r"""Spectral energy distribution peak energy (`~astropy.utils.Quantity`). This is the peak in E^2 x dN/dE and is given by: .. math:: E_{Peak} = (2 - \Gamma) / \lambda """ p = self.parameters reference = p["reference"].quantity index = p["index"].quantity lambda_ = p["lambda_"].quantity if index >= 2: return np.nan * reference.unit else: return (2 - index) / lambda_
[docs]class ExponentialCutoffPowerLaw3FGL(SpectralModel): r"""Spectral exponential cutoff power-law model used for 3FGL. Note that the parametrization is different from `ExponentialCutoffPowerLaw`: .. math:: \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma} \exp \left( \frac{E_0 - E}{E_{C}} \right) Parameters ---------- index : `~astropy.units.Quantity` :math:`\Gamma` amplitude : `~astropy.units.Quantity` :math:`\phi_0` reference : `~astropy.units.Quantity` :math:`E_0` ecut : `~astropy.units.Quantity` :math:`E_{C}` Examples -------- This is how to plot the default `ExponentialCutoffPowerLaw3FGL` model:: from astropy import units as u from gammapy.spectrum.models import ExponentialCutoffPowerLaw3FGL ecpl_3fgl = ExponentialCutoffPowerLaw3FGL() ecpl_3fgl.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__( self, index=1.5, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV", ecut="10 TeV", ): self.parameters = Parameters( [ Parameter("index", index), Parameter("amplitude", amplitude), Parameter("reference", reference, frozen=True), Parameter("ecut", ecut), ] )
[docs] @staticmethod def evaluate(energy, index, amplitude, reference, ecut): """Evaluate the model (static function).""" pwl = amplitude * (energy / reference) ** (-index) try: cutoff = np.exp((reference - energy) / ecut) except AttributeError: from uncertainties.unumpy import exp cutoff = exp((reference - energy) / ecut) return pwl * cutoff
[docs]class PLSuperExpCutoff3FGL(SpectralModel): r"""Spectral super exponential cutoff power-law model used for 3FGL. .. math:: \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma_1} \exp \left( \left(\frac{E_0}{E_{C}} \right)^{\Gamma_2} - \left(\frac{E}{E_{C}} \right)^{\Gamma_2} \right) Parameters ---------- index_1 : `~astropy.units.Quantity` :math:`\Gamma_1` index_2 : `~astropy.units.Quantity` :math:`\Gamma_2` amplitude : `~astropy.units.Quantity` :math:`\phi_0` reference : `~astropy.units.Quantity` :math:`E_0` ecut : `~astropy.units.Quantity` :math:`E_{C}` Examples -------- This is how to plot the default `PLSuperExpCutoff3FGL` model:: from astropy import units as u from gammapy.spectrum.models import PLSuperExpCutoff3FGL secpl_3fgl = PLSuperExpCutoff3FGL() secpl_3fgl.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__( self, index_1=1.5, index_2=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV", ecut="10 TeV", ): self.parameters = Parameters( [ Parameter("index_1", index_1), Parameter("index_2", index_2), Parameter("amplitude", amplitude), Parameter("reference", reference, frozen=True), Parameter("ecut", ecut), ] )
[docs] @staticmethod def evaluate(energy, amplitude, reference, ecut, index_1, index_2): """Evaluate the model (static function).""" pwl = amplitude * (energy / reference) ** (-index_1) try: cutoff = np.exp((reference / ecut) ** index_2 - (energy / ecut) ** index_2) except AttributeError: from uncertainties.unumpy import exp cutoff = exp((reference / ecut) ** index_2 - (energy / ecut) ** index_2) return pwl * cutoff
[docs]class LogParabola(SpectralModel): r"""Spectral log parabola model. .. math:: \phi(E) = \phi_0 \left( \frac{E}{E_0} \right) ^ { - \alpha - \beta \log{ \left( \frac{E}{E_0} \right) } } Note that :math:`log` refers to the natural logarithm. This is consistent with the `Fermi Science Tools <https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html>`_ and `ctools <http://cta.irap.omp.eu/ctools-devel/users/user_manual/getting_started/models.html#log-parabola>`_. The `Sherpa <http://cxc.harvard.edu/sherpa/ahelp/logparabola.html_ package>`_ package, however, uses :math:`log_{10}`. If you have parametrization based on :math:`log_{10}` you can use the :func:`~gammapy.spectrum.models.LogParabola.from_log10` method. Parameters ---------- amplitude : `~astropy.units.Quantity` :math:`\phi_0` reference : `~astropy.units.Quantity` :math:`E_0` alpha : `~astropy.units.Quantity` :math:`\alpha` beta : `~astropy.units.Quantity` :math:`\beta` Examples -------- This is how to plot the default `LogParabola` model:: from astropy import units as u from gammapy.spectrum.models import LogParabola log_parabola = LogParabola() log_parabola.plot(energy_range=[0.1, 100] * u.TeV) plt.show() """ def __init__( self, amplitude="1e-12 cm-2 s-1 TeV-1", reference="10 TeV", alpha=2, beta=1 ): self.parameters = Parameters( [ Parameter("amplitude", amplitude), Parameter("reference", reference, frozen=True), Parameter("alpha", alpha), Parameter("beta", beta), ] )
[docs] @classmethod def from_log10(cls, amplitude, reference, alpha, beta): """Construct LogParabola from :math:`log_{10}` parametrization""" beta_ = beta / np.log(10) return cls(amplitude=amplitude, reference=reference, alpha=alpha, beta=beta_)
[docs] @staticmethod def evaluate(energy, amplitude, reference, alpha, beta): """Evaluate the model (static function).""" try: xx = (energy / reference).to("") exponent = -alpha - beta * np.log(xx) except AttributeError: from uncertainties.unumpy import log xx = energy / reference exponent = -alpha - beta * log(xx) return amplitude * np.power(xx, exponent)
@property def e_peak(self): r"""Spectral energy distribution peak energy (`~astropy.units.Quantity`). This is the peak in E^2 x dN/dE and is given by: .. math:: E_{Peak} = E_{0} \exp{ (2 - \alpha) / (2 * \beta)} """ p = self.parameters reference = p["reference"].quantity alpha = p["alpha"].quantity beta = p["beta"].quantity return reference * np.exp((2 - alpha) / (2 * beta))
[docs]class TableModel(SpectralModel): """A model generated from a table of energy and value arrays. the units returned will be the units of the values array provided at initialization. The model will return values interpolated in log-space, returning 0 for energies outside of the limits of the provided energy array. Class implementation follows closely what has been done in `naima.models.TableModel` Parameters ---------- energy : `~astropy.units.Quantity` array Array of energies at which the model values are given values : array Array with the values of the model at energies ``energy``. norm : float Model scale that is multiplied to the supplied arrays. Defaults to 1. values_scale : {'log', 'lin', 'sqrt'} Interpolation scaling applied to values. If the values vary over many magnitudes a 'log' scaling is recommended. interp_kwargs : dict Interpolation keyword arguments pass to `scipy.interpolate.interp1d`. By default all values outside the interpolation range are set to zero. If you want to apply linear extrapolation you can pass `interp_kwargs={'fill_value': 'extrapolate', 'kind': 'linear'}` meta : dict, optional Meta information, meta['filename'] will be used for serialization """ def __init__( self, energy, values, norm=1, values_scale="log", interp_kwargs=None, meta=None ): self.parameters = Parameters([Parameter("norm", norm, unit="")]) self.energy = energy self.values = values self.meta = dict() if meta is None else meta interp_kwargs = interp_kwargs or {} interp_kwargs.setdefault("values_scale", "log") self._evaluate = ScaledRegularGridInterpolator( points=(np.log(energy.value),), values=values, **interp_kwargs )
[docs] @classmethod def read_xspec_model(cls, filename, param, **kwargs): """Read XSPEC table model The input is a table containing absorbed values from a XSPEC model as a function of energy. TODO: Format of the file should be described and discussed in https://gamma-astro-data-formats.readthedocs.io/en/latest/index.html Parameters ---------- filename : str File containing the XSPEC model param : float Model parameter value Examples -------- Fill table from an EBL model (Franceschini, 2008) >>> from gammapy.spectrum.models import TableModel >>> filename = '$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz' >>> table_model = TableModel.read_xspec_model(filename=filename, param=0.3) """ filename = str(make_path(filename)) # Check if parameter value is in range table_param = Table.read(filename, hdu="PARAMETERS") param_min = table_param["MINIMUM"] param_max = table_param["MAXIMUM"] if param < param_min or param > param_max: err = "Parameter out of range, param={}, param_min={}, param_max={}".format( param, param_min, param_max ) raise ValueError(err) # Get energy values table_energy = Table.read(filename, hdu="ENERGIES") energy_lo = table_energy["ENERG_LO"] energy_hi = table_energy["ENERG_HI"] # Hack while format is not fixed, energy values are in keV energy_bounds = EnergyBounds.from_lower_and_upper_bounds( lower=energy_lo, upper=energy_hi, unit=u.keV ) energy = energy_bounds.log_centers # Get spectrum values (no interpolation, take closest value for param) table_spectra = Table.read(filename, hdu="SPECTRA") idx = np.abs(table_spectra["PARAMVAL"] - param).argmin() values = u.Quantity(table_spectra[idx][1], "", copy=False) # no dimension kwargs.setdefault("values_scale", "lin") return cls(energy=energy, values=values, **kwargs)
[docs] @classmethod def read_fermi_isotropic_model(cls, filename, **kwargs): """Read Fermi isotropic diffuse model see `LAT Background models <https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html>`_ Parameters ---------- filename : str filename """ filename = str(make_path(filename)) vals = np.loadtxt(filename) energy = u.Quantity(vals[:, 0], "MeV", copy=False) values = u.Quantity(vals[:, 1], "MeV-1 s-1 cm-2 sr-1", copy=False) return cls(energy=energy, values=values, **kwargs)
[docs] def evaluate(self, energy, norm): """Evaluate the model (static function).""" x = np.log(energy.to_value(self.energy.unit)) values = self._evaluate((x,), clip=True) return norm * values
class ScaleModel(SpectralModel): """Wrapper to scale another spectral model by a norm factor. Parameters ---------- model : `SpectralModel` Spectral model to wrap. norm : float Multiplicative norm factor for the model value. """ def __init__(self, model, norm=1): self.parameters = Parameters([Parameter("norm", norm, unit="")]) self.model = model def evaluate(self, energy, norm): return norm * self.model(energy)
[docs]class Absorption(object): r"""Gamma-ray absorption models. Parameters ---------- energy_lo, energy_hi : `~astropy.units.Quantity` Lower and upper bin edges of energy axis param_lo, param_hi : `~astropy.units.Quantity` Lower and upper bin edges of parameter axis data : `~astropy.units.Quantity` Model value Examples -------- Create and plot EBL absorption models for a redshift of 0.5: .. plot:: :include-source: import matplotlib.pyplot as plt import astropy.units as u from gammapy.spectrum.models import Absorption # Load tables for z=0.5 redshift = 0.5 dominguez = Absorption.read_builtin('dominguez').table_model(redshift) franceschini = Absorption.read_builtin('franceschini').table_model(redshift) finke = Absorption.read_builtin('finke').table_model(redshift) # start customised plot energy_range = [0.08, 3] * u.TeV ax = plt.gca() opts = dict(energy_range=energy_range, energy_unit='TeV', ax=ax, flux_unit='') franceschini.plot(label='Franceschini 2008', **opts) finke.plot(label='Finke 2010', **opts) dominguez.plot(label='Dominguez 2011', **opts) # tune plot ax.set_ylabel(r'Absorption coefficient [$\exp{(-\tau(E))}$]') ax.set_xlim(energy_range.value) # we get ride of units ax.set_ylim([1.e-4, 2.]) ax.set_yscale('log') ax.set_title('EBL models (z=' + str(redshift) + ')') plt.grid(which='both') plt.legend(loc='best') # legend # show plot plt.show() """ def __init__(self, energy_lo, energy_hi, param_lo, param_hi, data): axes = [ BinnedDataAxis( param_lo, param_hi, interpolation_mode="linear", name="parameter" ), BinnedDataAxis( energy_lo, energy_hi, interpolation_mode="log", name="energy" ), ] self.data = NDDataArray(axes=axes, data=data) self.data.default_interp_kwargs["fill_value"] = None
[docs] @classmethod def read(cls, filename): """Build object from an XSPEC model. Todo: Format of XSPEC binary files should be referenced at https://gamma-astro-data-formats.readthedocs.io/en/latest/ Parameters ---------- filename : str File containing the model. """ # Create EBL data array filename = str(make_path(filename)) table_param = Table.read(filename, hdu="PARAMETERS") par_min = table_param["MINIMUM"] par_max = table_param["MAXIMUM"] par_array = table_param[0]["VALUE"] par_delta = np.diff(par_array) * 0.5 param_lo, param_hi = par_array, par_array # initialisation param_lo[0] = par_min - par_delta[0] param_lo[1:] -= par_delta param_hi[:-1] += par_delta param_hi[-1] = par_max # Get energy values table_energy = Table.read(filename, hdu="ENERGIES") energy_lo = u.Quantity( table_energy["ENERG_LO"], "keV", copy=False ) # unit not stored in file energy_hi = u.Quantity( table_energy["ENERG_HI"], "keV", copy=False ) # unit not stored in file # Get spectrum values table_spectra = Table.read(filename, hdu="SPECTRA") data = table_spectra["INTPSPEC"].data return cls( energy_lo=energy_lo, energy_hi=energy_hi, param_lo=param_lo, param_hi=param_hi, data=data, )
[docs] @classmethod def read_builtin(cls, name): """Read one of the built-in absorption models. Parameters ---------- name : {'franceschini', 'dominguez', 'finke'} name of one of the available model in gammapy-extra References ---------- .. [1] Franceschini et al., "Extragalactic optical-infrared background radiation, its time evolution and the cosmic photon-photon opacity", `Link <http://adsabs.harvard.edu/abs/2008A%26A...487..837F>`__ .. [2] Dominguez et al., " Extragalactic background light inferred from AEGIS galaxy-SED-type fractions" `Link <http://adsabs.harvard.edu/abs/2011MNRAS.410.2556D>`__ .. [3] Finke et al., "Modeling the Extragalactic Background Light from Stars and Dust" `Link <http://adsabs.harvard.edu/abs/2010ApJ...712..238F>`__ """ models = dict() models["franceschini"] = "$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz" models["dominguez"] = "$GAMMAPY_DATA/ebl/ebl_dominguez11.fits.gz" models["finke"] = "$GAMMAPY_DATA/ebl/frd_abs.fits.gz" return cls.read(models[name])
[docs] def table_model(self, parameter, unit="TeV"): """Table model for a given parameter (`~gammapy.spectrum.models.TableModel`). Parameters ---------- parameter : float Parameter value. unit : str, (optional) desired value for energy axis """ energy_axis = self.data.axes[1] energy = (energy_axis.log_center()).to(unit) values = self.evaluate(energy=energy, parameter=parameter) return TableModel(energy=energy, values=values, values_scale="lin")
[docs] def evaluate(self, energy, parameter): """Evaluate model for energy and parameter value.""" return self.data.evaluate(energy=energy, parameter=parameter)
[docs]class AbsorbedSpectralModel(SpectralModel): """Spectral model with EBL absorption. Parameters ---------- spectral_model : `~gammapy.spectrum.models.SpectralModel` Spectral model. absorption : `~gammapy.spectrum.models.Absorption` Absorption model. parameter : float parameter value for absorption model parameter_name : str, optional parameter name """ def __init__( self, spectral_model, absorption, parameter, parameter_name="redshift" ): self.spectral_model = spectral_model self.absorption = absorption self.parameter = parameter self.parameter_name = parameter_name # initialise list parameters from spectral model param_list = [] for param in spectral_model.parameters.parameters: param_list.append(param) # Add parameter to the list min_ = self.absorption.data.axes[0].lo[0] max_ = self.absorption.data.axes[0].lo[-1] par = Parameter(parameter_name, parameter, min=min_, max=max_, frozen=True) param_list.append(par) self.parameters = Parameters(param_list)
[docs] def evaluate(self, energy, **kwargs): """Evaluate the model at a given energy.""" # assign redshift value and remove it from dictionnary # since it does not belong to the spectral model parameter = kwargs[self.parameter_name] del kwargs[self.parameter_name] flux = self.spectral_model.evaluate(energy=energy, **kwargs) absorption = self.absorption.evaluate(energy=energy, parameter=parameter) return flux * absorption