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.

This tutorial also shows how to do a classical image analysis using the ring bakground.

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
[2]:
import yaml
import numpy as np
from scipy.stats import norm
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.convolution import Tophat2DKernel
from regions import CircleSkyRegion
from gammapy.data import DataStore
from gammapy.maps import Map, MapAxis, WcsGeom
from gammapy.cube import (
    MapMaker,
    MapDataset,
    PSFKernel,
    MapMakerRing,
    RingBackgroundEstimator,
)
from gammapy.modeling.models import SkyModel, BackgroundModel
from gammapy.modeling.models import PowerLaw
from gammapy.modeling.models import create_crab_spectral_model
from gammapy.modeling.models import SkyPointSource
from gammapy.detect import compute_lima_on_off_image
from gammapy.scripts import Analysis
from gammapy.modeling import Fit

Data access

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

[3]:
data_store = DataStore.from_file(
    "$GAMMAPY_DATA/hess-dl3-dr1/hess-dl3-dr3-with-background.fits.gz"
)
mask = data_store.obs_table["TARGET_NAME"] == "Crab"
obs_table = data_store.obs_table[mask]
observations = data_store.get_observations(obs_table["OBS_ID"])
[4]:
# pos_crab = SkyCoord.from_name('Crab')
pos_crab = SkyCoord(83.633, 22.014, unit="deg")

Maps

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

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

[5]:
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),
    binsz=0.02,
    width=(5, 5),
    coordsys="CEL",
    proj="TAN",
    axes=[energy_axis],
)
[6]:
%%time
maker = MapMaker(geom, offset_max="2.5 deg")
maps = maker.run(observations)
images = maker.run_images()
CPU times: user 3.86 s, sys: 256 ms, total: 4.12 s
Wall time: 4.13 s
[7]:
maps.keys()
[7]:
dict_keys(['counts', 'exposure', 'background'])
[8]:
images["counts"].smooth(3).plot(stretch="sqrt", vmax=2);
../_images/notebooks_hess_11_0.png

PSF

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

[9]:
from gammapy.irf import make_mean_psf

table_psf = make_mean_psf(observations, pos_crab)
[10]:
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, and a background with a flexible normalisation

[11]:
spatial_model = SkyPointSource(
    lon_0="83.6 deg", lat_0="22.0 deg", frame="icrs"
)
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)
background_model = BackgroundModel(maps["background"], norm=1.0)
background_model.parameters["tilt"].frozen = False
[12]:
%%time
dataset = MapDataset(
    model=model,
    counts=maps["counts"],
    exposure=maps["exposure"],
    background_model=background_model,
    psf=psf_kernel,
)
fit = Fit(dataset)
result = fit.run()
print(result)
OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 272
        total stat : 61231.07

CPU times: user 2.54 s, sys: 45 ms, total: 2.58 s
Wall time: 2.59 s

Best fit parameters:

[13]:
result.parameters.to_table()
[13]:
Table length=8
namevalueerrorunitminmaxfrozen
str9float64float64str14float64float64bool
lon_08.362e+012.978e-03deg-1.800e+021.800e+02False
lat_02.202e+012.789e-03deg-9.000e+019.000e+01False
index2.634e+009.310e-02nannanFalse
amplitude4.992e-114.171e-12cm-2 s-1 TeV-1nannanFalse
reference1.000e+000.000e+00TeVnannanTrue
norm2.919e+006.043e-020.000e+00nanFalse
tilt-1.993e-022.061e-02nannanFalse
reference1.000e+000.000e+00TeVnannanTrue

Parameters covariance:

[14]:
result.parameters.covariance_to_table()
[14]:
Table length=8
namelon_0lat_0indexamplitudereferencenormtilt
str9float64float64float64float64float64float64float64
lon_08.871e-06-3.309e-071.464e-06-2.359e-170.000e+001.754e-07-5.872e-08
lat_0-3.309e-077.778e-067.324e-071.038e-160.000e+00-3.875e-07-6.672e-08
index1.464e-067.324e-078.668e-033.104e-130.000e+00-2.040e-04-8.833e-05
amplitude-2.359e-171.038e-163.104e-131.740e-230.000e+00-1.315e-14-3.562e-15
reference0.000e+000.000e+000.000e+000.000e+000.000e+000.000e+000.000e+00
norm1.754e-07-3.875e-07-2.040e-04-1.315e-140.000e+003.651e-031.012e-03
tilt-5.872e-08-6.672e-08-8.833e-05-3.562e-150.000e+001.012e-034.247e-04
reference0.000e+000.000e+000.000e+000.000e+000.000e+000.000e+000.000e+00

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.

[15]:
npred = dataset.npred()
residual = Map.from_geom(maps["counts"].geom)
residual.data = maps["counts"].data - npred.data
[16]:
residual.sum_over_axes().smooth("0.1 deg").plot(
    cmap="coolwarm", vmin=-0.2, vmax=0.2, add_cbar=True
);
../_images/notebooks_hess_24_0.png

Spectrum

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.

We will use the high-level interface with the following configuration:

[17]:
config = """
model:
    components:
    - name: source
      spectral:
            parameters:
            - factor: 2.0
              frozen: false
              name: index
              unit: ''
              value: 2.6
            - factor: 1.0e-12
              frozen: false
              name: amplitude
              unit: cm-2 s-1 TeV-1
              value: 5.0e-11
            - factor: 1.0
              frozen: true
              name: reference
              unit: TeV
              value: 1.0
            type: PowerLaw
observations:
  datastore: $GAMMAPY_DATA/hess-dl3-dr1/hess-dl3-dr3-with-background.fits.gz
  filters:
  - filter_type: par_value
    variable: TARGET_NAME
    value_param: Crab
reduction:
    background:
        background_estimator: reflected
        on_region:
            center:
            - 83.633 deg
            - 22.014 deg
            frame: icrs
            radius: 0.11 deg
    containment_correction: true
    data_reducer: 1d
fit:
    fit_range:
        min: 1 TeV
        max: 10 TeV
flux:
    fp_binning:
        lo_bnd: 1
        hi_bnd: 10
        nbin: 10
        unit: TeV
        interp: log
"""
config = yaml.safe_load(config)
[18]:
%%time
analysis = Analysis(config)
analysis.get_observations()
analysis.reduce()
analysis.fit()
analysis.get_flux_points()
INFO:gammapy.scripts.analysis:Setting logging config: {'level': 'INFO'}
INFO:gammapy.scripts.analysis:Fetching observations.
INFO:gammapy.scripts.analysis:Info for OBS_ID = 23592
- Start time: 53347.91
- Pointing pos: RA 82.01 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.212 %

INFO:gammapy.scripts.analysis:Info for OBS_ID = 23523
- Start time: 53343.92
- Pointing pos: RA 83.63 deg / Dec 21.51 deg
- Observation duration: 1687.0 s
- Dead-time fraction: 6.240 %

INFO:gammapy.scripts.analysis:Info for OBS_ID = 23526
- Start time: 53343.95
- Pointing pos: RA 83.63 deg / Dec 22.51 deg
- Observation duration: 1683.0 s
- Dead-time fraction: 6.555 %

INFO:gammapy.scripts.analysis:Info for OBS_ID = 23559
- Start time: 53345.96
- Pointing pos: RA 85.25 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.398 %

INFO:gammapy.scripts.analysis:Reducing data sets.
INFO:gammapy.spectrum.reflected:Found 41 reflected regions for the Obs #23592
INFO:gammapy.spectrum.reflected:Found 12 reflected regions for the Obs #23523
INFO:gammapy.spectrum.reflected:Found 12 reflected regions for the Obs #23526
INFO:gammapy.spectrum.reflected:Found 41 reflected regions for the Obs #23559
INFO:gammapy.spectrum.extract:Running <gammapy.spectrum.extract.SpectrumExtraction object at 0x11bf50f28>
INFO:gammapy.spectrum.extract:Process observation
 Info for OBS_ID = 23592
- Start time: 53347.91
- Pointing pos: RA 82.01 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.212 %

