This is a fixed-text formatted version of a Jupyter notebook

Flux point fitting in Gammapy


In this tutorial we’re going to learn how to fit spectral models to combined Fermi-LAT and IACT flux points.

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:


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

%matplotlib inline
import numpy as np
from astropy import units as u
from gammapy.modeling.models import (
from gammapy.spectrum import FluxPointsDataset, FluxPoints
from gammapy.catalog import SOURCE_CATALOGS
from gammapy.modeling import Fit

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 = SOURCE_CATALOGS["3fgl"]()
catalog_3fhl = SOURCE_CATALOGS["3fhl"]()
catalog_gammacat = SOURCE_CATALOGS["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:

flux_points_gammacat = source_gammacat.flux_points
Table length=6
TeV1 / (cm2 s TeV)1 / (cm2 s TeV)1 / (cm2 s TeV)

In the Fermi-LAT catalogs, integral flux points are given. Currently the flux point fitter only works with differential flux points, so we apply the conversion here.

flux_points_3fgl = source_fermi_3fgl.flux_points.to_sed_type(
    sed_type="dnde", model=source_fermi_3fgl.spectral_model()
flux_points_3fhl = source_fermi_3fhl.flux_points.to_sed_type(
    sed_type="dnde", model=source_fermi_3fhl.spectral_model()

Finally we stack the flux points into a single gammapy.spectrum.FluxPoints object and drop the upper limit values, because currently we can’t handle them in the fit:

# Stack flux point tables
flux_points = FluxPoints.stack(
    [flux_points_gammacat, flux_points_3fhl, flux_points_3fgl]

t = flux_points.table
t["dnde_err"] = 0.5 * (t["dnde_errn"] + t["dnde_errp"])

# Remove upper limit points, where `dnde_errn = nan`
is_ul = np.isfinite(t["dnde_err"])
flux_points = FluxPoints(t[is_ul])
FluxPoints(sed_type='dnde', n_points=14)

Power Law Fit

First we start with fitting a simple gammapy.modeling.models.PowerLawSpectralModel.

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

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

dataset_pwl = FluxPointsDataset(model, flux_points)
fitter = Fit([dataset_pwl])
result_pwl =

And print the result:


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


   name     value   error      unit      min max frozen
--------- --------- ----- -------------- --- --- ------
    index 1.985e+00   nan                nan nan  False
amplitude 1.283e-12   nan cm-2 s-1 TeV-1 nan nan  False
reference 1.000e+00   nan            TeV nan nan   True

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

ax = flux_points.plot(energy_power=2)
pwl.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)

# assign covariance for plotting
pwl.parameters.covariance = result_pwl.parameters.covariance

pwl.plot_error(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)
ax.set_ylim(1e-13, 1e-11);

Exponential Cut-Off Powerlaw Fit

Next we fit an gammapy.modeling.models.ExpCutoffPowerLawSpectralModel law to the data.

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

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

dataset_ecpl = FluxPointsDataset(model, flux_points)
fitter = Fit([dataset_ecpl])
result_ecpl =

   name     value   error      unit      min max frozen
--------- --------- ----- -------------- --- --- ------
    index 1.894e+00   nan                nan nan  False
amplitude 1.962e-12   nan cm-2 s-1 TeV-1 nan nan  False
reference 1.000e+00   nan            TeV nan nan   True
  lambda_ 7.766e-02   nan          TeV-1 nan nan  False
    alpha 1.000e+00   nan                nan nan   True

We plot the data and best fit model:

ax = flux_points.plot(energy_power=2)
ecpl.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)

# assign covariance for plotting
ecpl.parameters.covariance = result_ecpl.parameters.covariance

ecpl.plot_error(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)
ax.set_ylim(1e-13, 1e-11)
(1e-13, 1e-11)

Log-Parabola Fit

Finally we try to fit a gammapy.modeling.models.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)
dataset_log_parabola = FluxPointsDataset(model, flux_points)
fitter = Fit([dataset_log_parabola])
result_log_parabola =

   name     value   error      unit      min max frozen
--------- --------- ----- -------------- --- --- ------
amplitude 1.877e-12   nan cm-2 s-1 TeV-1 nan nan  False
reference 1.000e+00   nan            TeV nan nan   True
    alpha 2.144e+00   nan                nan nan  False
     beta 4.936e-02   nan                nan nan  False
ax = flux_points.plot(energy_power=2)
log_parabola.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)

# assign covariance for plotting
log_parabola.parameters.covariance = result_log_parabola.parameters.covariance

    energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2
ax.set_ylim(1e-13, 1e-11);


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 spectrum analysis tutorial.

  • To learn more about other parts of Gammapy (e.g. Fermi-LAT and TeV data analysis), check out the other tutorial notebooks.

  • To see what’s available in Gammapy, browse the Gammapy docs or use the full-text search.

  • If you have any questions, ask on the mailing list .

[ ]: