TeV_GammaRay_Data_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PHYS 3009: TeV Gamma-Ray Data Analysis with GammaPy\n",
    "\n",
    "This jupyter notebook presents a quick analysis of H.E.S.S. observations of the Crab Nebula. The data set is part of the first public test data release. You can find more information here: https://www.mpi-hd.mpg.de/hfm/HESS/pages/dl3-dr1/\n",
    "\n",
    "The analysis is performed using gammapy, a community-developed, open-source Python package for gamma-ray astronomy (https://docs.gammapy.org/). More information on the individual data analysis steps can be found in the gammapy tutorials: https://docs.gammapy.org/1.1/tutorials/index.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports\n",
    "Let's start with importing the modules that we used in the last notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from numpy import sqrt\n",
    "\n",
    "import astropy.units as u\n",
    "\n",
    "from astropy.coordinates import (\n",
    "    SkyCoord, \n",
    "    Angle,\n",
    ")\n",
    "\n",
    "\n",
    "from gammapy.utils.check import check_tutorials_setup\n",
    "\n",
    "from gammapy.data import (\n",
    "    DataStore,\n",
    "    EventList,\n",
    ")\n",
    "\n",
    "from gammapy.stats import WStatCountsStatistic\n",
    "\n",
    "from gammapy.maps import Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data download\n",
    "The following lines check the setup and download the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "check_tutorials_setup()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import os\n",
    "#os.environ['GAMMAPY_DATA'] = './gammapy-data'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/1.1/hess-dl3-dr1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected output:\n",
    "\n",
    "```\n",
    "Data store:\n",
    "HDU index table:\n",
    "BASE_DIR: gammapy-data/1.1/hess-dl3-dr1\n",
    "Rows: 630\n",
    "OBS_ID: 20136 -- 47829\n",
    "HDU_TYPE: ['aeff', 'bkg', 'edisp', 'events', 'gti', 'psf']\n",
    "HDU_CLASS: ['aeff_2d', 'bkg_3d', 'edisp_2d', 'events', 'gti', 'psf_table']\n",
    "\n",
    "\n",
    "Observation table:\n",
    "Observatory name: 'N/A'\n",
    "Number of observations: 105\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the source position. frame='icrs' indicates that we are using coordinates in right ascension and declination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_pos = SkyCoord(83.63311446*u.deg, 22.01448714*u.deg, frame='icrs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selectradius = 2.5*u.deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conesearch = data_store.obs_table.select_sky_circle(source_pos, selectradius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "runlist = conesearch['OBS_ID'].value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(runlist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(runlist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observations = data_store.get_observations(runlist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a dictionary where we will store final results which we produce along the way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's add some information to our results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results['run list'] = runlist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sky Maps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last notebook we binned all events into a two-dimensional map with coordinates Right Ascension and Declination:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_crab = Map.create(binsz=0.01*u.deg, \n",
    "                      width=(5*u.deg, 5*u.deg), \n",
    "                      skydir=source_pos, \n",
    "                      frame='icrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The axes of our map are stored in a geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_crab.geom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now want to bin the events into RA, Dec and energy. We will need a three-dimensional geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# minimum and maximum energy for the analysis\n",
    "Emin =0.1*u.TeV\n",
    "Emax = 50*u.TeV\n",
    "\n",
    "# number of energy bins in the map\n",
    "map_nEbins = 6\n",
    "\n",
    "# width of the map\n",
    "map_width = 8*u.deg\n",
    "\n",
    "# bin size of the map\n",
    "map_binsz = 0.025*u.deg\n",
    "\n",
    "# maximumum offset\n",
    "offset_max = 2.5*u.deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.maps import MapAxis, WcsGeom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_energy_axis = MapAxis.from_energy_bounds(Emin, Emax,\n",
    "                                             nbin=map_nEbins, \n",
    "                                             name=\"energy\")\n",
    "\n",
    "map_geom = WcsGeom.create(skydir=source_pos,\n",
    "                          axes=[map_energy_axis],\n",
    "                          width=map_width,\n",
    "                          binsz=map_binsz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_geom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We will also need an axis for true energy:\n",
    "\n",
    "energy_axis_true = MapAxis.from_energy_bounds(0.08*u.TeV, 80*u.TeV,\n",
    "                                              nbin=8,\n",
    "                                              per_decade=True,\n",
    "                                              name=\"energy_true\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MapDatasetMaker\n",
    "We will now use gammapy code to do the binning of our event lists. Objects that manipulate data in gammapy are called Makers. We start with a DatasetMaker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.datasets import MapDataset\n",
    "\n",
    "from gammapy.makers import MapDatasetMaker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_empty = MapDataset.create(map_geom, \n",
    "                              name='empty', \n",
    "                              energy_axis_true=energy_axis_true\n",
    "                             )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_maker = MapDatasetMaker()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will test everything on the first run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = observations[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_dataset = map_maker.run(map_empty, obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_dataset.counts.plot_grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#map_dataset.counts_off"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SafeMaskMaker\n",
    "We want to reject badly reconstructed events. This could be events with a very large distance from the camera centre (outside the camera) or at low energies where the sensitivity of the instruments is too low. We will reject all events outside a radius of 2.5 deg and where the sensitivity drops below 10% of the maximum value. The SafeMaskMaker will do that for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_dataset.mask_image.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.makers import SafeMaskMaker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "safe_mask_maker = SafeMaskMaker(methods=[\"offset-max\", \"aeff-max\"], \n",
    "                                offset_max=offset_max,\n",
    "                                aeff_percent=10\n",
    "                               )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_dataset = safe_mask_maker.run(map_dataset, obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_dataset.counts.plot_grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_dataset.mask_image.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ring Background\n",
    "Finally we use the RingBackgroundMaker to create a background estimate. First we need an exclusion mask which flags all known and potential sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from regions import CircleSkyRegion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exclusion_regions = [CircleSkyRegion(center=source_pos, radius=0.4*u.deg),  ## exclude the Crab Nebula\n",
    "                     CircleSkyRegion(center=SkyCoord(183.604, -8.708,       ## RGB J0521+212, recommended in gammapy tutorial\n",
    "                                                     unit=\"deg\", \n",
    "                                                     frame=\"galactic\"),\n",
    "                                     radius=0.4*u.deg)\n",
    "                    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "exclusion_mask = map_geom.to_image().region_mask(exclusion_regions, inside = False) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exclusion_mask.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we create the RingBackgroundMaker. We can choose the inner radius of the ring, which must be larger than the largest exclusion region, and the width of the ring."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.makers import RingBackgroundMaker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ring_bkg_maker = RingBackgroundMaker(exclusion_mask=exclusion_mask,\n",
    "                                     r_in=0.4*u.deg,\n",
    "                                     width=0.2*u.deg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_dataset = ring_bkg_maker.run(map_dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_dataset.counts.plot_grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_dataset.counts_off.plot_grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These objects can get quite large. So we need to do a bit of memory management. We delete this object as we don't need it anymore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del map_dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combination of all runs\n",
    "We now loop over all runs and stack the individual data sets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this loop we will delete the individual Datasets and keep only the stacked map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.datasets import Datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_stacked = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_datasets = Datasets()\n",
    "\n",
    "for obs in observations:\n",
    "    map_dataset = map_maker.run(map_empty.copy(name = str(obs.obs_id)), obs)\n",
    "    map_dataset = safe_mask_maker.run(map_dataset, obs)\n",
    "    map_dataset = ring_bkg_maker.run(map_dataset)\n",
    "    \n",
    "    #map_datasets.append(map_dataset)\n",
    "    \n",
    "    if map_stacked is None :\n",
    "        map_stacked = map_dataset.copy(name = 'map stacked')\n",
    "    else :\n",
    "        map_stacked.stack(map_dataset)\n",
    "    \n",
    "    del map_dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#map_stacked = map_datasets.stack_reduce()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_stacked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_stacked.counts.plot_grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_on = map_stacked.counts.sum_over_axes()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_on.plot(add_cbar = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_on.smooth(0.05*u.deg).plot(add_cbar = True,\n",
    "                               cmap = 'coolwarm')\n",
    "\n",
    "#plt.savefig('MapOn.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_stacked.counts_off.plot_grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_stacked.alpha.plot_grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_off = (map_stacked.counts_off * map_stacked.alpha).sum_over_axes()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_off.plot(add_cbar = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_off.smooth(0.05*u.deg).plot(add_cbar = True,\n",
    "                                cmap = 'coolwarm')\n",
    "\n",
    "#plt.savefig('MapOff.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_excess = map_on - map_off"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_excess.plot(add_cbar = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_excess.smooth(0.05*u.deg).plot(add_cbar = True,\n",
    "                                   cmap = 'coolwarm')\n",
    "\n",
    "#plt.savefig('MapExcess.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excess and Significance Maps\n",
    "We could calculate the significance in each bin. There is an easy way to do that. We can use gammapy's ExcessMapEstimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.estimators import ExcessMapEstimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ExcessMapEstimator uses oversampling. So we need to decide our oversampling radius. For that we should have a look at our point-spread function (PSF)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_stacked.psf.plot_containment_radius_vs_energy(fraction=(0.34, 0.68, 0.95, 0.99))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 68% radius is often a good choice. But the Crab Nebula is so bright that we can go lower than that. Let's take 0.05 deg."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimator = ExcessMapEstimator(0.05*u.deg)\n",
    "\n",
    "fluxmaps = estimator.run(map_stacked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxmaps.npred_excess.plot(add_cbar = True, cmap = 'coolwarm')\n",
    "\n",
    "#plt.savefig('FinalMap_Excess.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxmaps.sqrt_ts.plot(add_cbar = True, cmap = 'coolwarm')\n",
    "\n",
    "#plt.savefig('FinalMap_Significance.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxmaps.flux.plot(add_cbar = True, cmap = 'coolwarm')\n",
    "\n",
    "#plt.savefig('FinalMap_Flux.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's keep the results for later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results['excess map'] = fluxmaps.npred_excess.copy()\n",
    "final_results['significance map'] = fluxmaps.sqrt_ts.copy()\n",
    "final_results['flux map'] = fluxmaps.flux.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### your playground\n",
    "Please try different convolution radii. You will see that smaller radii will lead to a more noisy image and larger radii will make the source appear bigger. This does not mean that the source is indeed bigger. You just smear out the emission over a larger area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluation of the Sky Map\n",
    "We need to check that our background estimate is reasonable and we would like to know if there are any other potential sources in the map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will make a histogram of the significance distribution. This will require some data manipulation.v We will make use of numpy, which is typically imported as np."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmapdata = final_results['significance map'].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmapdata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's remove anything that is not a number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.isfinite(sigmapdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmapdata[np.isfinite(sigmapdata)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_hist = plt.hist(sigmapdata[np.isfinite(sigmapdata)],\n",
    "                    bins = 50,\n",
    "                    density = True,\n",
    "                    color = 'tab:blue'\n",
    "                   )\n",
    "\n",
    "plt.yscale('log')\n",
    "\n",
    "#plt.savefig('SigDist_all.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We clearly see that there are bins with significances up to 30 sigma (this is the Crab Nebula). We also many bins with significances around 0 (this is in the parts of the image with no source, which is our baclground). We can also study the regions containing only background, we need to multiply the significance image with the exclusion mask. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_off = final_results['significance map'].reduce_over_axes()*exclusion_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_off.plot(add_cbar = True,\n",
    "             cmap = 'coolwarm'\n",
    "            )\n",
    "\n",
    "#plt.savefig('SigMap_masked.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmapdata_off = sig_off.data[np.isfinite(sig_off.data)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_hist_off = plt.hist(sigmapdata_off,\n",
    "                        bins = sig_hist[1],\n",
    "                        density = True\n",
    "                       )\n",
    "\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should be a Gaussian centred on 0 with a width of 1. Let's test it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu, std = norm.fit(sigmapdata_off)\n",
    "\n",
    "print(f'Fit results: mu = {mu:.2f}, std = {std:.2f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 5, 50)\n",
    "p = norm.pdf(x, mu, std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.stairs(sig_hist_off[0], sig_hist_off[1],\n",
    "           fill = True,\n",
    "           color = 'tab:red'\n",
    "          )\n",
    "\n",
    "plt.plot(x, p, lw=2, color=\"black\")\n",
    "\n",
    "plt.yscale('log')\n",
    "\n",
    "#plt.savefig('SigDist_off.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.stairs(sig_hist[0], sig_hist[1],\n",
    "           fill = True,\n",
    "           color = 'tab:blue',\n",
    "           label='all bins'\n",
    "          )\n",
    "\n",
    "plt.stairs(sig_hist_off[0], sig_hist_off[1],\n",
    "           fill = True,\n",
    "           color = 'tab:red',\n",
    "           label = 'off bins',\n",
    "           alpha = 0.5\n",
    "          )\n",
    "\n",
    "plt.plot(x, p, lw=2, color=\"black\", \n",
    "         label = f'mu = {mu:.2f}, std = {std:.2f}')\n",
    "\n",
    "plt.yscale('log')\n",
    "\n",
    "plt.xlabel('significance')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "#plt.savefig('SigDist_final.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy Spectrum\n",
    "We now want to study the energy distribution of the emission from our object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that we plotted the energy distribution of all the events in one observation run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.events.plot_energy()\n",
    "\n",
    "#plt.savefig('Counts_vs_Energy.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These event counts were recorded within a certain observation time. We need to divide the counts by the time to get a rate (counts per second)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.obs_info['LIVETIME']  ## observation time in seconds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we want to measure a flux. So we still need to divide by the effective area of the instrument. The effective area depends on the energy, but also on the conditions of the observations (like the distance of the source from the camera centre)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.aeff.plot_energy_dependence()\n",
    "\n",
    "plt.xlim(0.4)\n",
    "\n",
    "#plt.savefig('Aeff.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.edisp.peek()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.edisp.to_edisp_kernel(0.7*u.deg).plot_matrix()\n",
    "\n",
    "#plt.savefig('EDisp.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "obs.edisp.plot_migration(offset = 0*u.deg, energy_true = [1,10]*u.TeV)\n",
    "\n",
    "plt.legend(loc='upper right')\n",
    "\n",
    "#plt.savefig('Eresolution.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SpectrumDatasetMaker\n",
    "We will now use gammapy code to do the binning of our event lists. This is very similar to the generation of the MapDataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.datasets import SpectrumDataset\n",
    "\n",
    "from gammapy.makers import SpectrumDatasetMaker"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to use only the events from a small region around the source, the on region. We start with a region at the source position and with a radius of 0.1 degrees. We can change that later if needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "on_region_radius = Angle(0.1*u.deg)\n",
    "\n",
    "on_region = CircleSkyRegion(center=source_pos, radius=on_region_radius)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check on the sky map if the region really encompasses all of the emission. If not, we we need to increase the size. We will use a zoomed version of the significance map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skymap = final_results['significance map'].cutout(source_pos, 1*u.deg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skymap.plot(add_cbar = True, cmap = 'coolwarm')\n",
    "\n",
    "on_region.to_pixel(skymap.geom.wcs).plot(color = 'white')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Emin, Emax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_nEbins = 20\n",
    "\n",
    "spectrum_energy_axis = MapAxis.from_energy_bounds(Emin, Emax,\n",
    "                                                  nbin=spectrum_nEbins,\n",
    "                                                  name=\"energy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.maps import RegionGeom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_geom = RegionGeom.create(region=on_region, \n",
    "                                  axes=[spectrum_energy_axis]\n",
    "                                 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_empty =  SpectrumDataset.create(geom=spectrum_geom, \n",
    "                                         energy_axis_true=energy_axis_true\n",
    "                                        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.makers import SpectrumDatasetMaker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_maker = SpectrumDatasetMaker()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset = spectrum_maker.run(spectrum_empty, obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.counts.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#spectrum_dataset.counts_off.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.mask.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SafeMaskMaker\n",
    "As before, we need to select the good energy and offset range. We will use the same SafeMaskMaker as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset = safe_mask_maker.run(spectrum_dataset, obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "spectrum_dataset.mask.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.energy_range_total"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this run we will use only events with energies of more than 645 GeV."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reflected Background\n",
    "We will use the reflected background to estimate the off source counts. We will need an exclusion mask, we will use the same as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.makers import ReflectedRegionsBackgroundMaker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reflected_bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset = reflected_bkg_maker.run(spectrum_dataset, obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.counts_off.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check on the sky map the location of our off-source regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.visualization import plot_spectrum_datasets_off_regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "skymap = final_results['significance map']\n",
    "\n",
    "ax = skymap.plot()\n",
    "\n",
    "on_region.to_pixel(skymap.geom.wcs).plot(color = 'white')\n",
    "\n",
    "plot_spectrum_datasets_off_regions(ax=ax, datasets=[spectrum_dataset])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combination of all runs\n",
    "We now loop over all runs and stack the individual data sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "spectrum_datasets = Datasets()\n",
    "\n",
    "for obs in observations:\n",
    "    spectrum_dataset = spectrum_maker.run(spectrum_empty.copy(name = str(obs.obs_id)), obs)\n",
    "    spectrum_dataset = safe_mask_maker.run(spectrum_dataset, obs)\n",
    "    spectrum_dataset = reflected_bkg_maker.run(spectrum_dataset, obs)\n",
    "    \n",
    "    spectrum_datasets.append(spectrum_dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_stacked = spectrum_datasets.stack_reduce()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ax = skymap.plot(cmap = 'coolwarm', add_cbar = True)\n",
    "\n",
    "on_region.to_pixel(skymap.geom.wcs).plot(color = 'white')\n",
    "\n",
    "plot_spectrum_datasets_off_regions(ax=ax, datasets=spectrum_datasets)\n",
    "\n",
    "#plt.savefig('ReflectedRegions.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excess and Significance\n",
    "We can now use the SpectrumDatasets to study the excess and the significance of the source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info_table = spectrum_datasets.info_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info_table.colnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info_table['name', 'counts', 'counts_off', 'alpha', 'excess', 'sqrt_ts']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#info_table['name', 'counts', 'counts_off', 'alpha', 'excess', 'sqrt_ts'].write('info_table.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_table = spectrum_datasets.info_table(cumulative=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_table['name', 'counts', 'counts_off', 'alpha', 'excess', 'sqrt_ts']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sum_table['name', 'counts', 'counts_off', 'alpha', 'excess', 'sqrt_ts'].write('sum_table.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last line gives us the on and off counts, as well as excess and significance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "last_line = sum_table[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "last_line['name', 'counts', 'counts_off', 'alpha', 'excess', 'sqrt_ts']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "last_line['counts'] - last_line['alpha'] * last_line['counts_off']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stat = WStatCountsStatistic(n_on=last_line['counts'], \n",
    "                            n_off=last_line['counts_off'], \n",
    "                            alpha=last_line['alpha']\n",
    "                           )\n",
    "\n",
    "print('excess: {} \\nsignificance: {}'.format(stat.n_sig,stat.sqrt_ts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's keep these results for later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results['excess'] = sum_table[-1]['excess']\n",
    "final_results['significance'] = sum_table[-1]['sqrt_ts']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The excess should increase roughly linearly with time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(sum_table['livetime'].to(u.h),\n",
    "             sum_table['excess'],\n",
    "             np.sqrt(sum_table['excess']),\n",
    "             marker='o',\n",
    "             ls='none'\n",
    "            )\n",
    "\n",
    "plt.xlabel('Livetime [h]')\n",
    "plt.ylabel('Excess')\n",
    "\n",
    "#plt.savefig('ExcessVsT.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(sum_table['livetime'].to(u.h),\n",
    "             sum_table['sqrt_ts'],\n",
    "             1,\n",
    "             marker='o',\n",
    "             ls='none'\n",
    "            )\n",
    "\n",
    "plt.xlabel('Livetime [h]')\n",
    "plt.ylabel('Significance')\n",
    "\n",
    "#plt.savefig('SigVsT.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excess and Flux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can subtract the off-source counts from the on-source counts to get the excess:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_stacked.plot_counts()\n",
    "\n",
    "#plt.savefig('Spectrum_OnOffCounts.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we can use (on - alpha*off) to calculate the excess:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "excess = spectrum_stacked.counts.data - spectrum_stacked.counts_off.data * spectrum_stacked.alpha.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "excess = excess.reshape(excess.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_stacked.geoms['geom'].axes['energy'].center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(spectrum_stacked.geoms['geom'].axes['energy'].center,\n",
    "         excess,\n",
    "         marker = 'o',\n",
    "         ls = 'none'\n",
    "        )\n",
    "\n",
    "plt.xscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_stacked.plot_excess()\n",
    "\n",
    "#plt.savefig('Spectrum_ExcessCounts.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could divide the excess by the exposure to get a flux:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_stacked.exposure.plot()\n",
    "\n",
    "#plt.savefig('Spectrum_Exposure.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But the binning is different and the exposure is in true energy. So we need to apply the energy dispersion matrix and to interpolate. We don't need to do that by hand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if False :\n",
    "\n",
    "    flux = spectrum_dataset.excess.quantity / spectrum_dataset.exposure.quantity\n",
    "\n",
    "    plt.plot(spectrum_stacked.geoms['geom'].axes['energy'].center,\n",
    "             flux.reshape(20),\n",
    "             marker = 'o',\n",
    "             ls = 'none'\n",
    "            )\n",
    "\n",
    "    plt.xlabel('E [{}]'.format(spectrum_stacked.geoms['geom'].axes['energy'].center.unit))\n",
    "    plt.ylabel('flux [{}]'.format(flux.unit))\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "\n",
    "    #plt.savefig('SimpleFlux.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectrum Fit\n",
    "Now we want to describe the energy distribution of the gamma rays with a function. We will make a spectral fit. We will stack all the data of all runs into one data set and proceed with the fit. We could also fit the model to each run individually."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start with a simple power law. Remember, the power law is\n",
    "$$\n",
    "f(E) = A \\times \\left( \\frac{E}{E_0} \\right) ^{-\\Gamma}.\n",
    "$$\n",
    "The amplitude $A$ and the spectral index $\\Gamma$ are free parameters in the fit. $E_0$ is the reference energy, which is not fitted. You can freely chose the value of $E_0$, but it is best to keep it within the energy range of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.modeling.models import PowerLawSpectralModel, SkyModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_model = PowerLawSpectralModel(index=2, \n",
    "                                       amplitude=1e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), \n",
    "                                       reference=1 * u.TeV\n",
    "                                      )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the spectral model to a power law with some meaningful start parameters. The spectral model is only a part of a more general model, the SkyModel, which can also contain a spatial and a temporal model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = SkyModel(spectral_model=spectral_model, name=\"Crab\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we set this model as the model of our Dataset. We can also have a look how well our start parameters describe the data already:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.models = model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.plot_excess()\n",
    "\n",
    "plt.ylim(1e-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do the fit we use a gammapy object called Fit. We can reuse the same object later, we need to create it only once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.modeling import Fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit = Fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we run the fit of our Dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit.run(spectrum_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is very important to check the output. Success must be True and everything has to terminate succesfully. Otherwise our fit result is wrong. If the fit fails, try to run it again. Or change the start parameters and run again. Do not continue with a failed fit!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now take a look at our result and the fit parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax_spectrum, ax_residuals = spectrum_dataset.plot_fit()\n",
    "\n",
    "ax_spectrum.set_ylim(0.1, 40)\n",
    "\n",
    "#plt.savefig('PLfit.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.models.to_parameters_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will make a copy of the model for later use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bestmodel_PL = spectrum_dataset.models.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.modeling.models import create_crab_spectral_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "plot_kwargs = {\n",
    "    \"energy_bounds\": [0.1, 30] * u.TeV,\n",
    "    \"ax\": ax,\n",
    "}\n",
    "\n",
    "bestmodel_PL[0].spectral_model.plot(**plot_kwargs, label=\"this work\")\n",
    "bestmodel_PL[0].spectral_model.plot_error(facecolor=\"blue\", alpha=0.3, **plot_kwargs)\n",
    "\n",
    "create_crab_spectral_model(\"hess_pl\").plot(\n",
    "    **plot_kwargs,\n",
    "    label=\"Crab reference\",\n",
    ")\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "#plt.savefig('PLmodel.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also compare this model with a different model. In order to decide which model to be used we will save the test statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TS_PL = spectrum_dataset.stat_sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TS_PL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do the spectrum again, this time we want to fit a power law with an exponential cut-off at high energies. This function is defined as\n",
    "$$\n",
    "f(E) = A \\times \\left( \\frac{E}{E_0} \\right)^{-\\Gamma} \\times \\exp \\left(-\\frac{E}{E_c} \\right).\n",
    "$$\n",
    "The last term can be written as\n",
    "$$\n",
    "\\exp \\left(-\\frac{E}{E_c} \\right) = \\exp \\left(-\\lambda E \\right)\n",
    "$$\n",
    "with $\\lambda = 1/E_c$.\n",
    "\n",
    "Here it is a good idea to use our best fit parameters from the power law as starting parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bestmodel_PL.parameters['index'].quantity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bestmodel_PL.parameters['amplitude'].quantity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_model = ExpCutoffPowerLawSpectralModel(index=bestmodel_PL.parameters['index'].quantity, \n",
    "                                                amplitude=bestmodel_PL.parameters['amplitude'].quantity, \n",
    "                                                reference=1 * u.TeV, \n",
    "                                                lambda_ = 1./(10*u.TeV)\n",
    "                                               )\n",
    "\n",
    "model = SkyModel(spectral_model=spectral_model, name=\"Crab\")\n",
    "\n",
    "spectrum_dataset.models = model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "spectrum_dataset.plot_excess()\n",
    "\n",
    "plt.ylim(1e-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit.run(spectrum_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Don't forget to check for success!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax_spectrum, ax_residuals = spectrum_dataset.plot_fit()\n",
    "\n",
    "ax_spectrum.set_ylim(0.1, 40)\n",
    "\n",
    "#plt.savefig('ExpPLfit.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.models.to_parameters_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bestmodel_expPL = spectrum_dataset.models.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the best-fit cut-off energy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1/(bestmodel_expPL.parameters['lambda_'].quantity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "plot_kwargs = {\n",
    "    \"energy_bounds\": [0.1, 30] * u.TeV,\n",
    "    \"ax\": ax,\n",
    "}\n",
    "\n",
    "bestmodel_expPL[0].spectral_model.plot(**plot_kwargs, label=\"this work\")\n",
    "bestmodel_expPL[0].spectral_model.plot_error(facecolor=\"blue\", alpha=0.3, **plot_kwargs)\n",
    "\n",
    "create_crab_spectral_model(\"hess_ecpl\").plot(\n",
    "    **plot_kwargs,\n",
    "    label=\"Crab reference\",\n",
    ")\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "#plt.savefig('ExpPLmodel.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TS_expPL = spectrum_dataset.stat_sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TS_expPL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare the Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which model is better, with or without the cut-off? We will need to compare the test statistics of the fits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(TS_PL, TS_expPL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gammapy uses a log-likelihood fit. So we can use the likelihood and Wilk's theorem, provided that we are dealing with nested models. The TS value returned by gammapy is already 2xln(L). So only need to take the difference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TS = TS_PL-TS_expPL\n",
    "print(TS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = scipy.stats.chi2.sf(TS,1)\n",
    "\n",
    "print('probabilty: ',P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('significant?', P < 2.7e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model with the cut-off is not significantly better. So we will use the straight power law for our further analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bestmodel = bestmodel_PL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set this model as the current model of the Dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.models = bestmodel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we store the fit parameters in the final results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results['fit parameters'] = bestmodel.parameters.to_table()\n",
    "\n",
    "final_results['model type'] = bestmodel[0].spectral_model.tag[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Your playground\n",
    "You can try to fit yet another model. A LogParabolaSpectralModel could work as well. Check the documentation (https://docs.gammapy.org/0.18.2/api/gammapy.modeling.models.LogParabolaSpectralModel.html#gammapy.modeling.models.LogParabolaSpectralModel) for the parameters of this model.\n",
    "You can compare the log-parabola model with the power law. But you cannot use Wilk's theorem to compare it to the exponential cut-off power law, as they are not nested models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flux Points\n",
    "In the final step we want to create flux points which can be used for later analysis and astrophysical modelling. These points will depend on our best-fit model. Let's check first that we indeed have our best-fit model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.models[0].spectral_model.tag[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.models.parameters.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need only flux points within our energy range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_dataset.energy_range_total"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use 10 bins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_edges = np.geomspace(0.65, 50, 11) * u.TeV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The FluxPointsEstimator will do the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.estimators import FluxPointsEstimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fpe = FluxPointsEstimator(energy_edges = energy_edges,\n",
    "                          selection_optional = ['ul']\n",
    "                         )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points = fpe.run(spectrum_dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "flux_points.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's keep that for later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results['flux points'] = flux_points.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points.plot()\n",
    "\n",
    "#plt.savefig('Flux_dNdE.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flux points are often plotted in $E^2 \\times dN/dE$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points.plot(sed_type = 'e2dnde')\n",
    "\n",
    "#plt.savefig('Flux_E2dNdE.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can save the points for later use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points.write('Crab_spectrum.ecsv',\n",
    "                  overwrite = True\n",
    "                 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to plot the spectral points with our best-fit model then we best use a FluxPointsDataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.datasets import FluxPointsDataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points_dataset = FluxPointsDataset(data = flux_points,\n",
    "                                        models = bestmodel\n",
    "                                       )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points_dataset.plot_fit()\n",
    "\n",
    "#plt.savefig('FluxPoints_wModel.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your playground\n",
    "You can make a spectrum with more or less data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's all. Let's see what we have and summarise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('We have found {} runs'. format(len(final_results['run list'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('We have detected an excess of {:4.1f} gamma rays with a statistical significance of {:3.1f} sigma.'.format(final_results['excess'], final_results['significance']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The spectrum is best described by a {:s}.'.format(final_results['model type']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The best fit parameters are:\\n', final_results['fit parameters'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The spectral data points are:\\n', final_results['flux points'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results['excess map'].plot(cmap = 'coolwarm', add_cbar = 'true')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results['significance map'].plot(cmap = 'coolwarm', add_cbar = 'true')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}