Flux point fitting#

Fit spectral models to combined Fermi-LAT and IACT flux points tables.

Prerequisites#

  • Some knowledge about retrieving information from catalogs,

see Source catalogs tutorial.

Context#

Some high level studies do not rely on reduced datasets with their IRFs but directly on higher level products such as flux points. This is not ideal because flux points already contain some hypothesis for the underlying spectral shape and the uncertainties they carry are usually simplified (e.g. symmetric gaussian errors). Yet, this is an efficient way to combine heterogeneous data.

Objective: fit spectral models to combined Fermi-LAT and IACT flux points.

Proposed approach#

Here we will load, the spectral points from Fermi-LAT and TeV catalogs and fit them with various spectral models to find the best representation of the wide band spectrum.

The central class we’re going to use for this example analysis is:

In addition we will work with the following data classes:

And the following spectral model classes:

Setup#

Let us start with the usual IPython notebook and Python imports:

from astropy import units as u

# %matplotlib inline
import matplotlib.pyplot as plt
from gammapy.catalog import CATALOG_REGISTRY
from gammapy.datasets import Datasets, FluxPointsDataset
from gammapy.modeling import Fit
from gammapy.modeling.models import (
    ExpCutoffPowerLawSpectralModel,
    LogParabolaSpectralModel,
    PowerLawSpectralModel,
    SkyModel,
)

Load spectral points#

For this analysis we choose to work with the source ‘HESS J1507-622’ and the associated Fermi-LAT sources ‘3FGL J1506.6-6219’ and ‘3FHL J1507.9-6228e’. We load the source catalogs, and then access source of interest by name:

catalog_3fgl = CATALOG_REGISTRY.get_cls("3fgl")()
catalog_3fhl = CATALOG_REGISTRY.get_cls("3fhl")()
catalog_gammacat = CATALOG_REGISTRY.get_cls("gamma-cat")()

source_fermi_3fgl = catalog_3fgl["3FGL J1506.6-6219"]
source_fermi_3fhl = catalog_3fhl["3FHL J1507.9-6228e"]
source_gammacat = catalog_gammacat["HESS J1507-622"]

The corresponding flux points data can be accessed with .flux_points attribute:

dataset_gammacat = FluxPointsDataset(data=source_gammacat.flux_points, name="gammacat")
dataset_gammacat.data.to_table(sed_type="dnde", formatted=True)

dataset_3fgl = FluxPointsDataset(data=source_fermi_3fgl.flux_points, name="3fgl")
dataset_3fgl.data.to_table(sed_type="dnde", formatted=True)

dataset_3fhl = FluxPointsDataset(data=source_fermi_3fhl.flux_points, name="3fhl")
dataset_3fhl.data.to_table(sed_type="dnde", formatted=True)
Table length=5
e_refe_mine_maxdndednde_errpdnde_errndnde_ulsqrt_tsis_ul
GeVGeVGeV1 / (cm2 GeV s)1 / (cm2 GeV s)1 / (cm2 GeV s)1 / (cm2 GeV s)
float64float64float64float64float64float64float64float32bool
14.14210.00020.0009.288e-122.343e-122.128e-12nan5.660False
31.62320.00050.0002.777e-126.572e-135.818e-13nan6.940False
86.60350.000150.0002.335e-131.055e-138.554e-14nan3.835False
273.861150.000500.0006.411e-142.697e-142.133e-14nan5.697False
1000.000500.0002000.0009.188e-214.034e-15nan8.068e-150.000True


Power Law Fit#

First we start with fitting a simple PowerLawSpectralModel.

pwl = PowerLawSpectralModel(
    index=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV"
)
model = SkyModel(spectral_model=pwl, name="j1507-pl")

After creating the model we run the fit by passing the flux_points and model objects:

Datasets
--------

Dataset 0:

  Type       : FluxPointsDataset
  Name       : gammacat
  Instrument :
  Models     : ['j1507-pl']

Dataset 1:

  Type       : FluxPointsDataset
  Name       : 3fgl
  Instrument :
  Models     : ['j1507-pl']

Dataset 2:

  Type       : FluxPointsDataset
  Name       : 3fhl
  Instrument :
  Models     : ['j1507-pl']

