<jupyter_start><jupyter_text>H.E.S.S. with GammapyThis tutorial explains how to analyse [H.E.S.S.](https://www.mpi-hd.mpg.de/hfm/HESS) 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](https://www.mpi-hd.mpg.de/hfm/HESS/pages/dl3-dr1/). In this tutorial we will make an image and a spectrum. The [light_curve.ipynb](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](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.<jupyter_code>%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 gammapy.data 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<jupyter_output><empty_output><jupyter_text>Data accessTo access the data, we use the `DataStore`, and we use the ``obs_table`` to select the Crab runs.<jupyter_code>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"])
# pos_crab = SkyCoord.from_name('Crab')
pos_crab = SkyCoord(83.633, 22.014, unit="deg")<jupyter_output><empty_output><jupyter_text>MapsLet's make some 3D cubes, as well as 2D images.For the energy, we make 5 bins from 1 TeV to 10 TeV.<jupyter_code>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],
)
%%time
maker = MapMaker(geom, offset_max="2.5 deg")
maps = maker.run(observations)
images = maker.make_images()
maps.keys()
images["counts"].smooth(3).plot(stretch="sqrt", vmax=2);<jupyter_output><empty_output><jupyter_text>PSFCompute the mean PSF for these observations at the Crab position.<jupyter_code>from gammapy.irf import make_mean_psf

table_psf = make_mean_psf(observations, pos_crab)
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)<jupyter_output><empty_output><jupyter_text>Map fitLet's fit this source assuming a Gaussian spatial shape and a power-law spectral shape<jupyter_code>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)
%%time
fit = MapFit(
    model=model,
    counts=maps["counts"],
    exposure=maps["exposure"],
    background=maps["background"],
    psf=psf_kernel,
)
result = fit.run()
print(result)
print(result.model.parameters.to_table())<jupyter_output><empty_output><jupyter_text>Residual imageWe 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.<jupyter_code>npred = fit.evaluator.compute_npred()
residual = Map.from_geom(maps["counts"].geom)
residual.data = maps["counts"].data - npred
residual.sum_over_axes().smooth(3).plot(
    cmap="coolwarm", vmin=-0.5, vmax=0.5, add_cbar=True
);<jupyter_output><empty_output><jupyter_text>SpectrumWe 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.<jupyter_code>%%time
on_region = CircleSkyRegion(pos_crab, 0.11 * u.deg)
exclusion_mask = images["counts"].copy()
exclusion_mask.data = np.ones_like(exclusion_mask.data, 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)
analysis.run()
print(analysis.fit.result[0])
opts = {
    "energy_range": analysis.fit.fit_range,
    "energy_power": 2,
    "flux_unit": "erg-1 cm-2 s-1",
}
axes = analysis.spectrum_result.plot(**opts)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11315789473684211, 'fraction_numerical': 0.018245614035087718, 'mean_word_length': 5.62137049941928, 'pattern_counts': {'":': 13, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '45186', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2040, 'n_tokens_neox': 1935, 'n_words': 543, 'path': 'tutorials/hess.ipynb', 'repo_name': 'fabiopintore/gammapy', 'size': '5700'}