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

Spectrum analysis with Gammapy (run pipeline)

In this tutorial we will learn how to perform a 1d spectral analysis.

We will use a “pipeline” or “workflow” class to run a standard analysis. If you’re interested in implementation detail of the analysis in order to create a custom analysis class, you should read (spectrum_analysis.ipynb) that executes the analysis using lower-level classes and methods in Gammapy.

In this tutorial we will use the folling Gammapy classes:

We use 4 Crab observations from H.E.S.S. for testing.


As usual, we’ll start with some setup for the notebook, and import the functionality we need.

%matplotlib inline
import matplotlib.pyplot as plt

import astropy.units as u
from astropy.coordinates import SkyCoord
from regions import CircleSkyRegion

from import EnergyBounds
from import DataStore
from gammapy.scripts import SpectrumAnalysisIACT
from gammapy.catalog import SourceCatalogGammaCat
from gammapy.maps import Map
from gammapy.spectrum.models import LogParabola
from gammapy.spectrum import CrabSpectrum

Select data

First, we select and load some H.E.S.S. data (simulated events for now). In real life you would do something fancy here, or just use the list of observations someone send you (and hope they have done something fancy before). We’ll just use the standard gammapy 4 crab runs.

data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/")
mask = data_store.obs_table["TARGET_NAME"] == "Crab"
obs_ids = data_store.obs_table["OBS_ID"][mask].data
observations = data_store.get_observations(obs_ids)
[23523 23526 23559 23592]

Configure the analysis

Now we’ll define the input for the spectrum analysis. It will be done the python way, i.e. by creating a config dict containing python objects. We plan to add also the convenience to configure the analysis using a plain text config file.

crab_pos = SkyCoord(83.633, 22.014, unit="deg")
on_region = CircleSkyRegion(crab_pos, 0.15 * u.deg)

model = LogParabola(
    amplitude=1e-11 * u.Unit("cm-2 s-1 TeV-1"),
    reference=1 * u.TeV,

flux_point_binning = EnergyBounds.equal_log_spacing(0.7, 30, 5, u.TeV)
exclusion_mask = Map.create(skydir=crab_pos, width=(10, 10), binsz=0.02)

gammacat = SourceCatalogGammaCat("$GAMMAPY_DATA/gamma-cat/gammacat.fits.gz")

regions = []
for source in gammacat:
    if not exclusion_mask.geom.contains(source.position):
    region = CircleSkyRegion(source.position, 0.15 * u.deg)
    regions.append(region) = exclusion_mask.geom.region_mask(regions, inside=False)
config = dict(
        min_distance=0.1 * u.rad,
        model=model, forward_folded=True, fit_range=flux_point_binning[[0, -1]]

Run the analysis

TODO: Clean up the log (partly done, get rid of remaining useless warnings)

analysis = SpectrumAnalysisIACT(observations=observations, config=config){"print_level": 1})

FCN = 108.59976656320704 TOTAL NCALL = 118 NCALLS = 118
EDM = 3.120012954320802e-06 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed?
0 par_000_amplitude 3.33756 0.222387 No
1 par_001_alpha 2.31962 0.193135 No
2 par_002_beta 18.9149 9.97232 No


Let’s look at the results, and also compare with a previously published Crab nebula spectrum for reference.


        backend    : minuit
        method     : minuit
        success    : True
        nfev       : 118
        total stat : 108.60
        message    : Optimization terminated successfully.

plt.figure(figsize=(8, 6))
opts = {
    "energy_range": config["fit"]["fit_range"],
    "energy_power": 2,
    "flux_unit": "erg-1 cm-2 s-1",

ax_spectrum, ax_residual = analysis.spectrum_result.peek()

CrabSpectrum().model.plot(ax=ax_spectrum, label="Crab reference", **opts)
ax_spectrum.set_ylim(1e-12, 1e-10)
<matplotlib.legend.Legend at 0x11ca0f940>


Rerun the analysis, changing some aspects of the analysis as you like:

  • only use one or two observations

  • a different spectral model

  • different config options for the spectral analysis

  • different energy binning for the spectral point computation

Observe how the measured spectrum changes.

[ ]: