# Licensed under a 3-clause BSD style license - see LICENSE.rstimportnumpyasnpimportastropy.unitsasufromastropy.visualizationimportquantity_supportimportmatplotlib.pyplotaspltfromgammapy.mapsimportMapAxes,MapAxisfromgammapy.maps.axesimportUNIT_STRING_FORMATfromgammapy.visualization.utilsimportadd_colorbarfrom.coreimportIRF__all__=["EffectiveAreaTable2D"]
[docs]classEffectiveAreaTable2D(IRF):"""2D effective area table. Data format specification: :ref:`gadf:aeff_2d`. Parameters ---------- energy_axis_true : `MapAxis` True energy axis. offset_axis : `MapAxis` Field of view offset axis. data : `~astropy.units.Quantity` Effective area. meta : dict Metadata dictionary. Examples -------- Here's an example you can use to learn about this class: >>> from gammapy.irf import EffectiveAreaTable2D >>> filename = '$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits' >>> aeff = EffectiveAreaTable2D.read(filename, hdu='EFFECTIVE AREA') >>> print(aeff) EffectiveAreaTable2D -------------------- <BLANKLINE> axes : ['energy_true', 'offset'] shape : (42, 6) ndim : 2 unit : m2 dtype : >f4 <BLANKLINE> Here's another one, created from scratch, without reading a file: >>> from gammapy.irf import EffectiveAreaTable2D >>> from gammapy.maps import MapAxis >>> energy_axis_true = MapAxis.from_energy_bounds( ... "0.1 TeV", "100 TeV", nbin=30, name="energy_true" ... ) >>> offset_axis = MapAxis.from_bounds(0, 5, nbin=4, name="offset") >>> aeff = EffectiveAreaTable2D(axes=[energy_axis_true, offset_axis], data=1e10, unit="cm2") >>> print(aeff) EffectiveAreaTable2D -------------------- <BLANKLINE> axes : ['energy_true', 'offset'] shape : (30, 4) ndim : 2 unit : cm2 dtype : float64 <BLANKLINE> """tag="aeff_2d"required_axes=["energy_true","offset"]default_unit=u.m**2
[docs]defplot_energy_dependence(self,ax=None,offset=None,**kwargs):"""Plot effective area versus energy for a given offset. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. Default is None. offset : `~astropy.coordinates.Angle`, optional Offset. Default is None. kwargs : dict Forwarded to plt.plot(). Returns ------- ax : `~matplotlib.axes.Axes` Matplotlib axes. """ax=plt.gca()ifaxisNoneelseaxifoffsetisNone:off_min,off_max=self.axes["offset"].boundsoffset=np.linspace(off_min,off_max,4)energy_axis=self.axes["energy_true"]foroffinoffset:area=self.evaluate(offset=off,energy_true=energy_axis.center)label=kwargs.pop("label",f"offset = {off:.1f}")withquantity_support():ax.plot(energy_axis.center,area,label=label,**kwargs)energy_axis.format_plot_xaxis(ax=ax)ax.set_ylabel(f"Effective Area [{ax.yaxis.units.to_string(UNIT_STRING_FORMAT)}]")ax.legend()returnax
[docs]defplot_offset_dependence(self,ax=None,energy=None,**kwargs):"""Plot effective area versus offset for a given energy. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. Default is None. energy : `~astropy.units.Quantity` Energy. **kwargs : dict Keyword argument passed to `~matplotlib.pyplot.plot`. Returns ------- ax : `~matplotlib.axes.Axes` Matplotlib axes. """ax=plt.gca()ifaxisNoneelseaxifenergyisNone:energy_axis=self.axes["energy_true"]e_min,e_max=energy_axis.center[[0,-1]]energy=np.geomspace(e_min,e_max,4)offset_axis=self.axes["offset"]foreeinenergy:area=self.evaluate(offset=offset_axis.center,energy_true=ee)area/=np.nanmax(area)ifnp.isnan(area).all():continuelabel=f"energy = {ee:.1f}"withquantity_support():ax.plot(offset_axis.center,area,label=label,**kwargs)offset_axis.format_plot_xaxis(ax=ax)ax.set_ylim(0,1.1)ax.set_ylabel("Relative Effective Area")ax.legend(loc="best")returnax
[docs]defplot(self,ax=None,add_cbar=True,axes_loc=None,kwargs_colorbar=None,**kwargs):"""Plot effective area image. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. Default is None. add_cbar : bool, optional Add a colorbar to the plot. Default is True. axes_loc : dict, optional Keyword arguments passed to `~mpl_toolkits.axes_grid1.axes_divider.AxesDivider.append_axes`. kwargs_colorbar : dict, optional Keyword arguments passed to `~matplotlib.pyplot.colorbar`. kwargs : dict Keyword arguments passed to `~matplotlib.pyplot.pcolormesh`. Returns ------- ax : `~matplotlib.axes.Axes` Matplotlib axes. """ax=plt.gca()ifaxisNoneelseaxenergy=self.axes["energy_true"]offset=self.axes["offset"]aeff=self.evaluate(offset=offset.center,energy_true=energy.center[:,np.newaxis])vmin,vmax=np.nanmin(aeff.value),np.nanmax(aeff.value)kwargs.setdefault("cmap","GnBu")kwargs.setdefault("edgecolors","face")kwargs.setdefault("vmin",vmin)kwargs.setdefault("vmax",vmax)kwargs_colorbar=kwargs_colorbaror{}withquantity_support():caxes=ax.pcolormesh(energy.edges,offset.edges,aeff.value.T,**kwargs)energy.format_plot_xaxis(ax=ax)offset.format_plot_yaxis(ax=ax)ifadd_cbar:label=f"Effective Area [{aeff.unit.to_string(UNIT_STRING_FORMAT)}]"kwargs_colorbar.setdefault("label",label)add_colorbar(caxes,ax=ax,axes_loc=axes_loc,**kwargs_colorbar)returnax
[docs]defpeek(self,figsize=(15,5)):"""Quick-look summary plots. Parameters ---------- figsize : tuple, optional Size of the figure. Default is (15, 5). """ncols=2ifself.is_pointlikeelse3fig,axes=plt.subplots(nrows=1,ncols=ncols,figsize=figsize)self.plot(ax=axes[ncols-1])self.plot_energy_dependence(ax=axes[0])ifself.is_pointlikeisFalse:self.plot_offset_dependence(ax=axes[1])plt.tight_layout()
[docs]@classmethoddeffrom_parametrization(cls,energy_axis_true=None,instrument="HESS"):r"""Create parametrized effective area. Parametrizations of the effective areas of different Cherenkov telescopes taken from Appendix B of Abramowski et al. (2010), see https://ui.adsabs.harvard.edu/abs/2010MNRAS.402.1342A . .. math:: A_{eff}(E) = g_1 \left(\frac{E}{\mathrm{MeV}}\right)^{-g_2}\exp{\left(-\frac{g_3}{E}\right)} This method does not model the offset dependence of the effective area, but just assumes that it is constant. Parameters ---------- energy_axis_true : `MapAxis`, optional Energy binning, analytic function is evaluated at log centers. Default is None. instrument : {'HESS', 'HESS2', 'CTA'} Instrument name. Default is 'HESS'. Returns ------- aeff : `EffectiveAreaTable2D` Effective area table. """# noqa: E501# Put the parameters g in a dictionary.# Units: g1 (cm^2), g2 (), g3 (MeV)pars={"HESS":[6.85e9,0.0891,5e5],"HESS2":[2.05e9,0.0891,1e5],"CTA":[1.71e11,0.0891,1e5],}ifinstrumentnotinpars.keys():ss=f"Unknown instrument: {instrument}\n"ss+=f"Valid instruments: {list(pars.keys())}"raiseValueError(ss)ifenergy_axis_trueisNone:energy_axis_true=MapAxis.from_energy_bounds("2 GeV","200 TeV",nbin=20,per_decade=True,name="energy_true")g1,g2,g3=pars[instrument]offset_axis=MapAxis.from_edges([0.0,5.0]*u.deg,name="offset")axes=MapAxes([energy_axis_true,offset_axis])coords=axes.get_coord()energy,offset=coords["energy_true"].to_value("MeV"),coords["offset"]data=np.ones_like(offset.value)*g1*energy**(-g2)*np.exp(-g3/energy)# TODO: fake offset dependence?meta={"TELESCOP":instrument}returncls(axes=axes,data=data,unit="cm2",meta=meta)