Source code for gammapy.irf.background

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
import astropy.units as u
from astropy.io import fits
from astropy.table import Table
from gammapy.maps import MapAxis
from gammapy.maps.utils import edges_from_lo_hi
from gammapy.utils.integrate import trapz_loglog
from gammapy.utils.nddata import NDDataArray
from gammapy.utils.scripts import make_path

__all__ = ["Background3D", "Background2D"]


[docs]class Background3D: """Background 3D. Data format specification: :ref:`gadf:bkg_3d` Parameters ---------- energy_lo, energy_hi : `~astropy.units.Quantity` Energy binning fov_lon_lo, fov_lon_hi : `~astropy.units.Quantity` FOV coordinate X-axis binning. fov_lat_lo, fov_lat_hi : `~astropy.units.Quantity` FOV coordinate Y-axis binning. data : `~astropy.units.Quantity` Background rate (usually: ``s^-1 MeV^-1 sr^-1``) Examples -------- Here's an example you can use to learn about this class: >>> from gammapy.irf import Background3D >>> filename = '$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits' >>> bkg_3d = Background3D.read(filename, hdu='BACKGROUND') >>> print(bkg_3d) Background3D NDDataArray summary info energy : size = 21, min = 0.016 TeV, max = 158.489 TeV fov_lon : size = 36, min = -5.833 deg, max = 5.833 deg fov_lat : size = 36, min = -5.833 deg, max = 5.833 deg Data : size = 27216, min = 0.000 1 / (MeV s sr), max = 0.421 1 / (MeV s sr) """ default_interp_kwargs = dict( bounds_error=False, fill_value=None, values_scale="log" ) """Default Interpolation kwargs for `~gammapy.utils.nddata.NDDataArray`. Extrapolate.""" def __init__( self, energy_lo, energy_hi, fov_lon_lo, fov_lon_hi, fov_lat_lo, fov_lat_hi, data, meta=None, interp_kwargs=None, ): if interp_kwargs is None: interp_kwargs = self.default_interp_kwargs e_edges = edges_from_lo_hi(energy_lo, energy_hi) energy_axis = MapAxis.from_edges(e_edges, interp="log", name="energy") fov_lon_edges = edges_from_lo_hi(fov_lon_lo, fov_lon_hi) fov_lon_axis = MapAxis.from_edges(fov_lon_edges, interp="lin", name="fov_lon") fov_lat_edges = edges_from_lo_hi(fov_lat_lo, fov_lat_hi) fov_lat_axis = MapAxis.from_edges(fov_lat_edges, interp="lin", name="fov_lat") self.data = NDDataArray( axes=[energy_axis, fov_lon_axis, fov_lat_axis], data=data, interp_kwargs=interp_kwargs, ) self.meta = meta or {} def __str__(self): ss = self.__class__.__name__ ss += f"\n{self.data}" return ss
[docs] @classmethod def from_table(cls, table): """Read from `~astropy.table.Table`.""" # Spec says key should be "BKG", but there are files around # (e.g. CTA 1DC) that use "BGD". For now we support both if "BKG" in table.colnames: bkg_name = "BKG" elif "BGD" in table.colnames: bkg_name = "BGD" else: raise ValueError('Invalid column names. Need "BKG" or "BGD".') # Currently some files (e.g. CTA 1DC) contain unit in the FITS file # '1/s/MeV/sr', which is invalid ( try: astropy.units.Unit('1/s/MeV/sr') # This should be corrected. # For now, we hard-code the unit here: data_unit = u.Unit("s-1 MeV-1 sr-1") return cls( energy_lo=table["ENERG_LO"].quantity[0], energy_hi=table["ENERG_HI"].quantity[0], fov_lon_lo=table["DETX_LO"].quantity[0], fov_lon_hi=table["DETX_HI"].quantity[0], fov_lat_lo=table["DETY_LO"].quantity[0], fov_lat_hi=table["DETY_HI"].quantity[0], data=table[bkg_name].data[0] * data_unit, meta=table.meta, )
[docs] @classmethod def from_hdulist(cls, hdulist, hdu="BACKGROUND"): """Create from `~astropy.io.fits.HDUList`.""" return cls.from_table(Table.read(hdulist[hdu]))
[docs] @classmethod def read(cls, filename, hdu="BACKGROUND"): """Read from file.""" with fits.open(make_path(filename), memmap=False) as hdulist: return cls.from_hdulist(hdulist, hdu=hdu)
[docs] def to_table(self): """Convert to `~astropy.table.Table`.""" meta = self.meta.copy() detx = self.data.axis("fov_lon").edges dety = self.data.axis("fov_lat").edges energy = self.data.axis("energy").edges table = Table(meta=meta) table["DETX_LO"] = detx[:-1][np.newaxis] table["DETX_HI"] = detx[1:][np.newaxis] table["DETY_LO"] = dety[:-1][np.newaxis] table["DETY_HI"] = dety[1:][np.newaxis] table["ENERG_LO"] = energy[:-1][np.newaxis] table["ENERG_HI"] = energy[1:][np.newaxis] table["BKG"] = self.data.data[np.newaxis] return table
[docs] def to_fits(self, name="BACKGROUND"): """Convert to `~astropy.io.fits.BinTableHDU`.""" return fits.BinTableHDU(self.to_table(), name=name)
[docs] def evaluate(self, fov_lon, fov_lat, energy_reco, method="linear", **kwargs): """Evaluate at given FOV position and energy. Parameters ---------- fov_lon, fov_lat : `~astropy.coordinates.Angle` FOV coordinates expecting in AltAz frame. energy_reco : `~astropy.units.Quantity` energy on which you want to interpolate. Same dimension than fov_lat and fov_lat method : str {'linear', 'nearest'}, optional Interpolation method kwargs : dict option for interpolation for `~scipy.interpolate.RegularGridInterpolator` Returns ------- array : `~astropy.units.Quantity` Interpolated values, axis order is the same as for the NDData array """ values = self.data.evaluate( fov_lon=fov_lon, fov_lat=fov_lat, energy=energy_reco, method=method, **kwargs, ) return values
[docs] def evaluate_integrate( self, fov_lon, fov_lat, energy_reco, method="linear", **kwargs ): """Integrate in a given energy band. Parameters ---------- fov_lon, fov_lat : `~astropy.coordinates.Angle` FOV coordinates expecting in AltAz frame. energy_reco: `~astropy.units.Quantity` Reconstructed energy edges. method : {'linear', 'nearest'}, optional Interpolation method Returns ------- array : `~astropy.units.Quantity` Returns 2D array with axes offset """ data = self.evaluate(fov_lon, fov_lat, energy_reco, method=method) return trapz_loglog(data, energy_reco, axis=0)
[docs] def to_2d(self): """Convert to `Background2D`. This takes the values at Y = 0 and X >= 0. """ idx_lon = self.data.axis("fov_lon").coord_to_idx(0 * u.deg)[0] idx_lat = self.data.axis("fov_lat").coord_to_idx(0 * u.deg)[0] data = self.data.data[:, idx_lon:, idx_lat].copy() energy = self.data.axis("energy").edges offset = self.data.axis("fov_lon").edges[idx_lon:] return Background2D( energy_lo=energy[:-1], energy_hi=energy[1:], offset_lo=offset[:-1], offset_hi=offset[1:], data=data, )
[docs] def peek(self, figsize=(10, 8)): return self.to_2d().peek(figsize)
[docs]class Background2D: """Background 2D. Data format specification: :ref:`gadf:bkg_2d` Parameters ---------- energy_lo, energy_hi : `~astropy.units.Quantity` Energy binning offset_lo, offset_hi : `~astropy.units.Quantity` FOV coordinate offset-axis binning data : `~astropy.units.Quantity` Background rate (usually: ``s^-1 MeV^-1 sr^-1``) """ default_interp_kwargs = dict(bounds_error=False, fill_value=None) """Default Interpolation kwargs for `~gammapy.utils.nddata.NDDataArray`. Extrapolate.""" def __init__( self, energy_lo, energy_hi, offset_lo, offset_hi, data, meta=None, interp_kwargs=None, ): if interp_kwargs is None: interp_kwargs = self.default_interp_kwargs e_edges = edges_from_lo_hi(energy_lo, energy_hi) energy_axis = MapAxis.from_edges(e_edges, interp="log", name="energy") offset_edges = edges_from_lo_hi(offset_lo, offset_hi) offset_axis = MapAxis.from_edges(offset_edges, interp="lin", name="offset") self.data = NDDataArray( axes=[energy_axis, offset_axis], data=data, interp_kwargs=interp_kwargs ) self.meta = meta or {} def __str__(self): ss = self.__class__.__name__ ss += f"\n{self.data}" return ss
[docs] @classmethod def from_table(cls, table): """Read from `~astropy.table.Table`.""" # Spec says key should be "BKG", but there are files around # (e.g. CTA 1DC) that use "BGD". For now we support both if "BKG" in table.colnames: bkg_name = "BKG" elif "BGD" in table.colnames: bkg_name = "BGD" else: raise ValueError('Invalid column names. Need "BKG" or "BGD".') # Currently some files (e.g. CTA 1DC) contain unit in the FITS file # '1/s/MeV/sr', which is invalid ( try: astropy.units.Unit('1/s/MeV/sr') # This should be corrected. # For now, we hard-code the unit here: data_unit = u.Unit("s-1 MeV-1 sr-1") return cls( energy_lo=table["ENERG_LO"].quantity[0], energy_hi=table["ENERG_HI"].quantity[0], offset_lo=table["THETA_LO"].quantity[0], offset_hi=table["THETA_HI"].quantity[0], data=table[bkg_name].data[0] * data_unit, meta=table.meta, )
[docs] @classmethod def from_hdulist(cls, hdulist, hdu="BACKGROUND"): """Create from `~astropy.io.fits.HDUList`.""" return cls.from_table(Table.read(hdulist[hdu]))
[docs] @classmethod def read(cls, filename, hdu="BACKGROUND"): """Read from file.""" with fits.open(make_path(filename), memmap=False) as hdulist: return cls.from_hdulist(hdulist, hdu=hdu)
[docs] def to_table(self): """Convert to `~astropy.table.Table`.""" meta = self.meta.copy() table = Table(meta=meta) theta = self.data.axis("offset").edges energy = self.data.axis("energy").edges table["THETA_LO"] = theta[:-1][np.newaxis] table["THETA_HI"] = theta[1:][np.newaxis] table["ENERG_LO"] = energy[:-1][np.newaxis] table["ENERG_HI"] = energy[1:][np.newaxis] table["BKG"] = self.data.data[np.newaxis] return table
[docs] def to_fits(self, name="BACKGROUND"): """Convert to `~astropy.io.fits.BinTableHDU`.""" return fits.BinTableHDU(self.to_table(), name=name)
[docs] def evaluate(self, fov_lon, fov_lat, energy_reco, method="linear", **kwargs): """Evaluate at a given FOV position and energy. The fov_lon, fov_lat, energy_reco has to have the same shape since this is a set of points on which you want to evaluate. To have the same API than background 3D for the background evaluation, the offset is ``fov_altaz_lon``. Parameters ---------- fov_lon, fov_lat : `~astropy.coordinates.Angle` FOV coordinates expecting in AltAz frame, same shape than energy_reco energy_reco : `~astropy.units.Quantity` Reconstructed energy, same dimension than fov_lat and fov_lat method : str {'linear', 'nearest'}, optional Interpolation method kwargs : dict option for interpolation for `~scipy.interpolate.RegularGridInterpolator` Returns ------- array : `~astropy.units.Quantity` Interpolated values, axis order is the same as for the NDData array """ offset = np.sqrt(fov_lon ** 2 + fov_lat ** 2) return self.data.evaluate( offset=offset, energy=energy_reco, method=method, **kwargs )
[docs] def evaluate_integrate(self, fov_lon, fov_lat, energy_reco, method="linear"): """Evaluate at given FOV position and energy, by integrating over the energy range. Parameters ---------- fov_lon, fov_lat : `~astropy.coordinates.Angle` FOV coordinates expecting in AltAz frame. energy_reco: `~astropy.units.Quantity` Reconstructed energy edges. method : {'linear', 'nearest'}, optional Interpolation method Returns ------- array : `~astropy.units.Quantity` Returns 2D array with axes offset """ data = self.evaluate(fov_lon, fov_lat, energy_reco, method=method) return trapz_loglog(data, energy_reco, axis=0)
[docs] def to_3d(self): """Convert to `Background3D`. Fill in a radially symmetric way. """ raise NotImplementedError
[docs] def plot(self, ax=None, add_cbar=True, **kwargs): """Plot energy offset dependence of the background model. """ import matplotlib.pyplot as plt from matplotlib.colors import LogNorm ax = plt.gca() if ax is None else ax x = self.data.axis("energy").edges.to_value("TeV") y = self.data.axis("offset").edges.to_value("deg") z = self.data.data.T.value kwargs.setdefault("cmap", "GnBu") kwargs.setdefault("edgecolors", "face") caxes = ax.pcolormesh(x, y, z, norm=LogNorm(), **kwargs) ax.set_xscale("log") ax.set_ylabel(f"Offset (deg)") ax.set_xlabel(f"Energy (TeV)") xmin, xmax = x.min(), x.max() ax.set_xlim(xmin, xmax) if add_cbar: label = f"Background rate ({self.data.data.unit})" ax.figure.colorbar(caxes, ax=ax, label=label)
[docs] def plot_offset_dependence(self, ax=None, offset=None, energy=None, **kwargs): """Plot background rate versus offset for a given energy. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Axis offset : `~astropy.coordinates.Angle` Offset axis energy : `~astropy.units.Quantity` Energy Returns ------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax if energy is None: e_min, e_max = np.log10(self.data.axis("energy").center.value[[0, -1]]) energy = np.logspace(e_min, e_max, 4) * self.data.axis("energy").unit if offset is None: offset = self.data.axis("offset").center for ee in energy: bkg = self.data.evaluate(offset=offset, energy=ee) if np.isnan(bkg).all(): continue label = f"energy = {ee:.1f}" ax.plot(offset, bkg.value, label=label, **kwargs) ax.set_xlabel(f"Offset ({self.data.axis('offset').unit})") ax.set_ylabel(f"Background rate ({self.data.data.unit})") ax.set_yscale("log") ax.legend(loc="upper right") return ax
[docs] def plot_energy_dependence(self, ax=None, offset=None, energy=None, **kwargs): """Plot background rate versus energy for a given offset. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Axis offset : `~astropy.coordinates.Angle` Offset energy : `~astropy.units.Quantity` Energy axis kwargs : dict Forwarded tp plt.plot() Returns ------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax if offset is None: off_min, off_max = self.data.axis("offset").center.value[[0, -1]] offset = np.linspace(off_min, off_max, 4) * self.data.axis("offset").unit if energy is None: energy = self.data.axis("energy").center for off in offset: bkg = self.data.evaluate(offset=off, energy=energy) label = f"offset = {off:.1f}" ax.plot(energy, bkg.value, label=label, **kwargs) ax.set_xscale("log") ax.set_yscale("log") ax.set_xlabel(f"Energy [{energy.unit}]") ax.set_ylabel(f"Background rate ({self.data.data.unit})") ax.set_xlim(min(energy.value), max(energy.value)) ax.legend(loc="best") return ax
[docs] def plot_spectrum(self, ax=None, **kwargs): """Plot angle integrated background rate versus energy. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Axis kwargs : dict Forwarded tp plt.plot() Returns ------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax offset = self.data.axis("offset").edges energy = self.data.axis("energy").center bkg = [] for ee in energy: data = self.data.evaluate(offset=offset, energy=ee) val = np.nansum(trapz_loglog(data, offset, axis=0)) bkg.append(val.value) ax.plot(energy, bkg, label="integrated spectrum", **kwargs) unit = self.data.data.unit * offset.unit * offset.unit ax.set_xscale("log") ax.set_yscale("log") ax.set_xlabel(f"Energy [{energy.unit}]") ax.set_ylabel(f"Background rate ({unit})") ax.set_xlim(min(energy.value), max(energy.value)) ax.legend(loc="best") return ax
[docs] def peek(self, figsize=(10, 8)): """Quick-look summary plots.""" import matplotlib.pyplot as plt fig, axes = plt.subplots(nrows=2, ncols=2, figsize=figsize) self.plot(ax=axes[1][1]) self.plot_offset_dependence(ax=axes[0][0]) self.plot_energy_dependence(ax=axes[1][0]) self.plot_spectrum(ax=axes[0][1]) plt.tight_layout()