INFO:gammapy.spectrum.extract:Update observation meta info
INFO:gammapy.spectrum.extract:Offset : 1.501568530069493 deg

INFO:gammapy.spectrum.extract:Fill events
INFO:gammapy.spectrum.extract:Extract IRFs
INFO:gammapy.spectrum.extract:Apply containment correction
/Users/adonath/github/adonath/gammapy/gammapy/spectrum/extract.py:232: RuntimeWarning: invalid value encountered in true_divide
  self.containment = new_aeff.data.data.value / self._aeff.data.data.value
INFO:gammapy.spectrum.extract:Process observation
 Info for OBS_ID = 23523
- Start time: 53343.92
- Pointing pos: RA 83.63 deg / Dec 21.51 deg
- Observation duration: 1687.0 s
- Dead-time fraction: 6.240 %

INFO:gammapy.spectrum.extract:Update observation meta info
INFO:gammapy.spectrum.extract:Offset : 0.4995557435558715 deg

INFO:gammapy.spectrum.extract:Fill events
INFO:gammapy.spectrum.extract:Extract IRFs
INFO:gammapy.spectrum.extract:Apply containment correction
INFO:gammapy.spectrum.extract:Process observation
 Info for OBS_ID = 23526
- Start time: 53343.95
- Pointing pos: RA 83.63 deg / Dec 22.51 deg
- Observation duration: 1683.0 s
- Dead-time fraction: 6.555 %

INFO:gammapy.spectrum.extract:Update observation meta info
INFO:gammapy.spectrum.extract:Offset : 0.5004444451150709 deg

INFO:gammapy.spectrum.extract:Fill events
INFO:gammapy.spectrum.extract:Extract IRFs
INFO:gammapy.spectrum.extract:Apply containment correction
INFO:gammapy.spectrum.extract:Process observation
 Info for OBS_ID = 23559
- Start time: 53345.96
- Pointing pos: RA 85.25 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.398 %

INFO:gammapy.spectrum.extract:Update observation meta info
INFO:gammapy.spectrum.extract:Offset : 1.5021898826765052 deg

INFO:gammapy.spectrum.extract:Fill events
INFO:gammapy.spectrum.extract:Extract IRFs
INFO:gammapy.spectrum.extract:Apply containment correction
INFO:gammapy.scripts.analysis:Reading model.
INFO:gammapy.scripts.analysis:PowerLaw

Parameters:

           name     value   error      unit      min max frozen
        --------- --------- ----- -------------- --- --- ------
            index 2.600e+00   nan                nan nan  False
        amplitude 5.000e-11   nan cm-2 s-1 TeV-1 nan nan  False
        reference 1.000e+00   nan            TeV nan nan   True
INFO:gammapy.scripts.analysis:Fitting data sets to model.
INFO:gammapy.scripts.analysis:OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 37
        total stat : 90.20

INFO:gammapy.scripts.analysis:Calculating flux points.
INFO:gammapy.scripts.analysis:
      e_ref               ref_flux        ...        dnde_err        is_ul
       TeV              1 / (cm2 s)       ...    1 / (cm2 s TeV)
------------------ ---------------------- ... ---------------------- -----
1.1364636663857248  9.306620013351773e-12 ... 3.3940432509161936e-12 False
  1.46779926762207 6.2076576499943984e-12 ... 2.2514420094570795e-12 False
1.7782794100389228  2.279176298773673e-12 ...  2.015595929194876e-12 False
2.1544346900318834  3.381671528979042e-12 ... 1.0047945436009904e-12 False
2.7825594022071245 2.2556265439567095e-12 ...  6.133726300311522e-13 False
3.5938136638046276 1.5045373455706862e-12 ...  3.156367662615518e-13 False
4.6415888336127775 1.0035493820028052e-12 ...  1.918514604421739e-13 False
 5.623413251903488 3.6845878027952315e-13 ...  1.615215476931349e-13 False
 6.812920690579611  5.466916129059527e-13 ... 1.0752395077167674e-13 False
  8.79922543569107  3.646516531431027e-13 ... 5.5627626797410457e-14 False