And print the result:

print(result_pwl)

print(model)
OptimizeResult

        backend    : minuit
        method     : migrad
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 40
        total stat : 28.29

CovarianceResult

        backend    : minuit
        method     : hesse
        success    : True
        message    : Hesse terminated successfully.

SkyModel

  Name                      : j1507-pl
  Datasets names            : None
  Spectral model type       : PowerLawSpectralModel
  Spatial  model type       :
  Temporal model type       :
  Parameters:
    index                         :      1.985   +/-    0.03
    amplitude                     :   1.28e-12   +/- 1.6e-13 1 / (cm2 s TeV)
    reference             (frozen):      1.000       TeV

Finally we plot the data points and the best fit model:

ax = plt.subplot()
ax.yaxis.set_units(u.Unit("TeV cm-2 s-1"))

kwargs = {"ax": ax, "sed_type": "e2dnde"}

for d in datasets:
    d.data.plot(label=d.name, **kwargs)

energy_bounds = [1e-4, 1e2] * u.TeV
pwl.plot(energy_bounds=energy_bounds, color="k", **kwargs)
pwl.plot_error(energy_bounds=energy_bounds, **kwargs)
ax.set_ylim(1e-13, 1e-11)
ax.set_xlim(energy_bounds)
ax.legend()
plt.show()
sed fitting

Exponential Cut-Off Powerlaw Fit#

Next we fit an ExpCutoffPowerLawSpectralModel law to the data.

ecpl = ExpCutoffPowerLawSpectralModel(
    index=1.8,
    amplitude="2e-12 cm-2 s-1 TeV-1",
    reference="1 TeV",
    lambda_="0.1 TeV-1",
)
model = SkyModel(spectral_model=ecpl, name="j1507-ecpl")

We run the fitter again by passing the flux points and the model instance:

SkyModel

  Name                      : j1507-ecpl
  Datasets names            : None
  Spectral model type       : ExpCutoffPowerLawSpectralModel
  Spatial  model type       :
  Temporal model type       :
  Parameters:
    index                         :      1.894   +/-    0.05
    amplitude                     :   1.96e-12   +/- 3.9e-13 1 / (cm2 s TeV)
    reference             (frozen):      1.000       TeV
    lambda_                       :      0.078   +/-    0.05 1 / TeV
    alpha                 (frozen):      1.000

We plot the data and best fit model:

sed fitting

Log-Parabola Fit#

Finally we try to fit a LogParabolaSpectralModel model:

log_parabola = LogParabolaSpectralModel(
    alpha=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV", beta=0.1
)
model = SkyModel(spectral_model=log_parabola, name="j1507-lp")

datasets.models = model
result_log_parabola = fitter.run(datasets=datasets)
print(model)

ax = plt.subplot()

kwargs = {"ax": ax, "sed_type": "e2dnde"}

ax.yaxis.set_units(u.Unit("TeV cm-2 s-1"))

for d in datasets:
    d.data.plot(label=d.name, **kwargs)

log_parabola.plot(energy_bounds=energy_bounds, color="k", **kwargs)
log_parabola.plot_error(energy_bounds=energy_bounds, **kwargs)
ax.set_ylim(1e-13, 1e-11)
ax.set_xlim(energy_bounds)
ax.legend()
plt.show()
sed fitting
SkyModel

  Name                      : j1507-lp
  Datasets names            : None
  Spectral model type       : LogParabolaSpectralModel
  Spatial  model type       :
  Temporal model type       :
  Parameters:
    amplitude                     :   1.88e-12   +/- 2.8e-13 1 / (cm2 s TeV)
    reference             (frozen):      1.000       TeV
    alpha                         :      2.144   +/-    0.07
    beta                          :      0.049   +/-    0.02

Exercises#

What next?#

This was an introduction to SED fitting in Gammapy.

  • If you would like to learn how to perform a full Poisson maximum likelihood spectral fit, please check out the Spectral analysis tutorial.

  • To learn how to combine heterogeneous datasets to perform a multi-instrument forward-folding fit see the Multi instrument joint 3D and 1D analysis tutorial.

Gallery generated by Sphinx-Gallery