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

H.E.S.S. with Gammapy

This tutorial explains how to analyse H.E.S.S. data with Gammapy.

We will analyse four observation runs of the Crab nebula, which are part of the H.E.S.S. first public test data release. In this tutorial we will make an image and a spectrum. The light_curve.ipynb notbook contains an example how to make a light curve.

To do a 3D analysis, one needs to do a 3D background estimate. In background_model.ipynb we have started to make a background model, and in this notebook we have a first look at a 3D analysis. But the results aren’t OK yet, the background model needs to be improved. In this analysis, we also don’t use the energy dispersion IRF yet, and we only analyse the data in the 1 TeV to 10 TeV range. The H.E.S.S. data was only released very recently, and 3D analysis in Gammapy is new. This tutorial will be improved soon.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord
from regions import CircleSkyRegion
from import DataStore
from gammapy.maps import Map, MapAxis, WcsGeom, WcsNDMap
from gammapy.cube import MapMaker, MapFit, PSFKernel
from gammapy.cube.models import SkyModel
from gammapy.spectrum.models import PowerLaw, ExponentialCutoffPowerLaw
from gammapy.image.models import SkyGaussian, SkyPointSource
from gammapy.detect import TSMapEstimator
from gammapy.scripts import SpectrumAnalysisIACT

Data access

To access the data, we use the DataStore, and we use the obs_table to select the Crab runs.

data_store = DataStore.from_file(
mask = data_store.obs_table["TARGET_NAME"] == "Crab"
obs_table = data_store.obs_table[mask]
observations = data_store.get_observations(obs_table["OBS_ID"])
# pos_crab = SkyCoord.from_name('Crab')
pos_crab = SkyCoord(83.633, 22.014, unit="deg")


Let’s make some 3D cubes, as well as 2D images.

For the energy, we make 5 bins from 1 TeV to 10 TeV.

energy_axis = MapAxis.from_edges(
    np.logspace(0, 1.0, 5), unit="TeV", name="energy", interp="log"
geom = WcsGeom.create(
    skydir=(83.633, 22.014),
    width=(5, 5),
maker = MapMaker(geom, offset_max="2.5 deg")
maps =
images = maker.make_images()
CPU times: user 3.27 s, sys: 305 ms, total: 3.58 s
Wall time: 3.58 s
dict_keys(['counts', 'exposure', 'background'])
images["counts"].smooth(3).plot(stretch="sqrt", vmax=2);


Compute the mean PSF for these observations at the Crab position.

from gammapy.irf import make_mean_psf

table_psf = make_mean_psf(observations, pos_crab)
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/ RuntimeWarning: invalid value encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
psf_kernel = PSFKernel.from_table_psf(table_psf, geom, max_radius="0.3 deg")
psf_kernel_array = psf_kernel.psf_kernel_map.sum_over_axes().data
# psf_kernel.psf_kernel_map.slice_by_idx({'energy': 0}).plot()
# plt.imshow(psf_kernel_array)

Map fit

Let’s fit this source assuming a Gaussian spatial shape and a power-law spectral shape

spatial_model = SkyPointSource(lon_0="83.6 deg", lat_0="22.0 deg")
spectral_model = PowerLaw(
    index=2.6, amplitude="5e-11 cm-2 s-1 TeV-1", reference="1 TeV"
model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)
fit = MapFit(
result =

        backend    : minuit
        method     : minuit
        success    : True
        nfev       : 116
        total stat : 67039.61
        message    : Optimization terminated successfully.

   name     value   error      unit      min max frozen
--------- --------- ----- -------------- --- --- ------
    lon_0 8.362e+01   nan            deg nan nan  False
    lat_0 2.203e+01   nan            deg nan nan  False
    index 2.616e+00   nan                nan nan  False
amplitude 6.012e-11   nan cm-2 s-1 TeV-1 nan nan  False
reference 1.000e+00   nan            TeV nan nan   True
CPU times: user 3.63 s, sys: 65.4 ms, total: 3.69 s
Wall time: 3.7 s

Residual image

We compute a residual image as residual = counts - model. Note that this is counts per pixel and our pixel size is 0.02 deg. Smoothing is counts-preserving. The residual image shows that currently both the source and the background modeling isn’t very good. The background model is underestimated (so residual is positive), and the source model is overestimated.

npred = fit.evaluator.compute_npred()
residual = Map.from_geom(maps["counts"].geom) = maps["counts"].data - npred
    cmap="coolwarm", vmin=-0.5, vmax=0.5, add_cbar=True


We could try to improve the background modeling and spatial model of the source. But let’s instead turn to one of the classic IACT analysis techniques: use a circular on region and reflected regions for background estimation, and derive a spectrum for the source without having to assume a spatial model, or without needing a 3D background model.

on_region = CircleSkyRegion(pos_crab, 0.11 * u.deg)
exclusion_mask = images["counts"].copy() = np.ones_like(, dtype=bool)

model = PowerLaw(
    index=2.6, amplitude="5e-11 cm-2 s-1 TeV-1", reference="1 TeV"

config = {
    "outdir": ".",
    "background": {"on_region": on_region, "exclusion_mask": exclusion_mask},
    "extraction": {"containment_correction": True},
    "fit": {"model": model, "fit_range": [1, 10] * u.TeV},
    "fp_binning": np.logspace(0, 1, 7) * u.TeV,
analysis = SpectrumAnalysisIACT(observations=observations, config=config)
CPU times: user 4.08 s, sys: 40.6 ms, total: 4.12 s
Wall time: 4.08 s

Fit result info
Model: PowerLaw


           name     value     error        unit      min max frozen
        --------- --------- --------- -------------- --- --- ------
            index 2.577e+00 1.078e-01                nan nan  False
        amplitude 4.354e-11 3.959e-12 cm-2 s-1 TeV-1 nan nan  False
        reference 1.000e+00 0.000e+00            TeV nan nan   True


           name     index   amplitude reference
        --------- --------- --------- ---------
            index 1.163e-02 3.288e-13 0.000e+00
        amplitude 3.288e-13 1.567e-23 0.000e+00
        reference 0.000e+00 0.000e+00 0.000e+00

Statistic: 27.477 (wstat)
Fit Range: [ 1. 10.] TeV

opts = {
    "energy_power": 2,
    "flux_unit": "erg-1 cm-2 s-1",
axes = analysis.spectrum_result.plot(**opts)

Again: please note that this tutorial notebook was put together quickly, the results obtained here are very preliminary. We will work on Gammapy and the analysis of data from the H.E.S.S. test release and update this tutorial soon.


  • Try analysing another source, e.g. MSH 15-52.
  • Try another model, e.g. a Gaussian spatial shape or exponential cutoff power-law spectrum.