# Licensed under a 3-clause BSD style license - see LICENSE.rst"""Spectral models for Gammapy."""importloggingimportoperatorimportosimportwarningsfrompathlibimportPathimportnumpyasnpimportscipy.optimizeimportscipy.specialimportastropy.unitsasufromastropyimportconstantsasconstfromastropy.tableimportTablefromastropy.utils.decoratorsimportclasspropertyfromastropy.visualizationimportquantity_supportimportmatplotlib.pyplotaspltfromgammapy.mapsimportMapAxis,RegionNDMapfromgammapy.maps.axesimportUNIT_STRING_FORMATfromgammapy.modelingimportParameter,Parametersfromgammapy.utils.deprecationimportGammapyDeprecationWarningfromgammapy.utils.integrateimporttrapz_loglogfromgammapy.utils.interpolationimport(ScaledRegularGridInterpolator,interpolation_scale,)fromgammapy.utils.rootsimportfind_rootsfromgammapy.utils.scriptsimportmake_pathfrom.coreimportModelBaselog=logging.getLogger(__name__)__all__=["BrokenPowerLawSpectralModel","CompoundSpectralModel","ConstantSpectralModel","EBLAbsorptionNormSpectralModel","ExpCutoffPowerLaw3FGLSpectralModel","ExpCutoffPowerLawNormSpectralModel","ExpCutoffPowerLawSpectralModel","GaussianSpectralModel","integrate_spectrum","LogParabolaNormSpectralModel","LogParabolaSpectralModel","NaimaSpectralModel","PiecewiseNormSpectralModel","PowerLaw2SpectralModel","PowerLawNormSpectralModel","PowerLawSpectralModel","scale_plot_flux","ScaleSpectralModel","SmoothBrokenPowerLawSpectralModel","SpectralModel","SuperExpCutoffPowerLaw3FGLSpectralModel","SuperExpCutoffPowerLaw4FGLDR3SpectralModel","SuperExpCutoffPowerLaw4FGLSpectralModel","TemplateSpectralModel","TemplateNDSpectralModel","EBL_DATA_BUILTIN",]EBL_DATA_BUILTIN={"franceschini":"$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz","dominguez":"$GAMMAPY_DATA/ebl/ebl_dominguez11.fits.gz","finke":"$GAMMAPY_DATA/ebl/frd_abs.fits.gz","franceschini17":"$GAMMAPY_DATA/ebl/ebl_franceschini_2017.fits.gz","saldana-lopez21":"$GAMMAPY_DATA/ebl/ebl_saldana-lopez_2021.fits.gz",}
[docs]defscale_plot_flux(flux,energy_power=0):"""Scale flux to plot. Parameters ---------- flux : `Map` Flux map. energy_power : int, optional Power of energy to multiply flux axis with. Default is 0. Returns ------- flux : `Map` Scaled flux map. """energy=flux.geom.get_coord(sparse=True)["energy"]try:eunit=[_for_influx.unit.basesif_.physical_type=="energy"][0]exceptIndexError:eunit=energy.unity=flux*np.power(energy,energy_power)returny.to_unit(flux.unit*eunit**energy_power)
[docs]defintegrate_spectrum(func,energy_min,energy_max,ndecade=100):"""Integrate one-dimensional function using the log-log trapezoidal rule. Internally an oversampling of the energy bins to "ndecade" is used. Parameters ---------- func : callable Function to integrate. energy_min : `~astropy.units.Quantity` Integration range minimum. energy_max : `~astropy.units.Quantity` Integration range minimum. ndecade : int, optional Number of grid points per decade used for the integration. Default is 100. """# Here we impose to duplicate the numbernum=np.maximum(np.max(ndecade*np.log10(energy_max/energy_min)),2)energy=np.geomspace(energy_min,energy_max,num=int(num),axis=-1)integral=trapz_loglog(func(energy),energy,axis=-1)returnintegral.sum(axis=0)
[docs]classSpectralModel(ModelBase):"""Spectral model base class."""_type="spectral"
@classpropertydefis_norm_spectral_model(cls):"""Whether model is a norm spectral model."""return"Norm"incls.__name__@staticmethoddef_convert_evaluate_unit(kwargs_ref,energy):kwargs={}forname,quantityinkwargs_ref.items():ifquantity.unit.physical_type=="energy":quantity=quantity.to(energy.unit)kwargs[name]=quantityreturnkwargsdef__add__(self,model):ifnotisinstance(model,SpectralModel):model=ConstantSpectralModel(const=model)returnCompoundSpectralModel(self,model,operator.add)def__mul__(self,other):ifisinstance(other,SpectralModel):returnCompoundSpectralModel(self,other,operator.mul)else:raiseTypeError(f"Multiplication invalid for type {other!r}")def__radd__(self,model):returnself.__add__(model)def__sub__(self,model):ifnotisinstance(model,SpectralModel):model=ConstantSpectralModel(const=model)returnCompoundSpectralModel(self,model,operator.sub)def__rsub__(self,model):returnself.__sub__(model)def_propagate_error(self,epsilon,fct,**kwargs):"""Evaluate error for a given function with uncertainty propagation. Parameters ---------- fct : `~astropy.units.Quantity` Function to estimate the error. epsilon : float Step size of the gradient evaluation. Given as a fraction of the parameter error. **kwargs : dict Keyword arguments. Returns ------- f_cov : `~astropy.units.Quantity` Error of the given function. """eps=np.sqrt(np.diag(self.covariance))*epsilonn,f_0=len(self.parameters),fct(**kwargs)shape=(n,len(np.atleast_1d(f_0)))df_dp=np.zeros(shape)foridx,parameterinenumerate(self.parameters):ifparameter.frozenoreps[idx]==0:continueparameter.value+=eps[idx]df=fct(**kwargs)-f_0df_dp[idx]=df.value/eps[idx]parameter.value-=eps[idx]f_cov=df_dp.T@self.covariance@df_dpf_err=np.sqrt(np.diagonal(f_cov))returnu.Quantity([np.atleast_1d(f_0.value),f_err],unit=f_0.unit).squeeze()
[docs]defevaluate_error(self,energy,epsilon=1e-4):"""Evaluate spectral model with error propagation. Parameters ---------- energy : `~astropy.units.Quantity` Energy at which to evaluate. epsilon : float, optional Step size of the gradient evaluation. Given as a fraction of the parameter error. Default is 1e-4. Returns ------- dnde, dnde_error : tuple of `~astropy.units.Quantity` Tuple of flux and flux error. """returnself._propagate_error(epsilon=epsilon,fct=self,energy=energy)
@propertydefpivot_energy(self):"""Pivot or decorrelation energy, for a given spectral model calculated numerically. It is defined as the energy at which the correlation between the spectral parameters is minimized. Returns ------- pivot energy : `~astropy.units.Quantity` The energy at which the statistical error in the computed flux is smallest. If no minimum is found, NaN will be returned. """x_unit=self.reference.unitdefmin_func(x):"""Function to minimise."""x=np.exp(x)dnde,dnde_error=self.evaluate_error(x*x_unit)returndnde_error/dndebounds=[np.log(self.reference.value)-3,np.log(self.reference.value)+3]std=np.std(min_func(x=np.linspace(bounds[0],bounds[1],100)))ifstd<1e-5:log.warning("The relative error on the flux does not depend on energy. No pivot energy found.")returnnp.nan*x_unitminimizer=scipy.optimize.minimize_scalar(min_func,bounds=bounds)ifnotminimizer.success:log.warning("No minima found in the relative error on the flux. Pivot energy computation failed.")returnnp.nan*x_unitelse:returnnp.exp(minimizer.x)*x_unit
[docs]defintegral(self,energy_min,energy_max,**kwargs):r"""Integrate spectral model numerically if no analytical solution defined. .. math:: F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} \phi(E) dE Parameters ---------- energy_min, energy_max : `~astropy.units.Quantity` Lower and upper bound of integration range. **kwargs : dict Keyword arguments passed to :func:`~gammapy.modeling.models.spectral.integrate_spectrum`. """ifhasattr(self,"evaluate_integral"):kwargs={par.name:par.quantityforparinself.parameters}kwargs=self._convert_evaluate_unit(kwargs,energy_min)returnself.evaluate_integral(energy_min,energy_max,**kwargs)else:returnintegrate_spectrum(self,energy_min,energy_max,**kwargs)
[docs]defintegral_error(self,energy_min,energy_max,epsilon=1e-4,**kwargs):"""Evaluate the error of the integral flux of a given spectrum in a given energy range. Parameters ---------- energy_min, energy_max : `~astropy.units.Quantity` Lower and upper bound of integration range. epsilon : float, optional Step size of the gradient evaluation. Given as a fraction of the parameter error. Default is 1e-4. Returns ------- flux, flux_err : tuple of `~astropy.units.Quantity` Integral flux and flux error between energy_min and energy_max. """returnself._propagate_error(epsilon=epsilon,fct=self.integral,energy_min=energy_min,energy_max=energy_max,**kwargs,)
[docs]defenergy_flux(self,energy_min,energy_max,**kwargs):r"""Compute energy flux in given energy range. .. math:: G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} E \phi(E) dE Parameters ---------- energy_min, energy_max : `~astropy.units.Quantity` Lower and upper bound of integration range. **kwargs : dict Keyword arguments passed to :func:`~gammapy.modeling.models.spectral.integrate_spectrum`. """deff(x):returnx*self(x)ifhasattr(self,"evaluate_energy_flux"):kwargs={par.name:par.quantityforparinself.parameters}kwargs=self._convert_evaluate_unit(kwargs,energy_min)returnself.evaluate_energy_flux(energy_min,energy_max,**kwargs)else:returnintegrate_spectrum(f,energy_min,energy_max,**kwargs)
[docs]defenergy_flux_error(self,energy_min,energy_max,epsilon=1e-4,**kwargs):"""Evaluate the error of the energy flux of a given spectrum in a given energy range. Parameters ---------- energy_min, energy_max : `~astropy.units.Quantity` Lower and upper bound of integration range. epsilon : float, optional Step size of the gradient evaluation. Given as a fraction of the parameter error. Default is 1e-4. Returns ------- energy_flux, energy_flux_err : tuple of `~astropy.units.Quantity` Energy flux and energy flux error between energy_min and energy_max. """returnself._propagate_error(epsilon=epsilon,fct=self.energy_flux,energy_min=energy_min,energy_max=energy_max,**kwargs,)
[docs]defreference_fluxes(self,energy_axis):"""Get reference fluxes for a given energy axis. Parameters ---------- energy_axis : `MapAxis` Energy axis. Returns ------- fluxes : dict of `~astropy.units.Quantity` Reference fluxes. """energy=energy_axis.centerenergy_min,energy_max=energy_axis.edges_min,energy_axis.edges_maxreturn{"e_ref":energy,"e_min":energy_min,"e_max":energy_max,"ref_dnde":self(energy),"ref_flux":self.integral(energy_min,energy_max),"ref_eflux":self.energy_flux(energy_min,energy_max),"ref_e2dnde":self(energy)*energy**2,}
def_get_plot_flux(self,energy,sed_type):flux=RegionNDMap.create(region=None,axes=[energy])flux_err=RegionNDMap.create(region=None,axes=[energy])ifsed_typein["dnde","norm"]:flux.quantity,flux_err.quantity=self.evaluate_error(energy.center)elifsed_type=="e2dnde":flux.quantity,flux_err.quantity=energy.center**2*self.evaluate_error(energy.center)elifsed_type=="flux":flux.quantity,flux_err.quantity=self.integral_error(energy.edges_min,energy.edges_max)elifsed_type=="eflux":flux.quantity,flux_err.quantity=self.energy_flux_error(energy.edges_min,energy.edges_max)else:raiseValueError(f"Not a valid SED type: '{sed_type}'")returnflux,flux_err
[docs]defplot(self,energy_bounds,ax=None,sed_type="dnde",energy_power=0,n_points=100,**kwargs,):"""Plot spectral model curve. By default a log-log scaling of the axes is used, if you want to change the y-axis scaling to linear you can use:: >>> from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel >>> from astropy import units as u >>> pwl = ExpCutoffPowerLawSpectralModel() >>> ax = pwl.plot(energy_bounds=(0.1, 100) * u.TeV) >>> ax.set_yscale('linear') Parameters ---------- energy_bounds : `~astropy.units.Quantity` Plot energy bounds passed to MapAxis.from_energy_bounds. ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. Default is None. sed_type : {"dnde", "flux", "eflux", "e2dnde"} Evaluation methods of the model. Default is "dnde". energy_power : int, optional Power of energy to multiply flux axis with. Default is 0. n_points : int, optional Number of evaluation nodes. Default is 100. **kwargs : dict Keyword arguments forwarded to `~matplotlib.pyplot.plot`. Returns ------- ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. """fromgammapy.estimators.map.coreimportDEFAULT_UNITax=plt.gca()ifaxisNoneelseaxifself.is_norm_spectral_model:sed_type="norm"energy_min,energy_max=energy_boundsenergy=MapAxis.from_energy_bounds(energy_min,energy_max,n_points,)ifax.yaxis.unitsisNone:ax.yaxis.set_units(DEFAULT_UNIT[sed_type]*energy.unit**energy_power)flux,_=self._get_plot_flux(sed_type=sed_type,energy=energy)flux=scale_plot_flux(flux,energy_power=energy_power)withquantity_support():ax.plot(energy.center,flux.quantity[:,0,0],**kwargs)self._plot_format_ax(ax,energy_power,sed_type)returnax
[docs]defplot_error(self,energy_bounds,ax=None,sed_type="dnde",energy_power=0,n_points=100,**kwargs,):"""Plot spectral model error band. .. note:: This method calls ``ax.set_yscale("log", nonpositive='clip')`` and ``ax.set_xscale("log", nonposx='clip')`` to create a log-log representation. The additional argument ``nonposx='clip'`` avoids artefacts in the plot, when the error band extends to negative values (see also https://github.com/matplotlib/matplotlib/issues/8623). When you call ``plt.loglog()`` or ``plt.semilogy()`` explicitly in your plotting code and the error band extends to negative values, it is not shown correctly. To circumvent this issue also use ``plt.loglog(nonposx='clip', nonpositive='clip')`` or ``plt.semilogy(nonpositive='clip')``. Parameters ---------- energy_bounds : `~astropy.units.Quantity` Plot energy bounds passed to `~gammapy.maps.MapAxis.from_energy_bounds`. ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. Default is None. sed_type : {"dnde", "flux", "eflux", "e2dnde"} Evaluation methods of the model. Default is "dnde". energy_power : int, optional Power of energy to multiply flux axis with. Default is 0. n_points : int, optional Number of evaluation nodes. Default is 100. **kwargs : dict Keyword arguments forwarded to `matplotlib.pyplot.fill_between`. Returns ------- ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. """fromgammapy.estimators.map.coreimportDEFAULT_UNITax=plt.gca()ifaxisNoneelseaxifself.is_norm_spectral_model:sed_type="norm"energy_min,energy_max=energy_boundsenergy=MapAxis.from_energy_bounds(energy_min,energy_max,n_points,)kwargs.setdefault("facecolor","black")kwargs.setdefault("alpha",0.2)kwargs.setdefault("linewidth",0)ifax.yaxis.unitsisNone:ax.yaxis.set_units(DEFAULT_UNIT[sed_type]*energy.unit**energy_power)flux,flux_err=self._get_plot_flux(sed_type=sed_type,energy=energy)y_lo=scale_plot_flux(flux-flux_err,energy_power).quantity[:,0,0]y_hi=scale_plot_flux(flux+flux_err,energy_power).quantity[:,0,0]withquantity_support():ax.fill_between(energy.center,y_lo,y_hi,**kwargs)self._plot_format_ax(ax,energy_power,sed_type)returnax
[docs]defspectral_index(self,energy,epsilon=1e-5):"""Compute spectral index at given energy. Parameters ---------- energy : `~astropy.units.Quantity` Energy at which to estimate the index. epsilon : float, optional Fractional energy increment to use for determining the spectral index. Default is 1e-5. Returns ------- index : float Estimated spectral index. """f1=self(energy)f2=self(energy*(1+epsilon))returnnp.log(f1/f2)/np.log(1+epsilon)
[docs]defspectral_index_error(self,energy,epsilon=1e-5):"""Evaluate the error on spectral index at the given energy. Parameters ---------- energy : `~astropy.units.Quantity` Energy at which to estimate the index. epsilon : float, optional Fractional energy increment to use for determining the spectral index. Default is 1e-5. Returns ------- index, index_error : tuple of float Estimated spectral index and its error. """returnself._propagate_error(epsilon=epsilon,fct=self.spectral_index,energy=energy)
[docs]definverse(self,value,energy_min=0.1*u.TeV,energy_max=100*u.TeV):"""Return energy for a given function value of the spectral model. Calls the `scipy.optimize.brentq` numerical root finding method. Parameters ---------- value : `~astropy.units.Quantity` Function value of the spectral model. energy_min : `~astropy.units.Quantity`, optional Lower energy bound of the roots finding. Default is 0.1 TeV. energy_max : `~astropy.units.Quantity`, optional Upper energy bound of the roots finding. Default is 100 TeV. Returns ------- energy : `~astropy.units.Quantity` Energies at which the model has the given ``value``. """eunit="TeV"energy_min=energy_min.to(eunit)energy_max=energy_max.to(eunit)deff(x):# scale by 1e12 to achieve better precisionenergy=u.Quantity(x,eunit,copy=False)y=self(energy).to_value(value.unit)return1e12*(y-value.value)roots,res=find_roots(f,energy_min,energy_max,points_scale="log")returnroots
[docs]definverse_all(self,values,energy_min=0.1*u.TeV,energy_max=100*u.TeV):"""Return energies for multiple function values of the spectral model. Calls the `scipy.optimize.brentq` numerical root finding method. Parameters ---------- values : `~astropy.units.Quantity` Function values of the spectral model. energy_min : `~astropy.units.Quantity`, optional Lower energy bound of the roots finding. Default is 0.1 TeV. energy_max : `~astropy.units.Quantity`, optional Upper energy bound of the roots finding. Default is 100 TeV. Returns ------- energy : list of `~astropy.units.Quantity` Each element contains the energies at which the model has corresponding value of ``values``. """energies=[]forvalinnp.atleast_1d(values):res=self.inverse(val,energy_min,energy_max)energies.append(res)returnenergies
[docs]classConstantSpectralModel(SpectralModel):r"""Constant model. For more information see :ref:`constant-spectral-model`. Parameters ---------- const : `~astropy.units.Quantity` :math:`k`. Default is 1e-12 cm-2 s-1 TeV-1. """tag=["ConstantSpectralModel","const"]const=Parameter("const","1e-12 cm-2 s-1 TeV-1")const._is_norm=True
[docs]@staticmethoddefevaluate(energy,const):"""Evaluate the model (static function)."""returnnp.ones(np.atleast_1d(energy).shape)*const
[docs]classCompoundSpectralModel(SpectralModel):"""Arithmetic combination of two spectral models. For more information see :ref:`compound-spectral-model`. """tag=["CompoundSpectralModel","compound"]def__init__(self,model1,model2,operator):self.model1=model1self.model2=model2self.operator=operatorsuper().__init__()@propertydefparameters(self):returnself.model1.parameters+self.model2.parametersdef__str__(self):return(f"{self.__class__.__name__}\n"f" Component 1 : {self.model1}\n"f" Component 2 : {self.model2}\n"f" Operator : {self.operator.__name__}\n")
[docs]defto_dict(self,full_output=False):dict1=self.model1.to_dict(full_output)dict2=self.model2.to_dict(full_output)return{self._type:{"type":self.tag[0],"model1":dict1["spectral"],# for cleaner output"model2":dict2["spectral"],"operator":self.operator.__name__,}}
[docs]classPowerLawSpectralModel(SpectralModel):r"""Spectral power-law model. For more information see :ref:`powerlaw-spectral-model`. Parameters ---------- index : `~astropy.units.Quantity` :math:`\Gamma`. Default is 2.0. amplitude : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1e-12 cm-2 s-1 TeV-1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 1 TeV. See Also -------- PowerLaw2SpectralModel, PowerLawNormSpectralModel """tag=["PowerLawSpectralModel","pl"]index=Parameter("index",2.0)amplitude=Parameter("amplitude","1e-12 cm-2 s-1 TeV-1",scale_method="scale10",interp="log",)amplitude._is_norm=Truereference=Parameter("reference","1 TeV",frozen=True)
[docs]@staticmethoddefevaluate(energy,index,amplitude,reference):"""Evaluate the model (static function)."""returnamplitude*np.power((energy/reference),-index)
[docs]@staticmethoddefevaluate_integral(energy_min,energy_max,index,amplitude,reference):r"""Integrate power law analytically (static function). .. math:: F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}\phi(E)dE = \left. \phi_0 \frac{E_0}{-\Gamma + 1} \left( \frac{E}{E_0} \right)^{-\Gamma + 1} \right \vert _{E_{min}}^{E_{max}} Parameters ---------- energy_min, energy_max : `~astropy.units.Quantity` Lower and upper bound of integration range. """val=-1*index+1prefactor=amplitude*reference/valupper=np.power((energy_max/reference),val)lower=np.power((energy_min/reference),val)integral=prefactor*(upper-lower)mask=np.isclose(val,0)ifmask.any():integral[mask]=(amplitude*reference*np.log(energy_max/energy_min))[mask]returnintegral
[docs]@staticmethoddefevaluate_energy_flux(energy_min,energy_max,index,amplitude,reference):r"""Compute energy flux in given energy range analytically (static function). .. math:: G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}E \phi(E)dE = \left. \phi_0 \frac{E_0^2}{-\Gamma + 2} \left( \frac{E}{E_0} \right)^{-\Gamma + 2} \right \vert _{E_{min}}^{E_{max}} Parameters ---------- energy_min, energy_max : `~astropy.units.Quantity` Lower and upper bound of integration range. """val=-1*index+2prefactor=amplitude*reference**2/valupper=(energy_max/reference)**vallower=(energy_min/reference)**valenergy_flux=prefactor*(upper-lower)mask=np.isclose(val,0)ifmask.any():# see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2)# for referenceenergy_flux[mask]=(amplitude*reference**2*np.log(energy_max/energy_min)[mask])returnenergy_flux
[docs]definverse(self,value,*args):"""Return energy for a given function value of the spectral model. Parameters ---------- value : `~astropy.units.Quantity` Function value of the spectral model. """base=value/self.amplitude.quantityreturnself.reference.quantity*np.power(base,-1.0/self.index.value)
@propertydefpivot_energy(self):r"""The pivot or decorrelation energy is defined as: .. math:: E_D = E_0 * \exp{cov(\phi_0, \Gamma) / (\phi_0 \Delta \Gamma^2)} Formula (1) in https://arxiv.org/pdf/0910.4881.pdf Returns ------- pivot energy : `~astropy.units.Quantity` If no minimum is found, NaN will be returned. """index_err=self.index.errorreference=self.reference.quantityamplitude=self.amplitude.quantitycov_index_ampl=self.covariance.data[0,1]*amplitude.unitreturnreference*np.exp(cov_index_ampl/(amplitude*index_err**2))
[docs]classPowerLawNormSpectralModel(SpectralModel):r"""Spectral power-law model with normalized amplitude parameter. Parameters ---------- tilt : `~astropy.units.Quantity` :math:`\Gamma`. Default is 0. norm : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 1 TeV. See Also -------- PowerLawSpectralModel, PowerLaw2SpectralModel """tag=["PowerLawNormSpectralModel","pl-norm"]norm=Parameter("norm",1,unit="",interp="log")norm._is_norm=Truetilt=Parameter("tilt",0,frozen=True)reference=Parameter("reference","1 TeV",frozen=True)
[docs]@staticmethoddefevaluate(energy,tilt,norm,reference):"""Evaluate the model (static function)."""returnnorm*np.power((energy/reference),-tilt)
[docs]@staticmethoddefevaluate_energy_flux(energy_min,energy_max,tilt,norm,reference):"""Evaluate the energy flux (static function)."""val=-1*tilt+2prefactor=norm*reference**2/valupper=(energy_max/reference)**vallower=(energy_min/reference)**valenergy_flux=prefactor*(upper-lower)mask=np.isclose(val,0)ifmask.any():# see https://www.wolframalpha.com/input/?i=a+*+x+*+(x%2Fb)+%5E+(-2)# for referenceenergy_flux[mask]=(norm*reference**2*np.log(energy_max/energy_min)[mask])returnenergy_flux
[docs]definverse(self,value,*args):"""Return energy for a given function value of the spectral model. Parameters ---------- value : `~astropy.units.Quantity` Function value of the spectral model. """base=value/self.norm.quantityreturnself.reference.quantity*np.power(base,-1.0/self.tilt.value)
@propertydefpivot_energy(self):r"""The pivot or decorrelation energy is defined as: .. math:: E_D = E_0 * \exp{cov(\phi_0, \Gamma) / (\phi_0 \Delta \Gamma^2)} Formula (1) in https://arxiv.org/pdf/0910.4881.pdf Returns ------- pivot energy : `~astropy.units.Quantity` If no minimum is found, NaN will be returned. """tilt_err=self.tilt.errorreference=self.reference.quantitynorm=self.norm.quantitycov_tilt_norm=self.covariance.data[0,1]*norm.unitreturnreference*np.exp(cov_tilt_norm/(norm*tilt_err**2))
[docs]classPowerLaw2SpectralModel(SpectralModel):r"""Spectral power-law model with integral as amplitude parameter. For more information see :ref:`powerlaw2-spectral-model`. Parameters ---------- index : `~astropy.units.Quantity` Spectral index :math:`\Gamma`. Default is 2. amplitude : `~astropy.units.Quantity` Integral flux :math:`F_0`. Default is 1e-12 cm-2 s-1. emin : `~astropy.units.Quantity` Lower energy limit :math:`E_{0, min}`. Default is 0.1 TeV. emax : `~astropy.units.Quantity` Upper energy limit :math:`E_{0, max}`. Default is 100 TeV. See Also -------- PowerLawSpectralModel, PowerLawNormSpectralModel """tag=["PowerLaw2SpectralModel","pl-2"]amplitude=Parameter(name="amplitude",value="1e-12 cm-2 s-1",scale_method="scale10",interp="log",)amplitude._is_norm=Trueindex=Parameter("index",2)emin=Parameter("emin","0.1 TeV",frozen=True)emax=Parameter("emax","100 TeV",frozen=True)
[docs]@staticmethoddefevaluate(energy,amplitude,index,emin,emax):"""Evaluate the model (static function)."""top=-index+1# to get the energies dimensionless we use a modified formulabottom=emax-emin*(emin/emax)**(-index)returnamplitude*(top/bottom)*np.power(energy/emax,-index)
[docs]definverse(self,value,*args):"""Return energy for a given function value of the spectral model. Parameters ---------- value : `~astropy.units.Quantity` Function value of the spectral model. """amplitude=self.amplitude.quantityindex=self.index.valueenergy_min=self.emin.quantityenergy_max=self.emax.quantity# to get the energies dimensionless we use a modified formulatop=-index+1bottom=energy_max-energy_min*(energy_min/energy_max)**(-index)term=(bottom/top)*(value/amplitude)returnnp.power(term.to_value(""),-1.0/index)*energy_max
[docs]classBrokenPowerLawSpectralModel(SpectralModel):r"""Spectral broken power-law model. For more information see :ref:`broken-powerlaw-spectral-model`. Parameters ---------- index1 : `~astropy.units.Quantity` :math:`\Gamma1`. Default is 2. index2 : `~astropy.units.Quantity` :math:`\Gamma2`. Default is 2. amplitude : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1e-12 cm-2 s-1 TeV-1. ebreak : `~astropy.units.Quantity` :math:`E_{break}`. Default is 1 TeV. See Also -------- SmoothBrokenPowerLawSpectralModel """tag=["BrokenPowerLawSpectralModel","bpl"]index1=Parameter("index1",2.0)index2=Parameter("index2",2.0)amplitude=Parameter(name="amplitude",value="1e-12 cm-2 s-1 TeV-1",scale_method="scale10",interp="log",)amplitude._is_norm=Trueebreak=Parameter("ebreak","1 TeV")
[docs]@staticmethoddefevaluate(energy,index1,index2,amplitude,ebreak):"""Evaluate the model (static function)."""energy=np.atleast_1d(energy)cond=energy<ebreakbpwl=amplitude*np.ones(energy.shape)bpwl[cond]*=(energy[cond]/ebreak)**(-index1)bpwl[~cond]*=(energy[~cond]/ebreak)**(-index2)returnbpwl
[docs]classSmoothBrokenPowerLawSpectralModel(SpectralModel):r"""Spectral smooth broken power-law model. For more information see :ref:`smooth-broken-powerlaw-spectral-model`. Parameters ---------- index1 : `~astropy.units.Quantity` :math:`\Gamma_1`. Default is 2. index2 : `~astropy.units.Quantity` :math:`\Gamma_2`. Default is 2. amplitude : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1e-12 cm-2 s-1 TeV-1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 1 TeV. ebreak : `~astropy.units.Quantity` :math:`E_{break}`. Default is 1 TeV. beta : `~astropy.units.Quantity` :math:`\beta`. Default is 1. See Also -------- BrokenPowerLawSpectralModel """tag=["SmoothBrokenPowerLawSpectralModel","sbpl"]index1=Parameter("index1",2.0)index2=Parameter("index2",2.0)amplitude=Parameter(name="amplitude",value="1e-12 cm-2 s-1 TeV-1",scale_method="scale10",interp="log",)amplitude._is_norm=Trueebreak=Parameter("ebreak","1 TeV")reference=Parameter("reference","1 TeV",frozen=True)beta=Parameter("beta",1,frozen=True)
[docs]@staticmethoddefevaluate(energy,index1,index2,amplitude,ebreak,reference,beta):"""Evaluate the model (static function)."""beta*=np.sign(index2-index1)pwl=amplitude*(energy/reference)**(-index1)brk=(1+(energy/ebreak)**((index2-index1)/beta))**(-beta)returnpwl*brk
[docs]classPiecewiseNormSpectralModel(SpectralModel):"""Piecewise spectral correction with a free normalization at each fixed energy nodes. For more information see :ref:`piecewise-norm-spectral`. Parameters ---------- energy : `~astropy.units.Quantity` Array of energies at which the model values are given (nodes). norms : `~numpy.ndarray` or list of `Parameter` Array with the initial norms of the model at energies ``energy``. Normalisation parameters are created for each value. Default is one at each node. interp : str Interpolation scaling in {"log", "lin"}. Default is "log". """tag=["PiecewiseNormSpectralModel","piecewise-norm"]def__init__(self,energy,norms=None,interp="log"):self._energy=energyself._interp=interpifnormsisNone:norms=np.ones(len(energy))iflen(norms)!=len(energy):raiseValueError("dimension mismatch")iflen(norms)<2:raiseValueError("Input arrays must contain at least 2 elements")parameters_list=[]ifnotisinstance(norms[0],Parameter):parameters_list+=[Parameter(f"norm_{k}",norm)fork,norminenumerate(norms)]else:parameters_list+=normsself.default_parameters=Parameters(parameters_list)super().__init__()@propertydefenergy(self):"""Energy nodes."""returnself._energy@propertydefnorms(self):"""Norm values"""returnu.Quantity([p.valueforpinself.parameters])
[docs]@classmethoddeffrom_dict(cls,data,**kwargs):"""Create model from dictionary."""data=data["spectral"]energy=u.Quantity(data["energy"]["data"],data["energy"]["unit"])parameters=Parameters.from_dict(data["parameters"])if"interp"indata:returncls.from_parameters(parameters,energy=energy,interp=data["interp"])else:returncls.from_parameters(parameters,energy=energy)
[docs]@classmethoddeffrom_parameters(cls,parameters,**kwargs):"""Create model from parameters."""returncls(norms=parameters,**kwargs)
[docs]classExpCutoffPowerLawSpectralModel(SpectralModel):r"""Spectral exponential cutoff power-law model. For more information see :ref:`exp-cutoff-powerlaw-spectral-model`. Parameters ---------- index : `~astropy.units.Quantity` :math:`\Gamma`. Default is 1.5. amplitude : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1e-12 cm-2 s-1 TeV-1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 1 TeV. lambda_ : `~astropy.units.Quantity` :math:`\lambda`. Default is 0.1 TeV-1. alpha : `~astropy.units.Quantity` :math:`\alpha`. Default is 1. See Also -------- ExpCutoffPowerLawNormSpectralModel """tag=["ExpCutoffPowerLawSpectralModel","ecpl"]index=Parameter("index",1.5)amplitude=Parameter(name="amplitude",value="1e-12 cm-2 s-1 TeV-1",scale_method="scale10",interp="log",)amplitude._is_norm=Truereference=Parameter("reference","1 TeV",frozen=True)lambda_=Parameter("lambda_","0.1 TeV-1")alpha=Parameter("alpha","1.0",frozen=True)
[docs]@staticmethoddefevaluate(energy,index,amplitude,reference,lambda_,alpha):"""Evaluate the model (static function)."""pwl=amplitude*(energy/reference)**(-index)cutoff=np.exp(-np.power(energy*lambda_,alpha))returnpwl*cutoff
@propertydefe_peak(self):r"""Spectral energy distribution peak energy (`~astropy.units.Quantity`). This is the peak in E^2 x dN/dE and is given by: .. math:: E_{Peak} = \left(\frac{2 - \Gamma}{\alpha}\right)^{1/\alpha} / \lambda """reference=self.reference.quantityindex=self.index.quantitylambda_=self.lambda_.quantityalpha=self.alpha.quantityifindex>=2orlambda_==0.0oralpha==0.0:returnnp.nan*reference.unitelse:returnnp.power((2-index)/alpha,1/alpha)/lambda_
[docs]classExpCutoffPowerLawNormSpectralModel(SpectralModel):r"""Norm spectral exponential cutoff power-law model. Parameters ---------- index : `~astropy.units.Quantity` :math:`\Gamma`. Default is 1.5. norm : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 1 TeV. lambda_ : `~astropy.units.Quantity` :math:`\lambda`. Default is 0.1 TeV-1. alpha : `~astropy.units.Quantity` :math:`\alpha`. Default is 1. See Also -------- ExpCutoffPowerLawSpectralModel """tag=["ExpCutoffPowerLawNormSpectralModel","ecpl-norm"]index=Parameter("index",1.5)norm=Parameter("norm",1,unit="",interp="log")norm._is_norm=Truereference=Parameter("reference","1 TeV",frozen=True)lambda_=Parameter("lambda_","0.1 TeV-1")alpha=Parameter("alpha","1.0",frozen=True)def__init__(self,index=None,norm=None,reference=None,lambda_=None,alpha=None,**kwargs):ifindexisNone:warnings.warn("The default index value changed from 1.5 to 0 since v1.2",GammapyDeprecationWarning,)ifnormisnotNone:kwargs.update({"norm":norm})ifindexisnotNone:kwargs.update({"index":index})ifreferenceisnotNone:kwargs.update({"reference":reference})iflambda_isnotNone:kwargs.update({"lambda_":lambda_})ifalphaisnotNone:kwargs.update({"alpha":alpha})super().__init__(**kwargs)
[docs]@staticmethoddefevaluate(energy,index,norm,reference,lambda_,alpha):"""Evaluate the model (static function)."""pwl=norm*(energy/reference)**(-index)cutoff=np.exp(-np.power(energy*lambda_,alpha))returnpwl*cutoff
[docs]classExpCutoffPowerLaw3FGLSpectralModel(SpectralModel):r"""Spectral exponential cutoff power-law model used for 3FGL. For more information see :ref:`exp-cutoff-powerlaw-3fgl-spectral-model`. Parameters ---------- index : `~astropy.units.Quantity` :math:`\Gamma`. Default is 1.5. amplitude : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1e-12 cm-2 s-1 TeV-1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 1 TeV. ecut : `~astropy.units.Quantity` :math:`E_{C}`. Default is 10 TeV. """tag=["ExpCutoffPowerLaw3FGLSpectralModel","ecpl-3fgl"]index=Parameter("index",1.5)amplitude=Parameter("amplitude","1e-12 cm-2 s-1 TeV-1",scale_method="scale10",interp="log",)amplitude._is_norm=Truereference=Parameter("reference","1 TeV",frozen=True)ecut=Parameter("ecut","10 TeV")
[docs]@staticmethoddefevaluate(energy,index,amplitude,reference,ecut):"""Evaluate the model (static function)."""pwl=amplitude*(energy/reference)**(-index)cutoff=np.exp((reference-energy)/ecut)returnpwl*cutoff
[docs]classSuperExpCutoffPowerLaw3FGLSpectralModel(SpectralModel):r"""Spectral super exponential cutoff power-law model used for 3FGL. For more information see :ref:`super-exp-cutoff-powerlaw-3fgl-spectral-model`. .. math:: \phi(E) = \phi_0 \cdot \left(\frac{E}{E_0}\right)^{-\Gamma_1} \exp \left( \left(\frac{E_0}{E_{C}} \right)^{\Gamma_2} - \left(\frac{E}{E_{C}} \right)^{\Gamma_2} \right) Parameters ---------- index_1 : `~astropy.units.Quantity` :math:`\Gamma_1`. Default is 1.5. index_2 : `~astropy.units.Quantity` :math:`\Gamma_2`. Default is 2. amplitude : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1e-12 cm-2 s-1 TeV-1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 1 TeV. ecut : `~astropy.units.Quantity` :math:`E_{C}`. Default is 10 TeV. """tag=["SuperExpCutoffPowerLaw3FGLSpectralModel","secpl-3fgl"]amplitude=Parameter("amplitude","1e-12 cm-2 s-1 TeV-1",scale_method="scale10",interp="log",)amplitude._is_norm=Truereference=Parameter("reference","1 TeV",frozen=True)ecut=Parameter("ecut","10 TeV")index_1=Parameter("index_1",1.5)index_2=Parameter("index_2",2)
[docs]@staticmethoddefevaluate(energy,amplitude,reference,ecut,index_1,index_2):"""Evaluate the model (static function)."""pwl=amplitude*(energy/reference)**(-index_1)cutoff=np.exp((reference/ecut)**index_2-(energy/ecut)**index_2)returnpwl*cutoff
[docs]classSuperExpCutoffPowerLaw4FGLSpectralModel(SpectralModel):r"""Spectral super exponential cutoff power-law model used for 4FGL-DR1 (and DR2). See equation (4) of https://arxiv.org/pdf/1902.10045.pdf For more information see :ref:`super-exp-cutoff-powerlaw-4fgl-spectral-model`. Parameters ---------- index_1 : `~astropy.units.Quantity` :math:`\Gamma_1`. Default is 1.5. index_2 : `~astropy.units.Quantity` :math:`\Gamma_2`. Default is 2. amplitude : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1e-12 cm-2 s-1 TeV-1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 1 TeV. expfactor : `~astropy.units.Quantity` :math:`a`, given as dimensionless value but internally assumes unit of :math:`E_0^{-\Gamma_2}`. Default is 1e-2. """tag=["SuperExpCutoffPowerLaw4FGLSpectralModel","secpl-4fgl"]amplitude=Parameter("amplitude","1e-12 cm-2 s-1 TeV-1",scale_method="scale10",interp="log",)amplitude._is_norm=Truereference=Parameter("reference","1 TeV",frozen=True)expfactor=Parameter("expfactor","1e-2")index_1=Parameter("index_1",1.5)index_2=Parameter("index_2",2)
[docs]@staticmethoddefevaluate(energy,amplitude,reference,expfactor,index_1,index_2):"""Evaluate the model (static function)."""pwl=amplitude*(energy/reference)**(-index_1)cutoff=np.exp(expfactor/reference.unit**index_2*(reference**index_2-energy**index_2))returnpwl*cutoff
[docs]classSuperExpCutoffPowerLaw4FGLDR3SpectralModel(SpectralModel):r"""Spectral super exponential cutoff power-law model used for 4FGL-DR3. See equations (2) and (3) of https://arxiv.org/pdf/2201.11184.pdf For more information see :ref:`super-exp-cutoff-powerlaw-4fgl-dr3-spectral-model`. Parameters ---------- index_1 : `~astropy.units.Quantity` :math:`\Gamma_1`. Default is 1.5. index_2 : `~astropy.units.Quantity` :math:`\Gamma_2`. Default is 2. amplitude : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1e-12 cm-2 s-1 TeV-1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 1 TeV. expfactor : `~astropy.units.Quantity` :math:`a`, given as dimensionless value. Default is 1e-2. """tag=["SuperExpCutoffPowerLaw4FGLDR3SpectralModel","secpl-4fgl-dr3"]amplitude=Parameter(name="amplitude",value="1e-12 cm-2 s-1 TeV-1",scale_method="scale10",interp="log",)amplitude._is_norm=Truereference=Parameter("reference","1 TeV",frozen=True)expfactor=Parameter("expfactor","1e-2")index_1=Parameter("index_1",1.5)index_2=Parameter("index_2",2)
[docs]@staticmethoddefevaluate(energy,amplitude,reference,expfactor,index_1,index_2):"""Evaluate the model (static function)."""# https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html#PLSuperExpCutoff4pwl=amplitude*(energy/reference)**(-index_1)cutoff=(energy/reference)**(expfactor/index_2)*np.exp(expfactor/index_2**2*(1-(energy/reference)**index_2))mask=np.abs(index_2*np.log(energy/reference))<1e-2ln_=np.log(energy[mask]/reference)power=expfactor*(ln_/2.0+index_2/6.0*ln_**2.0+index_2**2.0/24.0*ln_**3)cutoff[mask]=(energy[mask]/reference)**powerreturnpwl*cutoff
[docs]classLogParabolaSpectralModel(SpectralModel):r"""Spectral log parabola model. For more information see :ref:`logparabola-spectral-model`. Parameters ---------- amplitude : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1e-12 cm-2 s-1 TeV-1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 10 TeV. alpha : `~astropy.units.Quantity` :math:`\alpha`. Default is 2. beta : `~astropy.units.Quantity` :math:`\beta`. Default is 1. See Also -------- LogParabolaNormSpectralModel """tag=["LogParabolaSpectralModel","lp"]amplitude=Parameter("amplitude","1e-12 cm-2 s-1 TeV-1",scale_method="scale10",interp="log",)amplitude._is_norm=Truereference=Parameter("reference","10 TeV",frozen=True)alpha=Parameter("alpha",2)beta=Parameter("beta",1)
[docs]@classmethoddeffrom_log10(cls,amplitude,reference,alpha,beta):"""Construct from :math:`log_{10}` parametrization."""beta_=beta/np.log(10)returncls(amplitude=amplitude,reference=reference,alpha=alpha,beta=beta_)
[docs]@staticmethoddefevaluate(energy,amplitude,reference,alpha,beta):"""Evaluate the model (static function)."""xx=energy/referenceexponent=-alpha-beta*np.log(xx)returnamplitude*np.power(xx,exponent)
@propertydefe_peak(self):r"""Spectral energy distribution peak energy (`~astropy.units.Quantity`). This is the peak in E^2 x dN/dE and is given by: .. math:: E_{Peak} = E_{0} \exp{ (2 - \alpha) / (2 * \beta)} """reference=self.reference.quantityalpha=self.alpha.quantitybeta=self.beta.quantityreturnreference*np.exp((2-alpha)/(2*beta))
[docs]classLogParabolaNormSpectralModel(SpectralModel):r"""Norm spectral log parabola model. Parameters ---------- norm : `~astropy.units.Quantity` :math:`\phi_0`. Default is 1. reference : `~astropy.units.Quantity` :math:`E_0`. Default is 10 TeV. alpha : `~astropy.units.Quantity` :math:`\alpha`. Default is 0. beta : `~astropy.units.Quantity` :math:`\beta`. Default is 0. See Also -------- LogParabolaSpectralModel """tag=["LogParabolaNormSpectralModel","lp-norm"]norm=Parameter("norm",1,unit="",interp="log")norm._is_norm=Truereference=Parameter("reference","10 TeV",frozen=True)alpha=Parameter("alpha",0)beta=Parameter("beta",0)def__init__(self,norm=None,reference=None,alpha=None,beta=None,**kwargs):ifalphaisNone:warnings.warn("The default alpha value changed from 2 to 0 since v1.2",GammapyDeprecationWarning,)ifbetaisNone:warnings.warn("The default beta value changed from 1 to 0 since v1.2",GammapyDeprecationWarning,)ifnormisnotNone:kwargs.update({"norm":norm})ifbetaisnotNone:kwargs.update({"beta":beta})ifreferenceisnotNone:kwargs.update({"reference":reference})ifalphaisnotNone:kwargs.update({"alpha":alpha})super().__init__(**kwargs)
[docs]@classmethoddeffrom_log10(cls,norm,reference,alpha,beta):"""Construct from :math:`log_{10}` parametrization."""beta_=beta/np.log(10)returncls(norm=norm,reference=reference,alpha=alpha,beta=beta_)
[docs]@staticmethoddefevaluate(energy,norm,reference,alpha,beta):"""Evaluate the model (static function)."""xx=energy/referenceexponent=-alpha-beta*np.log(xx)returnnorm*np.power(xx,exponent)
[docs]classTemplateSpectralModel(SpectralModel):"""A model generated from a table of energy and value arrays. For more information see :ref:`template-spectral-model`. Parameters ---------- energy : `~astropy.units.Quantity` Array of energies at which the model values are given values : `~numpy.ndarray` Array with the values of the model at energies ``energy``. interp_kwargs : dict Interpolation keyword arguments passed to `scipy.interpolate.RegularGridInterpolator`. By default all values outside the interpolation range are set to zero. If you want to apply linear extrapolation you can pass `interp_kwargs={'fill_value': 'extrapolate', 'kind': 'linear'}`. If you want to choose the interpolation scaling applied to values, you can use `interp_kwargs={"values_scale": "log"}`. meta : dict, optional Meta information, meta['filename'] will be used for serialisation. """tag=["TemplateSpectralModel","template"]def__init__(self,energy,values,interp_kwargs=None,meta=None,):self.energy=energyself.values=u.Quantity(values,copy=False)self.meta={}ifmetaisNoneelsemetainterp_kwargs=interp_kwargsor{}interp_kwargs.setdefault("values_scale","log")interp_kwargs.setdefault("points_scale",("log",))iflen(energy)==1:interp_kwargs["method"]="nearest"self._evaluate=ScaledRegularGridInterpolator(points=(energy,),values=values,**interp_kwargs)super().__init__()
[docs]@classmethoddefread_xspec_model(cls,filename,param,**kwargs):"""Read XSPEC table model. The input is a table containing absorbed values from a XSPEC model as a function of energy. TODO: Format of the file should be described and discussed in https://gamma-astro-data-formats.readthedocs.io/en/latest/index.html Parameters ---------- filename : str File containing the XSPEC model. param : float Model parameter value. Examples -------- Fill table from an EBL model (Franceschini, 2008) >>> from gammapy.modeling.models import TemplateSpectralModel >>> filename = '$GAMMAPY_DATA/ebl/ebl_franceschini.fits.gz' >>> table_model = TemplateSpectralModel.read_xspec_model(filename=filename, param=0.3) """filename=make_path(filename)# Check if parameter value is in rangetable_param=Table.read(filename,hdu="PARAMETERS")pmin=table_param["MINIMUM"]pmax=table_param["MAXIMUM"]ifparam<pminorparam>pmax:raiseValueError(f"Out of range: param={param}, min={pmin}, max={pmax}")# Get energy valuestable_energy=Table.read(filename,hdu="ENERGIES")energy_lo=table_energy["ENERG_LO"]energy_hi=table_energy["ENERG_HI"]# set energy to log-centersenergy=np.sqrt(energy_lo*energy_hi)# Get spectrum values (no interpolation, take closest value for param)table_spectra=Table.read(filename,hdu="SPECTRA")idx=np.abs(table_spectra["PARAMVAL"]-param).argmin()values=u.Quantity(table_spectra[idx][1],"",copy=False)# no dimensionkwargs.setdefault("interp_kwargs",{"values_scale":"lin"})returncls(energy=energy,values=values,**kwargs)
[docs]defevaluate(self,energy):"""Evaluate the model (static function)."""returnself._evaluate((energy,),clip=True)
[docs]@classmethoddeffrom_region_map(cls,map,**kwargs):"""Create model from region map."""energy=map.geom.axes["energy_true"].centervalues=map.quantity[:,0,0]returncls(energy=energy,values=values,**kwargs)
[docs]classTemplateNDSpectralModel(SpectralModel):"""A model generated from a ND map where extra dimensions define the parameter space. Parameters ---------- map : `~gammapy.maps.RegionNDMap` Map template. meta : dict, optional Meta information, meta['filename'] will be used for serialisation. interp_kwargs : dict Interpolation keyword arguments passed to `gammapy.maps.Map.interp_by_pix`. Default arguments are {'method': 'linear', 'fill_value': 0}. """tag=["TemplateNDSpectralModel","templateND"]def__init__(self,map,interp_kwargs=None,meta=None,filename=None):self._map=map.copy()self.meta=dict()ifmetaisNoneelsemetaiffilenameisnotNone:filename=str(make_path(filename))self.filename=filenameparameters=[]has_energy=Falseforaxisinmap.geom.axes:ifaxis.namenotin["energy_true","energy"]:unit=axis.unitcenter=(axis.bounds[1]+axis.bounds[0])/2parameter=Parameter(name=axis.name,value=center.to_value(unit),unit=unit,scale_method="scale10",min=axis.bounds[0].to_value(unit),max=axis.bounds[-1].to_value(unit),interp=axis.interp,)parameters.append(parameter)else:has_energy|=Trueifnothas_energy:raiseValueError("Invalid map, no energy axis found")self.default_parameters=Parameters(parameters)interp_kwargs=interp_kwargsor{}interp_kwargs.setdefault("values_scale","log")self._interp_kwargs=interp_kwargssuper().__init__()@propertydefmap(self):"""Template map as a `~gammapy.maps.RegionNDMap`."""returnself._map
[docs]defwrite(self,overwrite=False):""" Write the map. Parameters ---------- overwrite: bool, optional Overwrite existing file. Default is False, which will raise a warning if the template file exists already. """ifself.filenameisNone:raiseIOError("Missing filename")elifos.path.isfile(self.filename)andnotoverwrite:log.warning("Template file already exits, and overwrite is False")else:self.map.write(self.filename)
[docs]defto_dict(self,full_output=False):"""Create dictionary for YAML serilisation."""data=super().to_dict(full_output)data["spectral"]["filename"]=self.filenamedata["spectral"]["unit"]=str(self.map.unit)returndata
[docs]classScaleSpectralModel(SpectralModel):"""Wrapper to scale another spectral model by a norm factor. Parameters ---------- model : `SpectralModel` Spectral model to wrap. norm : float Multiplicative norm factor for the model value. Default is 1. """tag=["ScaleSpectralModel","scale"]norm=Parameter("norm",1,unit="",interp="log")def__init__(self,model,norm=norm.quantity):self.model=modelself._covariance=Nonesuper().__init__(norm=norm)
[docs]classEBLAbsorptionNormSpectralModel(SpectralModel):r"""Gamma-ray absorption models. For more information see :ref:`absorption-spectral-model`. Parameters ---------- energy : `~astropy.units.Quantity` Energy node values. param : `~astropy.units.Quantity` Parameter node values. data : `~astropy.units.Quantity` Model value. redshift : float Redshift of the absorption model. Default is 0.1. alpha_norm: float Norm of the EBL model. Default is 1. interp_kwargs : dict Interpolation option passed to `~gammapy.utils.interpolation.ScaledRegularGridInterpolator`. By default the models are extrapolated outside the range. To prevent this and raise an error instead use interp_kwargs = {"extrapolate": False}. """tag=["EBLAbsorptionNormSpectralModel","ebl-norm"]alpha_norm=Parameter("alpha_norm",1.0,frozen=True)redshift=Parameter("redshift",0.1,frozen=True)def__init__(self,energy,param,data,redshift,alpha_norm,interp_kwargs=None):self.filename=None# set values log centersself.param=paramself.energy=energyself.data=u.Quantity(data,copy=False)interp_kwargs=interp_kwargsor{}interp_kwargs.setdefault("points_scale",("lin","log"))interp_kwargs.setdefault("values_scale","log")interp_kwargs.setdefault("extrapolate",True)self._evaluate_table_model=ScaledRegularGridInterpolator(points=(self.param,self.energy),values=self.data,**interp_kwargs)super().__init__(redshift=redshift,alpha_norm=alpha_norm)
[docs]@classmethoddeffrom_dict(cls,data,**kwargs):data=data["spectral"]redshift=[p["value"]forpindata["parameters"]ifp["name"]=="redshift"][0]alpha_norm=[p["value"]forpindata["parameters"]ifp["name"]=="alpha_norm"][0]if"filename"indata:ifos.path.exists(data["filename"]):returncls.read(data["filename"],redshift=redshift,alpha_norm=alpha_norm)else:forreference,filenameinEBL_DATA_BUILTIN.items():ifPath(filename).stemindata["filename"]:returncls.read_builtin(reference,redshift=redshift,alpha_norm=alpha_norm)raiseIOError(f'File {data["filename"]} not found')else:energy=u.Quantity(data["energy"]["data"],data["energy"]["unit"])param=u.Quantity(data["param"]["data"],data["param"]["unit"])values=u.Quantity(data["values"]["data"],data["values"]["unit"])returncls(energy=energy,param=param,data=values,redshift=redshift,alpha_norm=alpha_norm,)
[docs]@classmethoddefread(cls,filename,redshift=0.1,alpha_norm=1,interp_kwargs=None):"""Build object from an XSPEC model. Parameters ---------- filename : str File containing the model. redshift : float, optional Redshift of the absorption model. Default is 0.1. alpha_norm: float, optional Norm of the EBL model. Default is 1. interp_kwargs : dict, optional Interpolation option passed to `~gammapy.utils.interpolation.ScaledRegularGridInterpolator`. Default is None. """# Create EBL data arrayfilename=make_path(filename)table_param=Table.read(filename,hdu="PARAMETERS")# TODO: for some reason the table contain duplicated valuesparam,idx=np.unique(table_param[0]["VALUE"],return_index=True)# Get energy valuestable_energy=Table.read(filename,hdu="ENERGIES")energy_lo=u.Quantity(table_energy["ENERG_LO"],"keV",copy=False)# unit not stored in fileenergy_hi=u.Quantity(table_energy["ENERG_HI"],"keV",copy=False)# unit not stored in fileenergy=np.sqrt(energy_lo*energy_hi)# Get spectrum valuestable_spectra=Table.read(filename,hdu="SPECTRA")data=table_spectra["INTPSPEC"].data[idx,:]model=cls(energy=energy,param=param,data=data,redshift=redshift,alpha_norm=alpha_norm,interp_kwargs=interp_kwargs,)model.filename=filenamereturnmodel
[docs]@classmethoddefread_builtin(cls,reference="dominguez",redshift=0.1,alpha_norm=1,interp_kwargs=None):"""Read from one of the built-in absorption models. Parameters ---------- reference : {'franceschini', 'dominguez', 'finke'}, optional Name of one of the available model in gammapy-data. Default is 'dominquez'. redshift : float, optional Redshift of the absorption model. Default is 0.1. alpha_norm : float, optional Norm of the EBL model. Default is 1. interp_kwargs : dict, optional Interpolation keyword arguments. Default is None. References ---------- .. [1] Franceschini et al. (2008), "Extragalactic optical-infrared background radiation, its time evolution and the cosmic photon-photon opacity", # noqa: E501 `Link <https://ui.adsabs.harvard.edu/abs/2008A%26A...487..837F>`__ .. [2] Dominguez et al. (2011), " Extragalactic background light inferred from AEGIS galaxy-SED-type fractions" # noqa: E501 `Link <https://ui.adsabs.harvard.edu/abs/2011MNRAS.410.2556D>`__ .. [3] Finke et al. (2010), "Modeling the Extragalactic Background Light from Stars and Dust" `Link <https://ui.adsabs.harvard.edu/abs/2010ApJ...712..238F>`__ .. [4] Franceschini et al. (2017), "The extragalactic background light revisited and the cosmic photon-photon opacity" `Link <https://ui.adsabs.harvard.edu/abs/2017A%26A...603A..34F/abstract>`__ .. [5] Saldana-Lopez et al. (2021) "An observational determination of the evolving extragalactic background light from the multiwavelength HST/CANDELS survey in the Fermi and CTA era" `Link <https://ui.adsabs.harvard.edu/abs/2021MNRAS.507.5144S/abstract>`__ """returncls.read(EBL_DATA_BUILTIN[reference],redshift,alpha_norm,interp_kwargs=interp_kwargs,)
[docs]defevaluate(self,energy,redshift,alpha_norm):"""Evaluate model for energy and parameter value."""absorption=np.clip(self._evaluate_table_model((redshift,energy)),0,1)returnnp.power(absorption,alpha_norm)
[docs]classNaimaSpectralModel(SpectralModel):r"""A wrapper for Naima models. For more information see :ref:`naima-spectral-model`. Parameters ---------- radiative_model : `~naima.models.BaseRadiative` An instance of a radiative model defined in `~naima.models`. distance : `~astropy.units.Quantity`, optional Distance to the source. If set to 0, the intrinsic differential luminosity will be returned. Default is 1 kpc. seed : str or list of str, optional Seed photon field(s) to be considered for the `radiative_model` flux computation, in case of a `~naima.models.InverseCompton` model. It can be a subset of the `seed_photon_fields` list defining the `radiative_model`. Default is the whole list of photon fields. nested_models : dict Additional parameters for nested models not supplied by the radiative model, for now this is used only for synchrotron self-compton model. """tag=["NaimaSpectralModel","naima"]def__init__(self,radiative_model,distance=1.0*u.kpc,seed=None,nested_models=None,use_cache=False,):importnaimaself.radiative_model=radiative_modelself.radiative_model._memoize=use_cacheself.distance=u.Quantity(distance)self.seed=seedifnested_modelsisNone:nested_models={}self.nested_models=nested_modelsifisinstance(self.particle_distribution,naima.models.TableModel):param_names=["amplitude"]else:param_names=self.particle_distribution.param_namesparameters=[]fornameinparam_names:value=getattr(self.particle_distribution,name)parameter=Parameter(name,value)parameters.append(parameter)# In case of a synchrotron radiative model, append B to the fittable parametersif"B"inself.radiative_model.param_names:value=getattr(self.radiative_model,"B")parameter=Parameter("B",value)parameters.append(parameter)# In case of a synchrotron self compton model, append B and Rpwn to the fittable parametersifself.include_ssc:B=self.nested_models["SSC"]["B"]radius=self.nested_models["SSC"]["radius"]parameters.append(Parameter("B",B))parameters.append(Parameter("radius",radius,frozen=True))self.default_parameters=Parameters(parameters)self.ssc_energy=np.logspace(-7,9,100)*u.eVsuper().__init__()@propertydefinclude_ssc(self):"""Whether the model includes an SSC component."""importnaimais_ic_model=isinstance(self.radiative_model,naima.models.InverseCompton)returnis_ic_modeland"SSC"inself.nested_models@propertydefssc_model(self):"""Synchrotron model."""importnaimaifself.include_ssc:returnnaima.models.Synchrotron(self.particle_distribution,B=self.B.quantity,Eemax=self.radiative_model.Eemax,Eemin=self.radiative_model.Eemin,)@propertydefparticle_distribution(self):"""Particle distribution."""returnself.radiative_model.particle_distributiondef_evaluate_ssc(self,energy,):""" Compute photon density spectrum from synchrotron emission for synchrotron self-compton model, assuming uniform synchrotron emissivity inside a sphere of radius R (see Section 4.1 of Atoyan & Aharonian 1996). Based on : https://naima.readthedocs.io/en/latest/examples.html#crab-nebula-ssc-model """Lsy=self.ssc_model.flux(self.ssc_energy,distance=0*u.cm)# use distance 0 to get luminosityphn_sy=Lsy/(4*np.pi*self.radius.quantity**2*const.c)*2.24# The factor 2.24 comes from the assumption on uniform synchrotron# emissivity inside a sphereif"SSC"notinself.radiative_model.seed_photon_fields:self.radiative_model.seed_photon_fields["SSC"]={"isotropic":True,"type":"array","energy":self.ssc_energy,"photon_density":phn_sy,}else:self.radiative_model.seed_photon_fields["SSC"]["photon_density"]=phn_sydnde=self.radiative_model.flux(energy,seed=self.seed,distance=self.distance)+self.ssc_model.flux(energy,distance=self.distance)returndndedef_update_naima_parameters(self,**kwargs):"""Update Naima model parameters."""forname,valueinkwargs.items():setattr(self.particle_distribution,name,value)if"B"inself.radiative_model.param_names:self.radiative_model.B=self.B.quantity
[docs]defevaluate(self,energy,**kwargs):"""Evaluate the model. Parameters ---------- energy : `~astropy.units.Quantity` Energy to evaluate the model at. Returns ------- dnde : `~astropy.units.Quantity` Differential flux at given energy. """self._update_naima_parameters(**kwargs)ifself.include_ssc:dnde=self._evaluate_ssc(energy.flatten())elifself.seedisnotNone:dnde=self.radiative_model.flux(energy.flatten(),seed=self.seed,distance=self.distance)else:dnde=self.radiative_model.flux(energy.flatten(),distance=self.distance)dnde=dnde.reshape(energy.shape)unit=1/(energy.unit*u.cm**2*u.s)returndnde.to(unit)
[docs]defto_dict(self,full_output=True):# for full_output to True otherwise brokenreturnsuper().to_dict(full_output=True)
[docs]@classmethoddeffrom_dict(cls,data,**kwargs):raiseNotImplementedError("Currently the NaimaSpectralModel cannot be read from YAML")
[docs]@classmethoddeffrom_parameters(cls,parameters,**kwargs):raiseNotImplementedError("Currently the NaimaSpectralModel cannot be built from a list of parameters.")
[docs]classGaussianSpectralModel(SpectralModel):r"""Gaussian spectral model. For more information see :ref:`gaussian-spectral-model`. Parameters ---------- amplitude : `~astropy.units.Quantity` :math:`N_0`. Default is 1e-12 cm-2 s-1. mean : `~astropy.units.Quantity` :math:`\bar{E}`. Default is 1 TeV. sigma : `~astropy.units.Quantity` :math:`\sigma`. Default is 2 TeV. """tag=["GaussianSpectralModel","gauss"]amplitude=Parameter("amplitude",1e-12*u.Unit("cm-2 s-1"),interp="log")amplitude._is_norm=Truemean=Parameter("mean",1*u.TeV)sigma=Parameter("sigma",2*u.TeV)
[docs]defintegral(self,energy_min,energy_max,**kwargs):r"""Integrate Gaussian analytically. .. math:: F(E_{min}, E_{max}) = \frac{N_0}{2} \left[ erf(\frac{E - \bar{E}}{\sqrt{2} \sigma})\right]_{E_{min}}^{E_{max}} Parameters ---------- energy_min, energy_max : `~astropy.units.Quantity` Lower and upper bound of integration range """# noqa: E501# kwargs are passed to this function but not used# this is to get a consistent API with SpectralModel.integral()u_min=((energy_min-self.mean.quantity)/(np.sqrt(2)*self.sigma.quantity)).to_value("")u_max=((energy_max-self.mean.quantity)/(np.sqrt(2)*self.sigma.quantity)).to_value("")return(self.amplitude.quantity/2*(scipy.special.erf(u_max)-scipy.special.erf(u_min)))
[docs]defenergy_flux(self,energy_min,energy_max):r"""Compute energy flux in given energy range analytically. .. math:: G(E_{min}, E_{max}) = \frac{N_0 \sigma}{\sqrt{2*\pi}}* \left[ - \exp(\frac{E_{min}-\bar{E}}{\sqrt{2} \sigma}) \right]_{E_{min}}^{E_{max}} + \frac{N_0 * \bar{E}}{2} \left[ erf(\frac{E - \bar{E}}{\sqrt{2} \sigma}) \right]_{E_{min}}^{E_{max}} Parameters ---------- energy_min, energy_max : `~astropy.units.Quantity` Lower and upper bound of integration range. """# noqa: E501u_min=((energy_min-self.mean.quantity)/(np.sqrt(2)*self.sigma.quantity)).to_value("")u_max=((energy_max-self.mean.quantity)/(np.sqrt(2)*self.sigma.quantity)).to_value("")a=self.amplitude.quantity*self.sigma.quantity/np.sqrt(2*np.pi)b=self.amplitude.quantity*self.mean.quantity/2returna*(np.exp(-(u_min**2))-np.exp(-(u_max**2)))+b*(scipy.special.erf(u_max)-scipy.special.erf(u_min))