error plot spectrum model with 4FGL

[jlzhang001]

wrong spectrum curve:

script use :
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from gammapy.catalog import SourceCatalog4FGL

catalog_4fgl = SourceCatalog4FGL()
source = catalog_4fgl["4FGL J0534.5+2200"]
flux_points = source.flux_points
flux_points.to_table(sed_type="flux")

model = source.sky_model()
print(model)
print(model.spectral_model)
print(model.spatial_model)
energy_bounds = (50 * u.MeV, 1000 * u.GeV)
opts = dict(sed_type="e2dnde", yunits=u.Unit("TeV cm-2 s-1"))
model.spectral_model.plot(energy_bounds, **opts)
model.spectral_model.plot_error(energy_bounds, **opts);
flux_points.plot(sed_type="e2dnde");
plt.show()

output info:
~/anaconda3/envs/pipthreeML/lib/python3.7/site-packages/gammapy/catalog/fermi.py:429: UserWarning: Warning: converting a masked element to nan.
"index_2": np.nan_to_num(float(self.data["Unc_PLEC_Exp_Index"])),
~/anaconda3/envs/pipthreeML/lib/python3.7/site-packages/gammapy/catalog/fermi.py:429: UserWarning: Warning: converting a masked element to nan.
"index_2": np.nan_to_num(float(self.data["Unc_PLEC_Exp_Index"])),
SkyModel

Name : 4FGL J0534.5+2200
Datasets names : None
Spectral model type : SuperExpCutoffPowerLaw4FGLSpectralModel
Spatial model type : PointSpatialModel
Temporal model type :
Parameters:
amplitude : 7.77e-11 +/- 5.0e-13 1 / (cm2 MeV s)
reference (frozen) : 1572.767 MeV
expfactor : 0.004 +/- 0.00
index_1 : 1.883 +/- 0.03
index_2 : 0.667 +/- 0.00
lon_0 : 83.637 +/- 0.00 deg
lat_0 : 22.015 +/- 0.00 deg

SuperExpCutoffPowerLaw4FGLSpectralModel

type name value unit error min max frozen link

---

spectral amplitude 7.7732e-11 cm-2 MeV-1 s-1 4.981e-13 nan nan False
spectral reference 1.5728e+03 MeV 0.000e+00 nan nan True
spectral expfactor 4.4347e-03 2.108e-04 nan nan False
spectral index_1 1.8833e+00 2.578e-02 nan nan False
spectral index_2 6.6667e-01 0.000e+00 nan nan False
PointSpatialModel

type name value unit error min max frozen link

---

spatial lon_0 8.3637e+01 deg 3.173e-03 nan nan False
spatial lat_0 2.2015e+01 deg 2.982e-03 -9.000e+01 9.000e+01 False

[adonath]

Thanks @jlzhang001! From a first look I can't see any incorrect behavior in Gammapy. It seems this is just a failed fit / bad measurement from the Fermi-LAT catalog. The warnings are from the errors values, which are masked / missing. This again points to a bad measurement. I'm not sure we can do anything about this in Gammapy...

[jlzhang001]

Thanks for your reply @adonath! I check the spectrum fit model at https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermilpsc.html, it fit the flux points well, see below plot:

I also check the model info of this source with gammapy, see the ouput info above. It is correct, the same values as in 4FGL catalog. But when plot the model, it will not fit the flux points well.

[adonath]

Thanks @jlzhang001, indeed this looks much more consistent. Maybe @QRemy can help?

[QRemy]

This plot is for the DR3 for now gammapy only support DR1-2 as the parametrization for the PLSuperExpCutoff changed in DR3.
The shape for DR2 is implemented in the SuperExpCutoffPowerLaw4FGLSpectralModel as described by equation (4) of https://arxiv.org/pdf/1902.10045.pdf (it seems correctly done to me). There was convergence issues with this parametrisation so it has been replaced by the one given by equation 2 and 3 of https://arxiv.org/pdf/2201.11184.pdf

Updating the SourceCatalog4FGL and SuperExpCutoffPowerLaw4FGLSpectralModel classes to support the latest version of the catalogue should fix the issue.

[jlzhang001]

Let's look forward to the updating!