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

H.E.S.S. data with Gammapy

H.E.S.S. is an array of gamma-ray telescopes located in Namibia. Gammapy is regularly used and fully supports H.E.S.S. high-level data analysis, after export to the current open data level 3 format.

The H.E.S.S. data is private, and H.E.S.S. analysis is mostly documented and discussed at https://hess-confluence.desy.de/ and in H.E.S.S.-internal communication channels. However, in 2018, a small sub-set of archival H.E.S.S. data was publicly released, called the H.E.S.S. DL3 DR1, the data level 3, data release number 1. This dataset is 50 MB in size and is used in many Gammapy analysis tutorials, and can be downloaded via gammapy download (TODO: add link to description).

This notebook is a quick introduction to H.E.S.S. data and instrument responses and contains some specifics that are important for H.E.S.S. users:

  • IRF formats and shapes

  • How to handle safe energy and max offset

  • EVENTS and GTI formats (e.g. how HESS 1, 2, configs, … are handled)

  • Link to HESS Confluence where data and help is available (Slack channel)?

Then at the end, link to other analysis tutorials that are likely of interest for H.E.S.S. people, and add a few exercises. This can be short, a 5-10 min read. It’s just supposed to be the “landing page” for someone new in H.E.S.S. that has never used Gammapy.

DL3 DR1

This is how to access data and IRFs from the H.E.S.S. data level 3, data release 1.

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
import numpy as np
import astropy.units as u
[2]:
from gammapy.data import DataStore
from gammapy.maps import MapAxis
from gammapy.makers.utils import make_theta_squared_table
from gammapy.visualization import plot_theta_squared_table
[3]:
data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1")
[4]:
data_store.info()
Data store:
HDU index table:
BASE_DIR: /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/hess-dl3-dr1
Rows: 630
OBS_ID: 20136 -- 47829
HDU_TYPE: ['aeff', 'bkg', 'edisp', 'events', 'gti', 'psf']
HDU_CLASS: ['aeff_2d', 'bkg_3d', 'edisp_2d', 'events', 'gti', 'psf_table']


Observation table:
Observatory name: 'N/A'
Number of observations: 105

[5]:
data_store.obs_table[:2][["OBS_ID", "DATE-OBS", "RA_PNT", "DEC_PNT", "OBJECT"]]
[5]:
ObservationTable length=2
OBS_IDDATE-OBSRA_PNTDEC_PNTOBJECT
degdeg
int64bytes10float32float32bytes18
201362004-03-26228.6125-58.771667MSH15-52
201372004-03-26228.6125-59.771667MSH15-52
[6]:
obs = data_store.obs(23523)
[7]:
obs.aeff.peek()
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/gammapy/irf/effective_area.py:547: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  caxes = ax.pcolormesh(energy.value, offset.value, aeff.value.T, **kwargs)
/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: invalid value encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
../_images/tutorials_hess_8_1.png
[8]:
obs.edisp.peek()
../_images/tutorials_hess_9_0.png
[9]:
obs.psf.peek()
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/gammapy/irf/psf_3d.py:398: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  caxes = ax.pcolormesh(x, y, containment.value.T, **kwargs)
../_images/tutorials_hess_10_1.png
[10]:
obs.bkg.to_2d().plot()
../_images/tutorials_hess_11_0.png

Theta sqared event distribution

As a quick look plot it can be helpful to plot the quadratic offset (theta squared) distribution of the events.

[11]:
position = SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs")
theta2_axis = MapAxis.from_bounds(0, 0.2, nbin=20, interp="lin", unit="deg2")

observations = data_store.get_observations([23523, 23526])
theta2_table = make_theta_squared_table(
    observations=observations,
    position=position,
    theta_squared_axis=theta2_axis,
)
[12]:
plt.figure(figsize=(10, 5))
plot_theta_squared_table(theta2_table)
../_images/tutorials_hess_14_0.png

Exercises

  • Find the OBS_ID for the runs of the Crab nebula

  • Compute the expected number of background events in the whole FOV for OBS_ID=23523 in the 1 TeV to 3 TeV energy band, from the background IRF.

Next steps

Now you know how to access and work with H.E.S.S. data. All other tutorials and documentation apply to H.E.S.S. and CTA or any other IACT that provides DL3 data and IRFs in the standard format.