# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
import scipy.special
from astropy import units as u
from astropy.coordinates import Angle, SkyCoord
from astropy.visualization import quantity_support
import matplotlib.pyplot as plt
from matplotlib.colors import PowerNorm
from gammapy.maps import MapAxes, MapAxis, RegionGeom
from ..core import IRF
__all__ = ["EnergyDispersion2D"]
[docs]class EnergyDispersion2D(IRF):
"""Offset-dependent energy dispersion matrix.
Data format specification: :ref:`gadf:edisp_2d`
Parameters
----------
energy_axis_true : `MapAxis`
True energy axis
migra_axis : `MapAxis`
Energy migration axis
offset_axis : `MapAxis`
Field of view offset axis
data : `~numpy.ndarray`
Energy dispersion probability density
Examples
--------
Read energy dispersion IRF from disk:
>>> from gammapy.maps import MapAxis
>>> from gammapy.irf import EnergyDispersion2D
>>> filename = '$GAMMAPY_DATA/hess-dl3-dr1/data/hess_dl3_dr1_obs_id_020136.fits.gz'
>>> edisp2d = EnergyDispersion2D.read(filename, hdu="EDISP")
Create energy dispersion matrix (`~gammapy.irf.EnergyDispersion`)
for a given field of view offset and energy binning:
>>> energy = MapAxis.from_bounds(0.1, 20, nbin=60, unit="TeV", interp="log").edges
>>> edisp = edisp2d.to_edisp_kernel(offset='1.2 deg', energy=energy, energy_true=energy)
See Also
--------
EnergyDispersion
"""
tag = "edisp_2d"
required_axes = ["energy_true", "migra", "offset"]
default_unit = u.one
@property
def _default_offset(self):
if self.axes["offset"].nbin == 1:
default_offset = self.axes["offset"].center
else:
default_offset = [1.0] * u.deg
return default_offset
def _mask_out_bounds(self, invalid):
return (
invalid[self.axes.index("energy_true")] & invalid[self.axes.index("migra")]
) | invalid[self.axes.index("offset")]
[docs] @classmethod
def from_gauss(
cls, energy_axis_true, migra_axis, offset_axis, bias, sigma, pdf_threshold=1e-6
):
"""Create Gaussian energy dispersion matrix (`EnergyDispersion2D`).
The output matrix will be Gaussian in (energy_true / energy).
The ``bias`` and ``sigma`` should be either floats or arrays of same dimension than
``energy_true``. ``bias`` refers to the mean value of the ``migra``
distribution minus one, i.e. ``bias=0`` means no bias.
Note that, the output matrix is flat in offset.
Parameters
----------
energy_axis_true : `MapAxis`
True energy axis
migra_axis : `~astropy.units.Quantity`
Migra axis
offset_axis : `~astropy.units.Quantity`
Bin edges of offset
bias : float or `~numpy.ndarray`
Center of Gaussian energy dispersion, bias
sigma : float or `~numpy.ndarray`
RMS width of Gaussian energy dispersion, resolution.
pdf_threshold : float, optional
Zero suppression threshold
"""
axes = MapAxes([energy_axis_true, migra_axis, offset_axis])
coords = axes.get_coord(mode="edges", axis_name="migra")
migra_min = coords["migra"][:, :-1, :]
migra_max = coords["migra"][:, 1:, :]
# Analytical formula for integral of Gaussian
s = np.sqrt(2) * sigma
t1 = (migra_max - 1 - bias) / s
t2 = (migra_min - 1 - bias) / s
pdf = (scipy.special.erf(t1) - scipy.special.erf(t2)) / 2
pdf = pdf / (migra_max - migra_min)
# no offset dependence
data = pdf.T * np.ones(axes.shape)
data[data < pdf_threshold] = 0
return cls(
axes=axes,
data=data.value,
)
[docs] def to_edisp_kernel(self, offset, energy_true=None, energy=None):
"""Detector response R(Delta E_reco, Delta E_true)
Probability to reconstruct an energy in a given true energy band
in a given reconstructed energy band
Parameters
----------
offset : `~astropy.coordinates.Angle`
Offset
energy_true : `~astropy.units.Quantity`, None
True energy axis
energy : `~astropy.units.Quantity`
Reconstructed energy axis
Returns
-------
edisp : `~gammapy.irf.EDispKernel`
Energy dispersion matrix
"""
from gammapy.makers.utils import make_edisp_kernel_map
offset = Angle(offset)
# TODO: expect directly MapAxis here?
if energy is None:
energy_axis = self.axes["energy_true"].copy(name="energy")
else:
energy_axis = MapAxis.from_energy_edges(energy)
if energy_true is None:
energy_axis_true = self.axes["energy_true"]
else:
energy_axis_true = MapAxis.from_energy_edges(
energy_true,
name="energy_true",
)
pointing = SkyCoord("0d", "0d")
center = pointing.directional_offset_by(
position_angle=0 * u.deg, separation=offset
)
geom = RegionGeom.create(region=center, axes=[energy_axis, energy_axis_true])
edisp = make_edisp_kernel_map(geom=geom, edisp=self, pointing=pointing)
return edisp.get_edisp_kernel()
[docs] def normalize(self):
"""Normalise energy dispersion"""
super().normalize(axis_name="migra")
[docs] def plot_migration(self, ax=None, offset=None, energy_true=None, **kwargs):
"""Plot energy dispersion for given offset and true energy.
Parameters
----------
ax : `~matplotlib.axes.Axes`, optional
Axis
offset : `~astropy.coordinates.Angle`, optional
Offset
energy_true : `~astropy.units.Quantity`, optional
True energy
**kwargs : dict
Keyword arguments forwarded to `~matplotlib.pyplot.plot`
Returns
-------
ax : `~matplotlib.axes.Axes`
Axis
"""
ax = plt.gca() if ax is None else ax
if offset is None:
offset = self._default_offset
else:
offset = np.atleast_1d(Angle(offset))
if energy_true is None:
energy_true = u.Quantity([0.1, 1, 10], "TeV")
else:
energy_true = np.atleast_1d(u.Quantity(energy_true))
migra = self.axes["migra"]
with quantity_support():
for ener in energy_true:
for off in offset:
disp = self.evaluate(
offset=off, energy_true=ener, migra=migra.center
)
label = f"offset = {off:.1f}\nenergy = {ener:.1f}"
ax.plot(migra.center, disp, label=label, **kwargs)
migra.format_plot_xaxis(ax=ax)
ax.set_ylabel("Probability density")
ax.legend(loc="upper left")
return ax
[docs] def plot_bias(self, ax=None, offset=None, add_cbar=False, **kwargs):
"""Plot migration as a function of true energy for a given offset.
Parameters
----------
ax : `~matplotlib.axes.Axes`, optional
Axis
offset : `~astropy.coordinates.Angle`, optional
Offset
add_cbar : bool
Add a colorbar to the plot.
kwargs : dict
Keyword arguments passed to `~matplotlib.pyplot.pcolormesh`.
Returns
-------
ax : `~matplotlib.axes.Axes`
Axis
"""
kwargs.setdefault("cmap", "GnBu")
kwargs.setdefault("norm", PowerNorm(gamma=0.5))
ax = plt.gca() if ax is None else ax
if offset is None:
offset = self._default_offset
energy_true = self.axes["energy_true"]
migra = self.axes["migra"]
z = self.evaluate(
offset=offset,
energy_true=energy_true.center.reshape(1, -1, 1),
migra=migra.center.reshape(1, 1, -1),
).value[0]
with quantity_support():
caxes = ax.pcolormesh(energy_true.edges, migra.edges, z.T, **kwargs)
energy_true.format_plot_xaxis(ax=ax)
migra.format_plot_yaxis(ax=ax)
if add_cbar:
label = "Probability density [A.U]."
ax.figure.colorbar(caxes, ax=ax, label=label)
return ax
[docs] def peek(self, figsize=(15, 5)):
"""Quick-look summary plots.
Parameters
----------
figsize : (float, float)
Size of the resulting plot
"""
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=figsize)
self.plot_bias(ax=axes[0])
self.plot_migration(ax=axes[1])
edisp = self.to_edisp_kernel(offset=self._default_offset[0])
edisp.plot_matrix(ax=axes[2])
plt.tight_layout()