# Licensed under a 3-clause BSD style license - see LICENSE.rstimportabcimportloggingimportnumpyasnpfromastropyimportunitsasufromgammapy.mapsimportMapAxes,MapAxisfromgammapy.utils.gaussimportMultiGauss2Dfromgammapy.utils.interpolationimportScaledRegularGridInterpolatorfrom.coreimportPSF__all__=["ParametricPSF","EnergyDependentMultiGaussPSF","PSFKing"]log=logging.getLogger(__name__)
[docs]classParametricPSF(PSF):"""Parametric PSF base class. Parameters ---------- axes : list of `MapAxis` or `MapAxes` Axes. data : dict of `~numpy.ndarray` or `~numpy.recarray` Data. unit : dict of str or `~astropy.units.Unit` Unit. meta : dict Metadata dictionary. """@property@abc.abstractmethoddefrequired_parameters(self):return[]
@propertydefquantity(self):"""Quantity."""quantity={}fornameinself.required_parameters:quantity[name]=self.data[name]*self.unit[name]returnquantity@propertydefunit(self):"""Map unit as a `~astropy.units.Unit`."""returnself._unit
[docs]defto_unit(self,unit):"""Convert IRF to unit."""raiseNotImplementedError
[docs]defevaluate_parameters(self,energy_true,offset):"""Evaluate analytic PSF parameters at a given energy and offset. Uses nearest-neighbor interpolation. Parameters ---------- energy_true : `~astropy.units.Quantity` Energy value. offset : `~astropy.coordinates.Angle` Offset in the field of view. Returns ------- values : `~astropy.units.Quantity` Interpolated value. """pars={}fornameinself.required_parameters:value=self._interpolators[name]((energy_true,offset))pars[name]=valuereturnpars
[docs]defto_table(self,format="gadf-dl3"):"""Convert PSF table data to table. Parameters ---------- format : {"gadf-dl3"} Format specification. Default is "gadf-dl3". Returns ------- hdu_list : `~astropy.io.fits.HDUList` PSF in HDU list format. """fromgammapy.irf.ioimportIRF_DL3_HDU_SPECIFICATIONtable=self.axes.to_table(format="gadf-dl3")spec=IRF_DL3_HDU_SPECIFICATION[self.tag]["column_name"]fornameinself.required_parameters:column_name=spec[name]table[column_name]=self.data[name].T[np.newaxis]table[column_name].unit=self.unit[name]# Create hdu and hdu listreturntable
[docs]@classmethoddeffrom_table(cls,table,format="gadf-dl3"):"""Create parametric PSF from `~astropy.table.Table`. Parameters ---------- table : `~astropy.table.Table` Table information. format : {"gadf-dl3"}, optional Format specification. Default is "gadf-dl3". Returns ------- psf : `~ParametricPSF` PSF class. """fromgammapy.irf.ioimportIRF_DL3_HDU_SPECIFICATIONaxes=MapAxes.from_table(table,format=format)[cls.required_axes]dtype={"names":cls.required_parameters,"formats":len(cls.required_parameters)*(np.float32,),}data=np.empty(axes.shape,dtype=dtype)unit={}spec=IRF_DL3_HDU_SPECIFICATION[cls.tag]["column_name"]fornameincls.required_parameters:column=table[spec[name]]values=column.data[0].transpose()# TODO: this fixes some files where sigma is written as zeroif"sigma"inname:values[values==0]=1.0data[name]=values.reshape(axes.shape)unit[name]=column.unitor""unit={key:u.Unit(val)forkey,valinunit.items()}returncls(axes=axes,data=data,meta=table.meta.copy(),unit=unit)
[docs]defto_psf3d(self,rad=None):"""Create a PSF3D from a parametric PSF. It will be defined on the same energy and offset values than the input PSF. Parameters ---------- rad : `~astropy.units.Quantity` Rad values. Returns ------- psf3d : `~gammapy.irf.PSF3D` 3D PSF. """fromgammapy.datasets.mapimportRAD_AXIS_DEFAULTfromgammapy.irfimportPSF3Doffset_axis=self.axes["offset"]energy_axis_true=self.axes["energy_true"]ifradisNone:rad_axis=RAD_AXIS_DEFAULTelse:rad_axis=MapAxis.from_edges(rad,name="rad")axes=MapAxes([energy_axis_true,offset_axis,rad_axis])data=self.evaluate(**axes.get_coord())returnPSF3D(axes=axes,data=data.value,unit=data.unit,meta=self.meta.copy())
[docs]defcontainment(self,rad,**kwargs):"""Containment of the PSF at given axes coordinates. Parameters ---------- rad : `~astropy.units.Quantity` Rad value. **kwargs : dict Other coordinates. Returns ------- containment : `~numpy.ndarray` Containment. """pars=self.evaluate_parameters(**kwargs)containment=self.evaluate_containment(rad=rad,**pars)returncontainment
[docs]defevaluate(self,rad,**kwargs):"""Evaluate the PSF model. Parameters ---------- rad : `~astropy.coordinates.Angle` Offset from PSF center used for evaluating the PSF on a grid. **kwargs : dict Other coordinates. Returns ------- psf_value : `~astropy.units.Quantity` PSF value. """pars=self.evaluate_parameters(**kwargs)value=self.evaluate_direct(rad=rad,**pars)returnvalue
[docs]defis_allclose(self,other,rtol_axes=1e-3,atol_axes=1e-6,**kwargs):"""Compare two data IRFs for equivalency. Parameters ---------- other : `gammapy.irfs.ParametricPSF` The PSF to compare against. rtol_axes : float, optional Relative tolerance for the axis comparison. Default is 1e-3. atol_axes : float, optional Relative tolerance for the axis comparison. Default is 1e-6. **kwargs : dict Keywords passed to `numpy.allclose`. Returns ------- is_allclose : bool Whether the IRF is all close. """ifnotisinstance(other,self.__class__):returnTypeError(f"Cannot compare {type(self)} and {type(other)}")data_eq=Trueforkeyinself.quantity.keys():ifself.quantity[key].shape!=other.quantity[key].shape:returnFalsedata_eq&=np.allclose(self.quantity[key],other.quantity[key],**kwargs)axes_eq=self.axes.is_allclose(other.axes,rtol=rtol_axes,atol=atol_axes)returnaxes_eqanddata_eq
defget_sigmas_and_norms(**kwargs):"""Convert scale and amplitude to norms."""sigmas=u.Quantity([kwargs[f"sigma_{idx}"]foridxin[1,2,3]])scale=kwargs["scale"]ones=np.ones(scale.shape)amplitudes=u.Quantity([ones,kwargs["ampl_2"],kwargs["ampl_3"]])norms=2*scale*amplitudes*sigmas**2returnsigmas,norms
[docs]classEnergyDependentMultiGaussPSF(ParametricPSF):"""Triple Gauss analytical PSF depending on true energy and offset. Parameters ---------- axes : list of `MapAxis` Required axes are ["energy_true", "offset"]. data : `~numpy.recarray` Data array. meta : dict Metadata dictionary. Examples -------- Plot R68 of the PSF vs. offset and true energy: .. plot:: :include-source: import matplotlib.pyplot as plt from gammapy.irf import EnergyDependentMultiGaussPSF filename = '$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits' psf = EnergyDependentMultiGaussPSF.read(filename, hdu='POINT SPREAD FUNCTION') psf.plot_containment_radius(fraction=0.68) plt.show() """tag="psf_3gauss"required_axes=["energy_true","offset"]required_parameters=["sigma_1","sigma_2","sigma_3","scale","ampl_2","ampl_3"]
[docs]@staticmethoddefevaluate_containment(rad,**kwargs):"""Containment of the PSF at given axes coordinates. Parameters ---------- rad : `~astropy.units.Quantity` Rad value. **kwargs : dict Parameters, see `required_parameters`. Returns ------- containment : `~numpy.ndarray` Containment. """sigmas,norms=get_sigmas_and_norms(**kwargs)m=MultiGauss2D(sigmas=sigmas,norms=norms)m.normalize()containment=m.containment_fraction(rad)returncontainment
[docs]@staticmethoddefevaluate_direct(rad,**kwargs):"""Evaluate PSF model. Parameters ---------- rad : `~astropy.units.Quantity` Rad value. **kwargs : dict Parameters, see `required_parameters`. Returns ------- value : `~numpy.ndarray` PSF value. """sigmas,norms=get_sigmas_and_norms(**kwargs)m=MultiGauss2D(sigmas=sigmas,norms=norms)m.normalize()returnm(rad)
[docs]classPSFKing(ParametricPSF):"""King profile analytical PSF depending on energy and offset. This PSF parametrisation and FITS data format is described here: :ref:`gadf:psf_king`. Parameters ---------- axes : list of `MapAxis` or `MapAxes` Data axes, required are ["energy_true", "offset"]. meta : dict Metadata dictionary. """tag="psf_king"required_axes=["energy_true","offset"]required_parameters=["gamma","sigma"]default_interp_kwargs=dict(bounds_error=False,fill_value=None)
[docs]@staticmethoddefevaluate_containment(rad,gamma,sigma):"""Containment of the PSF at given axes coordinates. Parameters ---------- rad : `~astropy.units.Quantity` Rad value. gamma : `~astropy.units.Quantity` Gamma parameter. sigma : `~astropy.units.Quantity` Sigma parameter. Returns ------- containment : `~numpy.ndarray` Containment. """withnp.errstate(divide="ignore",invalid="ignore"):powterm=1-gammaterm=(1+rad**2/(2*gamma*sigma**2))**powtermcontainment=1-termreturncontainment
[docs]@staticmethoddefevaluate_direct(rad,gamma,sigma):"""Evaluate the PSF model. Formula is given here: :ref:`gadf:psf_king`. Parameters ---------- rad : `~astropy.coordinates.Angle` Offset from PSF center used for evaluating the PSF on a grid. Returns ------- psf_value : `~astropy.units.Quantity` PSF value. """withnp.errstate(divide="ignore"):term1=1/(2*np.pi*sigma**2)term2=1-1/gammaterm3=(1+rad**2/(2*gamma*sigma**2))**(-gamma)returnterm1*term2*term3