{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectrum analysis with Gammapy (run pipeline)\n",
    "\n",
    "In this tutorial we will learn how to perform a 1d spectral analysis.\n",
    "\n",
    "We will use a \"pipeline\" or \"workflow\" class to run a standard analysis. If you're interested in implementation detail of the analysis in order to create a custom analysis class, you should read ([spectrum_analysis.ipynb](spectrum_analysis.ipynb)) that executes the analysis using lower-level classes and methods in Gammapy. \n",
    "\n",
    "In this tutorial we will use the folling Gammapy classes:\n",
    "\n",
    "- [gammapy.data.DataStore](https://docs.gammapy.org/dev/api/gammapy.data.DataStore.html) to load the data to \n",
    "- [gammapy.scripts.SpectrumAnalysisIACT](https://docs.gammapy.org/dev/api/gammapy.scripts.SpectrumAnalysisIACT.html) to run the analysis\n",
    "\n",
    "We use 4 Crab observations from H.E.S.S. for testing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "As usual, we'll start with some setup for the notebook, and import the functionality we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from regions import CircleSkyRegion\n",
    "\n",
    "from gammapy.utils.energy import EnergyBounds\n",
    "from gammapy.data import DataStore\n",
    "from gammapy.scripts import SpectrumAnalysisIACT\n",
    "from gammapy.catalog import SourceCatalogGammaCat\n",
    "from gammapy.maps import Map\n",
    "from gammapy.spectrum.models import LogParabola\n",
    "from gammapy.spectrum import CrabSpectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select data\n",
    "\n",
    "First, we select and load some H.E.S.S. data (simulated events for now). In real life you would do something fancy here, or just use the list of observations someone send you (and hope they have done something fancy before). We'll just use the standard gammapy 4 crab runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\n",
    "mask = data_store.obs_table[\"TARGET_NAME\"] == \"Crab\"\n",
    "obs_ids = data_store.obs_table[\"OBS_ID\"][mask].data\n",
    "observations = data_store.get_observations(obs_ids)\n",
    "print(obs_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Configure the analysis\n",
    "\n",
    "Now we'll define the input for the spectrum analysis. It will be done the python way, i.e. by creating a config dict containing python objects. We plan to add also the convenience to configure the analysis using a plain text config file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crab_pos = SkyCoord.from_name(\"crab\")\n",
    "on_region = CircleSkyRegion(crab_pos, 0.15 * u.deg)\n",
    "\n",
    "model = LogParabola(\n",
    "    alpha=2.3,\n",
    "    beta=0.01,\n",
    "    amplitude=1e-11 * u.Unit(\"cm-2 s-1 TeV-1\"),\n",
    "    reference=1 * u.TeV,\n",
    ")\n",
    "\n",
    "flux_point_binning = EnergyBounds.equal_log_spacing(0.7, 30, 5, u.TeV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exclusion_mask = Map.create(skydir=crab_pos, width=(10, 10), binsz=0.02)\n",
    "\n",
    "gammacat = SourceCatalogGammaCat(\"$GAMMAPY_DATA/gamma-cat/gammacat.fits.gz\")\n",
    "\n",
    "regions = []\n",
    "for source in gammacat:\n",
    "    if not exclusion_mask.geom.contains(source.position):\n",
    "        continue\n",
    "    region = CircleSkyRegion(source.position, 0.15 * u.deg)\n",
    "    regions.append(region)\n",
    "\n",
    "exclusion_mask.data = exclusion_mask.geom.region_mask(regions, inside=False)\n",
    "exclusion_mask.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "config = dict(\n",
    "    outdir=\".\",\n",
    "    background=dict(\n",
    "        on_region=on_region,\n",
    "        exclusion_mask=exclusion_mask,\n",
    "        min_distance=0.1 * u.rad,\n",
    "    ),\n",
    "    extraction=dict(containment_correction=False),\n",
    "    fit=dict(\n",
    "        model=model,\n",
    "        stat=\"wstat\",\n",
    "        forward_folded=True,\n",
    "        fit_range=flux_point_binning[[0, -1]],\n",
    "    ),\n",
    "    fp_binning=flux_point_binning,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the analysis\n",
    "\n",
    "TODO: Clean up the log (partly done, get rid of remaining useless warnings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis = SpectrumAnalysisIACT(observations=observations, config=config)\n",
    "analysis.run(optimize_opts={\"print_level\": 1})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results\n",
    "\n",
    "Let's look at the results, and also compare with a previously published Crab nebula spectrum for reference."
   ]
  },
  {
   "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)\n",
    "CrabSpectrum().model.plot(ax=axes[0], **opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "Rerun the analysis, changing some aspects of the analysis as you like:\n",
    "\n",
    "* only use one or two observations\n",
    "* a different spectral model\n",
    "* different config options for the spectral analysis\n",
    "* different energy binning for the spectral point computation\n",
    "\n",
    "Observe how the measured spectrum changes."
   ]
  },
  {
   "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
}