CPU times: user 3.68 s, sys: 77.3 ms, total: 3.76 s
Wall time: 3.74 s
[19]:
analysis.flux_points_dataset.model.parameters.covariance = (
    analysis.fit_result.parameters.covariance
)
print(analysis.flux_points_dataset.model)
PowerLaw

Parameters:

           name     value     error        unit      min max frozen
        --------- --------- --------- -------------- --- --- ------
            index 2.583e+00 1.142e-01                nan nan  False
        amplitude 4.424e-11 4.074e-12 cm-2 s-1 TeV-1 nan nan  False
        reference 1.000e+00 0.000e+00            TeV nan nan   True

Covariance:

           name     index   amplitude reference
        --------- --------- --------- ---------
            index 1.304e-02 3.598e-13 0.000e+00
        amplitude 3.598e-13 1.660e-23 0.000e+00
        reference 0.000e+00 0.000e+00 0.000e+00
[20]:
plt.figure(figsize=(10, 8))
crab_ref = create_crab_spectral_model("hess_pl")

dataset_fp = analysis.flux_points_dataset

plot_kwargs = {
    "energy_range": [1, 10] * u.TeV,
    "flux_unit": "erg-1 cm-2 s-1",
    "energy_power": 2,
}

model_kwargs = {"label": "1D best fit model"}
model_kwargs.update(plot_kwargs)
ax_spectrum, ax_residuals = dataset_fp.peek(model_kwargs=model_kwargs)

crab_ref.plot(ax=ax_spectrum, label="H.E.S.S. 2006 PWL", **plot_kwargs)
model.spectral_model.plot(
    ax=ax_spectrum, label="3D best fit model", **plot_kwargs
)

ax_spectrum.set_ylim(1e-11, 1e-10)
ax_spectrum.legend();
../_images/notebooks_hess_29_0.png

Classical Ring Background Analysis

In this section, we do a classical ring background analysis on the MSH 15-52 region. Let us first select these runs in the datastore

[21]:
data_store = DataStore.from_file(
    "$GAMMAPY_DATA/hess-dl3-dr1/hess-dl3-dr3-with-background.fits.gz"
)
data_sel = data_store.obs_table["TARGET_NAME"] == "MSH 15-52"
obs_table = data_store.obs_table[data_sel]
observations = data_store.get_observations(obs_table["OBS_ID"])
[22]:
# pos_msh1552 = SkyCoord.from_name('MSH15-52')
pos_msh1552 = SkyCoord(228.32, -59.08, unit="deg")

We first have to define the geometry on which we make our 2D map.

[23]:
energy_axis = MapAxis.from_edges(
    np.logspace(0, 5.0, 5), unit="TeV", name="energy", interp="log"
)
geom = WcsGeom.create(
    skydir=pos_msh1552,
    binsz=0.02,
    width=(5, 5),
    coordsys="CEL",
    proj="TAN",
    axes=[energy_axis],
)

Choose an exclusion mask.

We choose an exclusion mask on the position of MSH 1552.

[24]:
regions = CircleSkyRegion(center=pos_msh1552, radius=0.3 * u.deg)
mask = Map.from_geom(geom)
mask.data = mask.geom.region_mask([regions], inside=False)
mask.get_image_by_idx([0]).plot();
../_images/notebooks_hess_36_0.png

Now, we instantiate the ring background

[25]:
ring_bkg = RingBackgroundEstimator(r_in="0.5 deg", width="0.3 deg")

