Source code for gammapy.cube.exposure

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
from gammapy.maps import WcsNDMap
from gammapy.modeling.models import PowerLawSpectralModel

__all__ = ["make_map_exposure_true_energy"]


[docs]def make_map_exposure_true_energy(pointing, livetime, aeff, geom): """Compute exposure map. This map has a true energy axis, the exposure is not combined with energy dispersion. Parameters ---------- pointing : `~astropy.coordinates.SkyCoord` Pointing direction livetime : `~astropy.units.Quantity` Livetime aeff : `~gammapy.irf.EffectiveAreaTable2D` Effective area geom : `~gammapy.maps.WcsGeom` Map geometry (must have an energy axis) Returns ------- map : `~gammapy.maps.WcsNDMap` Exposure map """ offset = geom.separation(pointing) energy = geom.get_axis_by_name("energy").center exposure = aeff.data.evaluate( offset=offset, energy=energy[:, np.newaxis, np.newaxis] ) # TODO: Improve IRF evaluate to preserve energy axis if length 1 # For now, we handle that case via this hack: if len(exposure.shape) < 3: exposure = np.expand_dims(exposure.value, 0) * exposure.unit exposure = (exposure * livetime).to("m2 s") return WcsNDMap(geom, exposure.value.reshape(geom.data_shape), unit=exposure.unit)
def _map_spectrum_weight(map, spectrum=None): """Weight a map with a spectrum. This requires map to have an "energy" axis. The weights are normalised so that they sum to 1. The mean and unit of the output image is the same as of the input cube. At the moment this is used to get a weighted exposure image. Parameters ---------- map : `~gammapy.maps.Map` Input map with an "energy" axis. spectrum : `~gammapy.modeling.models.SpectralModel` Spectral model to compute the weights. Default is power-law with spectral index of 2. Returns ------- map_weighted : `~gammapy.maps.Map` Weighted image """ if spectrum is None: spectrum = PowerLawSpectralModel(index=2.0) # Compute weights vector energy_edges = map.geom.get_axis_by_name("energy").edges weights = spectrum.integral( emin=energy_edges[:-1], emax=energy_edges[1:], intervals=True ) weights /= weights.sum() shape = np.ones(len(map.geom.data_shape)) shape[0] = -1 return map * weights.reshape(shape.astype(int))