{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# H.E.S.S. with Gammapy\n",
    "\n",
    "This tutorial explains how to analyse [H.E.S.S.](https://www.mpi-hd.mpg.de/hfm/HESS) data with Gammapy.\n",
    "\n",
    "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.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from regions import CircleSkyRegion\n",
    "from gammapy.data import DataStore\n",
    "from gammapy.maps import Map, MapAxis, WcsGeom, WcsNDMap\n",
    "from gammapy.cube import MapMaker, MapFit, PSFKernel\n",
    "from gammapy.cube.models import SkyModel\n",
    "from gammapy.spectrum.models import PowerLaw, ExponentialCutoffPowerLaw\n",
    "from gammapy.image.models import SkyGaussian, SkyPointSource\n",
    "from gammapy.detect import TSMapEstimator\n",
    "from gammapy.scripts import SpectrumAnalysisIACT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data access\n",
    "\n",
    "To access the data, we use the `DataStore`, and we use the ``obs_table`` to select the Crab runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store = DataStore.from_file(\n",
    "    \"$GAMMAPY_DATA/hess-dl3-dr1/hess-dl3-dr3-with-background.fits.gz\"\n",
    ")\n",
    "mask = data_store.obs_table[\"TARGET_NAME\"] == \"Crab\"\n",
    "obs_table = data_store.obs_table[mask]\n",
    "observations = data_store.get_observations(obs_table[\"OBS_ID\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pos_crab = SkyCoord.from_name('Crab')\n",
    "pos_crab = SkyCoord(83.633, 22.014, unit=\"deg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maps\n",
    "\n",
    "Let's make some 3D cubes, as well as 2D images.\n",
    "\n",
    "For the energy, we make 5 bins from 1 TeV to 10 TeV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_axis = MapAxis.from_edges(\n",
    "    np.logspace(0, 1.0, 5), unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
    ")\n",
    "geom = WcsGeom.create(\n",
    "    skydir=(83.633, 22.014),\n",
    "    binsz=0.02,\n",
    "    width=(5, 5),\n",
    "    coordsys=\"CEL\",\n",
    "    proj=\"TAN\",\n",
    "    axes=[energy_axis],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "maker = MapMaker(geom, offset_max=\"2.5 deg\")\n",
    "maps = maker.run(observations)\n",
    "images = maker.make_images()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maps.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images[\"counts\"].smooth(3).plot(stretch=\"sqrt\", vmax=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PSF\n",
    "\n",
    "Compute the mean PSF for these observations at the Crab position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.irf import make_mean_psf\n",
    "\n",
    "table_psf = make_mean_psf(observations, pos_crab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psf_kernel = PSFKernel.from_table_psf(table_psf, geom, max_radius=\"0.3 deg\")\n",
    "psf_kernel_array = psf_kernel.psf_kernel_map.sum_over_axes().data\n",
    "# psf_kernel.psf_kernel_map.slice_by_idx({'energy': 0}).plot()\n",
    "# plt.imshow(psf_kernel_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map fit\n",
    "\n",
    "Let's fit this source assuming a Gaussian spatial shape and a power-law spectral shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_model = SkyPointSource(lon_0=\"83.6 deg\", lat_0=\"22.0 deg\")\n",
    "spectral_model = PowerLaw(\n",
    "    index=2.6, amplitude=\"5e-11 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
    ")\n",
    "model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "fit = MapFit(\n",
    "    model=model,\n",
    "    counts=maps[\"counts\"],\n",
    "    exposure=maps[\"exposure\"],\n",
    "    background=maps[\"background\"],\n",
    "    psf=psf_kernel,\n",
    ")\n",
    "result = fit.run()\n",
    "print(result)\n",
    "print(result.model.parameters.to_table())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Residual image\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npred = fit.evaluator.compute_npred()\n",
    "residual = Map.from_geom(maps[\"counts\"].geom)\n",
    "residual.data = maps[\"counts\"].data - npred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "residual.sum_over_axes().smooth(3).plot(\n",
    "    cmap=\"coolwarm\", vmin=-0.5, vmax=0.5, add_cbar=True\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectrum\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "on_region = CircleSkyRegion(pos_crab, 0.11 * u.deg)\n",
    "exclusion_mask = images[\"counts\"].copy()\n",
    "exclusion_mask.data = np.ones_like(exclusion_mask.data, dtype=bool)\n",
    "\n",
    "model = PowerLaw(\n",
    "    index=2.6, amplitude=\"5e-11 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
    ")\n",
    "\n",
    "config = {\n",
    "    \"outdir\": \".\",\n",
    "    \"background\": {\"on_region\": on_region, \"exclusion_mask\": exclusion_mask},\n",
    "    \"extraction\": {\"containment_correction\": True},\n",
    "    \"fit\": {\"model\": model, \"fit_range\": [1, 10] * u.TeV},\n",
    "    \"fp_binning\": np.logspace(0, 1, 7) * u.TeV,\n",
    "}\n",
    "analysis = SpectrumAnalysisIACT(observations=observations, config=config)\n",
    "analysis.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(analysis.fit.result[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opts = {\n",
    "    \"energy_range\": analysis.fit.fit_range,\n",
    "    \"energy_power\": 2,\n",
    "    \"flux_unit\": \"erg-1 cm-2 s-1\",\n",
    "}\n",
    "axes = analysis.spectrum_result.plot(**opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Try analysing another source, e.g. MSH 15-52.\n",
    "- Try another model, e.g. a Gaussian spatial shape or exponential cutoff power-law spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}