Source code for gammapy.irf.psf_map

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
import astropy.units as u
from gammapy.maps import Map, MapAxis, MapCoord, WcsGeom
from gammapy.modeling.models import PowerLawSpectralModel
from gammapy.utils.random import InverseCDFSampler, get_random_state
from .irf_map import IRFMap
from .psf_kernel import PSFKernel
from .psf_table import EnergyDependentTablePSF

__all__ = ["PSFMap"]


[docs]class PSFMap(IRFMap): """Class containing the Map of PSFs and allowing to interact with it. Parameters ---------- psf_map : `~gammapy.maps.Map` the input PSF Map. Should be a Map with 2 non spatial axes. rad and true energy axes should be given in this specific order. exposure_map : `~gammapy.maps.Map` Associated exposure map. Needs to have a consistent map geometry. Examples -------- :: import numpy as np from astropy import units as u from astropy.coordinates import SkyCoord from gammapy.maps import Map, WcsGeom, MapAxis from gammapy.irf import EnergyDependentMultiGaussPSF, EffectiveAreaTable2D from gammapy.cube import make_psf_map, PSFMap, make_map_exposure_true_energy # Define energy axis. Note that the name is fixed. energy_axis = MapAxis.from_edges(np.logspace(-1., 1., 4), unit='TeV', name='energy') # Define rad axis. Again note the axis name rads = np.linspace(0., 0.5, 100) * u.deg rad_axis = MapAxis.from_edges(rads, unit='deg', name='theta') # Define parameters pointing = SkyCoord(0., 0., unit='deg') max_offset = 4 * u.deg # Create WcsGeom geom = WcsGeom.create(binsz=0.25*u.deg, width=10*u.deg, skydir=pointing, axes=[rad_axis, energy_axis]) # Extract EnergyDependentTablePSF from CTA 1DC IRF filename = '$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits' psf = EnergyDependentMultiGaussPSF.read(filename, hdu='POINT SPREAD FUNCTION') psf3d = psf.to_psf3d(rads) aeff2d = EffectiveAreaTable2D.read(filename, hdu='EFFECTIVE AREA') # Create the exposure map exposure_geom = geom.to_image().to_cube([energy_axis]) exposure_map = make_map_exposure_true_energy(pointing, "1 h", aeff2d, exposure_geom) # create the PSFMap for the specified pointing psf_map = make_psf_map(psf3d, pointing, geom, max_offset, exposure_map) # Get an EnergyDependentTablePSF at any position in the image psf_table = psf_map.get_energy_dependent_table_psf(SkyCoord(2., 2.5, unit='deg')) # Write map to disk psf_map.write('psf_map.fits') """ _hdu_name = "psf" @property def psf_map(self): return self._irf_map @psf_map.setter def psf_map(self, value): self._irf_map = value def __init__(self, psf_map, exposure_map=None): if psf_map.geom.axes[1].name.upper() != "ENERGY_TRUE": raise ValueError("Incorrect energy axis position in input Map") if psf_map.geom.axes[0].name.upper() != "THETA": raise ValueError("Incorrect theta axis position in input Map") super().__init__(irf_map=psf_map, exposure_map=exposure_map)
[docs] def get_energy_dependent_table_psf(self, position): """Get energy-dependent PSF at a given position. Parameters ---------- position : `~astropy.coordinates.SkyCoord` the target position. Should be a single coordinates Returns ------- psf_table : `~gammapy.irf.EnergyDependentTablePSF` the table PSF """ if position.size != 1: raise ValueError( "EnergyDependentTablePSF can be extracted at one single position only." ) energy = self.psf_map.geom.get_axis_by_name("energy_true").center rad = self.psf_map.geom.get_axis_by_name("theta").center coords = { "skycoord": position, "energy_true": energy.reshape((-1, 1, 1, 1)), "theta": rad.reshape((1, -1, 1, 1)), } data = self.psf_map.interp_by_coord(coords) psf_values = u.Quantity(data[:, :, 0, 0], unit=self.psf_map.unit, copy=False) if self.exposure_map is not None: coords = { "skycoord": position, "energy_true": energy.reshape((-1, 1, 1)), "theta": 0 * u.deg, } data = self.exposure_map.interp_by_coord(coords).squeeze() exposure = data * self.exposure_map.unit else: exposure = None # Beware. Need to revert rad and energies to follow the TablePSF scheme. return EnergyDependentTablePSF( energy=energy, rad=rad, psf_value=psf_values, exposure=exposure )
[docs] def get_psf_kernel(self, position, geom, max_radius=None, factor=4): """Returns a PSF kernel at the given position. The PSF is returned in the form a WcsNDMap defined by the input Geom. Parameters ---------- position : `~astropy.coordinates.SkyCoord` the target position. Should be a single coordinate geom : `~gammapy.maps.Geom` the target geometry to use max_radius : `~astropy.coordinates.Angle` maximum angular size of the kernel map factor : int oversampling factor to compute the PSF Returns ------- kernel : `~gammapy.cube.PSFKernel` the resulting kernel """ table_psf = self.get_energy_dependent_table_psf(position) if max_radius is None: max_radius = np.max(table_psf.rad) return PSFKernel.from_table_psf(table_psf, geom, max_radius, factor)
[docs] def containment_radius_map(self, energy, fraction=0.68): """Containment radius map. Parameters ---------- energy : `~astropy.units.Quantity` Scalar energy at which to compute the containment radius fraction : float the containment fraction (range: 0 to 1) Returns ------- containment_radius_map : `~gammapy.maps.Map` Containment radius map """ coords = self.psf_map.geom.to_image().get_coord().skycoord.flatten() m = Map.from_geom(self.psf_map.geom.to_image(), unit="deg") for coord in coords: psf_table = self.get_energy_dependent_table_psf(coord) containment_radius = psf_table.containment_radius(energy, fraction) m.fill_by_coord(coord, containment_radius) return m
[docs] @classmethod def from_geom(cls, geom): """Create psf map from geom. Parameters ---------- geom : `Geom` PSF map geometry. Returns ------- psf_map : `PSFMap` Point spread function map. """ geom_exposure_psf = geom.squash(axis="theta") exposure_psf = Map.from_geom(geom_exposure_psf, unit="m2 s") psf_map = Map.from_geom(geom, unit="sr-1") return cls(psf_map, exposure_psf)
[docs] def sample_coord(self, map_coord, random_state=0): """Apply PSF corrections on the coordinates of a set of simulated events. Parameters ---------- map_coord : `~gammapy.maps.MapCoord` object. Sequence of coordinates and energies of sampled events. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Returns ------- corr_coord : `~gammapy.maps.MapCoord` object. Sequence of PSF-corrected coordinates of the input map_coord map. """ random_state = get_random_state(random_state) rad_axis = self.psf_map.geom.get_axis_by_name("theta") coord = { "skycoord": map_coord.skycoord.reshape(-1, 1), "energy_true": map_coord["energy_true"].reshape(-1, 1), "theta": rad_axis.center, } pdf = ( self.psf_map.interp_by_coord(coord) * rad_axis.center.value * rad_axis.bin_width.value ) sample_pdf = InverseCDFSampler(pdf, axis=1, random_state=random_state) pix_coord = sample_pdf.sample_axis() separation = rad_axis.pix_to_coord(pix_coord) position_angle = random_state.uniform(360, size=len(map_coord.lon)) * u.deg event_positions = map_coord.skycoord.directional_offset_by( position_angle=position_angle, separation=separation ) return MapCoord.create( {"skycoord": event_positions, "energy_true": map_coord["energy_true"]} )
[docs] @classmethod def from_energy_dependent_table_psf(cls, table_psf): """Create PSF map from table PSF object. Helper function to create an allsky PSF map from table PSF, which does not depend on position. Parameters ---------- table_psf : `EnergyDependentTablePSF` Table PSF Returns ------- psf_map : `PSFMap` Point spread function map. """ energy_axis = MapAxis.from_nodes( table_psf.energy, name="energy_true", interp="log" ) rad_axis = MapAxis.from_nodes(table_psf.rad, name="theta") geom = WcsGeom.create( npix=(2, 1), proj="CAR", binsz=180, axes=[rad_axis, energy_axis] ) coords = geom.get_coord() # TODO: support broadcasting in .evaluate() data = table_psf._interpolate( (coords["energy_true"], coords["theta"]) ).to_value("sr-1") psf_map = Map.from_geom(geom, data=data, unit="sr-1") geom_exposure = geom.squash(axis="theta") data = table_psf.exposure.reshape((-1, 1, 1, 1)) exposure_map = Map.from_geom(geom_exposure, unit="cm2 s") exposure_map.quantity += data return cls(psf_map=psf_map, exposure_map=exposure_map)
[docs] def to_image(self, spectrum=None, keepdims=True): """Reduce to a 2-D map after weighing with the associated exposure and a spectrum Parameters ---------- spectrum : `~gammapy.modeling.models.SpectralModel`, optional Spectral model to compute the weights. Default is power-law with spectral index of 2. keepdims : bool, optional If True, the energy axis is kept with one bin. If False, the axis is removed Returns ------- psf_out : `PSFMap` `PSFMap` with the energy axis summed over """ from gammapy.makers.utils import _map_spectrum_weight if spectrum is None: spectrum = PowerLawSpectralModel(index=2.0) exp_weighed = _map_spectrum_weight(self.exposure_map, spectrum) exposure = exp_weighed.sum_over_axes(axes=["energy_true"], keepdims=keepdims) psf_data = exp_weighed.data * self.psf_map.data / exposure.data psf_map = Map.from_geom(geom=self.psf_map.geom, data=psf_data, unit="sr-1") psf = psf_map.sum_over_axes(axes=["energy_true"], keepdims=keepdims) return self.__class__(psf_map=psf, exposure_map=exposure)