Source code for gammapy.irf.effective_area

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
import astropy.units as u
from astropy.visualization import quantity_support
import matplotlib.pyplot as plt
from gammapy.maps import MapAxes, MapAxis
from .core import IRF

__all__ = ["EffectiveAreaTable2D"]


[docs]class EffectiveAreaTable2D(IRF): """2D effective area table. Data format specification: :ref:`gadf:aeff_2d` Parameters ---------- energy_axis_true : `MapAxis` True energy axis offset_axis : `MapAxis` Field of view offset axis. data : `~astropy.units.Quantity` Effective area meta : dict Meta data Examples -------- Here's an example you can use to learn about this class: >>> from gammapy.irf import EffectiveAreaTable2D >>> filename = '$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits' >>> aeff = EffectiveAreaTable2D.read(filename, hdu='EFFECTIVE AREA') >>> print(aeff) EffectiveAreaTable2D -------------------- <BLANKLINE> axes : ['energy_true', 'offset'] shape : (42, 6) ndim : 2 unit : m2 dtype : >f4 <BLANKLINE> Here's another one, created from scratch, without reading a file: >>> from gammapy.irf import EffectiveAreaTable2D >>> from gammapy.maps import MapAxis >>> energy_axis_true = MapAxis.from_energy_bounds( "0.1 TeV", "100 TeV", nbin=30, name="energy_true" ) >>> offset_axis = MapAxis.from_bounds(0, 5, nbin=4, name="offset") >>> aeff = EffectiveAreaTable2D(axes=[energy_axis_true, offset_axis], data=1e10, unit="cm2") >>> print(aeff) EffectiveAreaTable2D -------------------- <BLANKLINE> axes : ['energy_true', 'offset'] shape : (30, 4) ndim : 2 unit : cm2 dtype : float64 <BLANKLINE> """ tag = "aeff_2d" required_axes = ["energy_true", "offset"] default_unit = u.m**2
[docs] def plot_energy_dependence(self, ax=None, offset=None, **kwargs): """Plot effective area versus energy for a given offset. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Axis offset : `~astropy.coordinates.Angle` Offset kwargs : dict Forwarded tp plt.plot() Returns ------- ax : `~matplotlib.axes.Axes` Axis """ ax = plt.gca() if ax is None else ax if offset is None: off_min, off_max = self.axes["offset"].bounds offset = np.linspace(off_min, off_max, 4) energy_axis = self.axes["energy_true"] for off in offset: area = self.evaluate(offset=off, energy_true=energy_axis.center) label = kwargs.pop("label", f"offset = {off:.1f}") with quantity_support(): ax.plot(energy_axis.center, area, label=label, **kwargs) energy_axis.format_plot_xaxis(ax=ax) ax.set_ylabel(f"Effective Area [{ax.yaxis.units}]") ax.legend() return ax
[docs] def plot_offset_dependence(self, ax=None, energy=None, **kwargs): """Plot effective area versus offset for a given energy. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Axis energy : `~astropy.units.Quantity` Energy **kwargs : dict Keyword argument passed to `~matplotlib.pyplot.plot` Returns ------- ax : `~matplotlib.axes.Axes` Axis """ ax = plt.gca() if ax is None else ax if energy is None: energy_axis = self.axes["energy_true"] e_min, e_max = energy_axis.center[[0, -1]] energy = np.geomspace(e_min, e_max, 4) offset_axis = self.axes["offset"] for ee in energy: area = self.evaluate(offset=offset_axis.center, energy_true=ee) area /= np.nanmax(area) if np.isnan(area).all(): continue label = f"energy = {ee:.1f}" with quantity_support(): ax.plot(offset_axis.center, area, label=label, **kwargs) offset_axis.format_plot_xaxis(ax=ax) ax.set_ylim(0, 1.1) ax.set_ylabel("Relative Effective Area") ax.legend(loc="best") return ax
[docs] def plot(self, ax=None, add_cbar=True, **kwargs): """Plot effective area image.""" ax = plt.gca() if ax is None else ax energy = self.axes["energy_true"] offset = self.axes["offset"] aeff = self.evaluate( offset=offset.center, energy_true=energy.center[:, np.newaxis] ) vmin, vmax = np.nanmin(aeff.value), np.nanmax(aeff.value) kwargs.setdefault("cmap", "GnBu") kwargs.setdefault("edgecolors", "face") kwargs.setdefault("vmin", vmin) kwargs.setdefault("vmax", vmax) with quantity_support(): caxes = ax.pcolormesh(energy.edges, offset.edges, aeff.value.T, **kwargs) energy.format_plot_xaxis(ax=ax) offset.format_plot_yaxis(ax=ax) if add_cbar: label = f"Effective Area [{aeff.unit}]" ax.figure.colorbar(caxes, ax=ax, label=label) return ax
[docs] def peek(self, figsize=(15, 5)): """Quick-look summary plots. Parameters ---------- figsize : tuple Size of the figure. """ fig, axes = plt.subplots(nrows=1, ncols=3, figsize=figsize) self.plot(ax=axes[2]) self.plot_energy_dependence(ax=axes[0]) self.plot_offset_dependence(ax=axes[1]) plt.tight_layout()
[docs] @classmethod def from_parametrization(cls, energy_axis_true=None, instrument="HESS"): r"""Create parametrized effective area. Parametrizations of the effective areas of different Cherenkov telescopes taken from Appendix B of Abramowski et al. (2010), see https://ui.adsabs.harvard.edu/abs/2010MNRAS.402.1342A . .. math:: A_{eff}(E) = g_1 \left(\frac{E}{\mathrm{MeV}}\right)^{-g_2}\exp{\left(-\frac{g_3}{E}\right)} # noqa: E501 This method does not model the offset dependence of the effective area, but just assumes that it is constant. Parameters ---------- energy_axis_true : `MapAxis` Energy binning, analytic function is evaluated at log centers instrument : {'HESS', 'HESS2', 'CTA'} Instrument name Returns ------- aeff : `EffectiveAreaTable2D` Effective area table """ # Put the parameters g in a dictionary. # Units: g1 (cm^2), g2 (), g3 (MeV) pars = { "HESS": [6.85e9, 0.0891, 5e5], "HESS2": [2.05e9, 0.0891, 1e5], "CTA": [1.71e11, 0.0891, 1e5], } if instrument not in pars.keys(): ss = f"Unknown instrument: {instrument}\n" ss += f"Valid instruments: {list(pars.keys())}" raise ValueError(ss) if energy_axis_true is None: energy_axis_true = MapAxis.from_energy_bounds( "2 GeV", "200 TeV", nbin=20, per_decade=True, name="energy_true" ) g1, g2, g3 = pars[instrument] offset_axis = MapAxis.from_edges([0.0, 5.0] * u.deg, name="offset") axes = MapAxes([energy_axis_true, offset_axis]) coords = axes.get_coord() energy, offset = coords["energy_true"].to_value("MeV"), coords["offset"] data = np.ones_like(offset.value) * g1 * energy ** (-g2) * np.exp(-g3 / energy) # TODO: fake offset dependence? meta = {"TELESCOP": instrument} return cls(axes=axes, data=data, unit="cm2", meta=meta)