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
import matplotlib.pyplot as plt
from astropy import units as u
from gammapy.spectrum.models import (
PowerLaw,
ExponentialCutoffPowerLaw,
LogParabola,
)
from gammapy.spectrum import FluxPointsDataset, FluxPoints
from gammapy.catalog import (
SourceCatalog3FGL,
SourceCatalogGammaCat,
SourceCatalog3FHL,
)
from gammapy.utils.fitting import Fit
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:
fermi_3fgl = SourceCatalog3FGL()
fermi_3fhl = SourceCatalog3FHL()
gammacat = SourceCatalogGammaCat("$GAMMAPY_DATA/gamma-cat/gammacat.fits.gz")
source_gammacat = gammacat["HESS J1507-622"]
source_fermi_3fgl = fermi_3fgl["3FGL J1506.6-6219"]
source_fermi_3fhl = fermi_3fhl["3FHL J1507.9-6228e"]
The corresponding flux points data can be accessed with .flux_points
attribute:
flux_points_gammacat = source_gammacat.flux_points
flux_points_gammacat.table
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 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]
)
# drop the flux upper limit values
flux_points = flux_points.drop_ul()
pwl = PowerLaw(index=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV")
After creating the model we run the fit by passing the 'flux_points'
and 'pwl'
objects:
dataset_pwl = FluxPointsDataset(pwl, flux_points, likelihood="chi2assym")
fitter = Fit(dataset_pwl)
result_pwl = fitter.run()
And print the result:
print(result_pwl)
print(pwl)
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);
Next we fit an exponential cut-off power law to the data.
ecpl = ExponentialCutoffPowerLaw(
index=1.8,
amplitude="2e-12 cm-2 s-1 TeV-1",
reference="1 TeV",
lambda_="0.1 TeV-1",
)
We run the fitter again by passing the flux points and the ecpl
model instance:
dataset_ecpl = FluxPointsDataset(ecpl, flux_points, likelihood="chi2assym")
fitter = Fit(dataset_ecpl)
result_ecpl = fitter.run()
print(ecpl)
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)
Finally we try to fit a log-parabola model:
log_parabola = LogParabola(
alpha=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV", beta=0.1
)
dataset_log_parabola = FluxPointsDataset(log_parabola, flux_points, likelihood="chi2assym")
fitter = Fit(dataset_log_parabola)
result_log_parabola = fitter.run()
print(log_parabola)
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
log_parabola.plot_error(
energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2
)
ax.set_ylim(1e-13, 1e-11);
PowerLaw2
and ExponentialCutoffPowerLaw3FGL
to the same data.ExponentialCutoffPowerLaw
model to Vela X ('HESS J0835-455') only and check if the best fit values correspond to the values given in the Gammacat catalogThis was an introduction to SED fitting in Gammapy.