{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting gammapy spectra with sherpa\n",
    "\n",
    "Once we have exported the spectral files (PHA, ARF, RMF and BKG) in the OGIP format, it becomes possible to fit them later with gammapy or with any existing OGIP compliant tool such as XSpec or sherpa.\n",
    "\n",
    "We show here how to do so with sherpa using the high-level user interface. For a general view on how to use stand-alone sherpa, see https://sherpa.readthedocs.io."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data stack\n",
    "\n",
    "- We first need to import the user interface and load the data with [load_data](http://cxc.harvard.edu/sherpa/ahelp/load_data.html).\n",
    "- One can load files one by one, or more simply load them all at once through a [DataStack](http://cxc.harvard.edu/sherpa/ahelp/datastack.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import glob  # to list files\n",
    "from os.path import expandvars\n",
    "from sherpa.astro.datastack import DataStack\n",
    "import sherpa.astro.datastack as sh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sherpa\n",
    "\n",
    "sherpa.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = DataStack()\n",
    "ANALYSIS_DIR = expandvars(\"$GAMMAPY_DATA/joint-crab/spectra/hess/\")\n",
    "filenames = glob.glob(ANALYSIS_DIR + \"pha_obs*.fits\")\n",
    "for filename in filenames:\n",
    "    sh.load_data(ds, filename)\n",
    "\n",
    "# see what is stored\n",
    "ds.show_stack()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define source model\n",
    "\n",
    "We can now use sherpa models. We need to remember that they were designed for X-ray astronomy and energy is written in keV. \n",
    "\n",
    "Here we start with a simple PL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the source model\n",
    "ds.set_source(\"powlaw1d.p1\")\n",
    "\n",
    "# Change reference energy of the model\n",
    "p1.ref = 1e9  # 1 TeV = 1e9 keV\n",
    "p1.gamma = 2.0\n",
    "p1.ampl = 1e-20  # in cm**-2 s**-1 keV**-1\n",
    "\n",
    "# View parameters\n",
    "print(p1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit and error estimation\n",
    "\n",
    "We need to set the correct statistic: [WSTAT](http://cxc.harvard.edu/sherpa/ahelp/wstat.html). We use functions [set_stat](http://cxc.harvard.edu/sherpa/ahelp/set_stat.html) to define the fit statistic, [notice](http://cxc.harvard.edu/sherpa/ahelp/notice.html) to set the energy range, and [fit](http://cxc.harvard.edu/sherpa/ahelp/fit.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define the statistic\n",
    "sh.set_stat(\"wstat\")\n",
    "\n",
    "### Define the fit range\n",
    "ds.notice(0.6e9, 20e9)\n",
    "\n",
    "### Do the fit\n",
    "ds.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results plot\n",
    "\n",
    "Note that sherpa does not provide flux points. It also only provides plot for each individual spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sh.get_data_plot_prefs()[\"xlog\"] = True\n",
    "sh.get_data_plot_prefs()[\"ylog\"] = True\n",
    "ds.plot_fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Errors and confidence contours\n",
    "\n",
    "We use [conf](http://cxc.harvard.edu/sherpa/ahelp/conf.html) and [reg_proj](http://cxc.harvard.edu/sherpa/ahelp/reg_proj.html) functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute confidence intervals\n",
    "ds.conf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute confidence contours for amplitude and index\n",
    "sh.reg_unc(p1.gamma, p1.ampl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Change the energy range of the fit to be only 2 to 10 TeV\n",
    "- Fit the built-in Sherpa exponential cutoff powerlaw model\n",
    "- Define your own spectral model class (e.g. powerlaw again to practice) and fit it"
   ]
  },
  {
   "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
}