To facilitate classical image analysis, we have a special class called MapMakerRing. Here, we do the analysis over the integrated energy range. To do an analysis for each image slice of the map (eg: to investigate the significance of the source detection with energy, just call run instead of run_images

[26]:
%%time
im = MapMakerRing(
    geom=geom,
    offset_max=2.0 * u.deg,
    exclusion_mask=mask,
    background_estimator=ring_bkg,
)
images = im.run_images(observations)
INFO:gammapy.cube.make:Skipping obs_id: 20137 (partial map overlap)
INFO:gammapy.cube.make:Skipping obs_id: 20283 (partial map overlap)
INFO:gammapy.cube.make:Skipping obs_id: 20302 (partial map overlap)
INFO:gammapy.cube.make:Skipping obs_id: 20322 (partial map overlap)
INFO:gammapy.cube.make:Skipping obs_id: 20324 (partial map overlap)
INFO:gammapy.cube.make:Skipping obs_id: 20344 (partial map overlap)
INFO:gammapy.cube.make:Skipping obs_id: 20346 (partial map overlap)
INFO:gammapy.cube.make:Skipping obs_id: 20365 (partial map overlap)
INFO:gammapy.cube.make:Skipping obs_id: 20367 (partial map overlap)
CPU times: user 6.81 s, sys: 429 ms, total: 7.24 s
Wall time: 7.27 s

We will use compute_lima_on_off_image in gammapy.detect to compute significance and excess maps. A common debug plot during Ring Analysis is to make histograms of the off count rates, so we will plot that as well.

For this we will first create a Tophat2DKernel convolution kernel for the significance maps. To convert from angles to pixel scales, we use geom.pixel_scales:

[27]:
scale = geom.pixel_scales[0].to("deg")
# Using a convolution radius of 0.05 degrees
theta = 0.05 * u.deg / scale
tophat = Tophat2DKernel(theta)
tophat.normalize("peak")
[28]:
lima_maps = compute_lima_on_off_image(
    images["on"],
    images["off"],
    images["exposure_on"],
    images["exposure_off"],
    tophat,
)
[29]:
significance_map = lima_maps["significance"]
excess_map = lima_maps["excess"]
[30]:
plt.figure(figsize=(10, 10))
ax1 = plt.subplot(221, projection=significance_map.geom.wcs)
ax2 = plt.subplot(222, projection=excess_map.geom.wcs)

ax1.set_title("Significance map")
significance_map.plot(ax=ax1, add_cbar=True)

ax2.set_title("Excess map")
excess_map.plot(ax=ax2, add_cbar=True);
../_images/notebooks_hess_45_0.png
[31]:
# create a 2D mask for the images
image_mask = mask.slice_by_idx({"energy": 0})
significance_map_off = significance_map * image_mask
[32]:
significance_all = significance_map.data[np.isfinite(significance_map.data)]
significance_off = significance_map_off.data[
    np.isfinite(significance_map_off.data)
]

plt.hist(
    significance_all,
    density=True,
    alpha=0.5,
    color="red",
    label="all bins",
    bins=21,
)

plt.hist(
    significance_off,
    density=True,
    alpha=0.5,
    color="blue",
    label="off bins",
    bins=21,
)

# Now, fit the off distribution with a Gaussian
mu, std = norm.fit(significance_off)
x = np.linspace(-8, 8, 50)
p = norm.pdf(x, mu, std)
plt.plot(x, p, lw=2, color="black")
plt.legend()
plt.xlabel("Significance")
plt.yscale("log")
plt.ylim(1e-5, 1)
xmin, xmax = np.min(significance_all), np.max(significance_all)
plt.xlim(xmin, xmax)

print("Fit results: mu = {:.2f}, std = {:.2f}".format(mu, std))
Fit results: mu = -0.07, std = 1.03
../_images/notebooks_hess_47_1.png

The significance and excess maps clearly indicate a bright source at the position of MSH 1552. This is also evident from the significance distribution. The off distribution should ideally be a Gaussian with mu=0, sigma=1

[33]:
print(
    "Excess from entire map: {:.1f}".format(
        np.nansum(lima_maps["excess"].data)
    )
)
print(
    "Excess from off regions: {:.1f}".format(
        np.nansum((lima_maps["excess"] * image_mask).data)
    )
)
Excess from entire map: 4832.2
Excess from off regions: -1093.5

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.

Exercises

  • Try analysing another source, e.g. RX J1713.7−3946
  • Try another model, e.g. a Gaussian spatial shape or exponential cutoff power-law spectrum.
[ ]: