{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![CTA first data challenge logo](images/cta-1dc.png)\n",
    "\n",
    "# CTA data analysis with Gammapy\n",
    "\n",
    "## Introduction\n",
    "\n",
    "**This notebook shows an example how to make a sky image and spectrum for simulated CTA data with Gammapy.**\n",
    "\n",
    "The dataset we will use is three observation runs on the Galactic center. This is a tiny (and thus quick to process and play with and learn) subset of the simulated CTA dataset that was produced for the first data challenge in August 2017.\n",
    "\n",
    "**This notebook can be considered part 2 of the introduction to CTA 1DC analysis. Part one is here: [cta_1dc_introduction.ipynb](cta_1dc_introduction.ipynb)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "As usual, we'll start with some setup ..."
   ]
  },
  {
   "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": [
    "!gammapy info --no-envvar --no-system"
   ]
  },
  {
   "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, Angle\n",
    "from astropy.convolution import Gaussian2DKernel\n",
    "from regions import CircleSkyRegion\n",
    "from gammapy.utils.energy import EnergyBounds\n",
    "from gammapy.data import DataStore\n",
    "from gammapy.spectrum import (\n",
    "    SpectrumExtraction,\n",
    "    SpectrumFit,\n",
    "    SpectrumResult,\n",
    "    models,\n",
    "    SpectrumEnergyGroupMaker,\n",
    "    FluxPointEstimator,\n",
    ")\n",
    "from gammapy.maps import Map, MapAxis, WcsNDMap, WcsGeom\n",
    "from gammapy.cube import MapMaker\n",
    "from gammapy.background import ReflectedRegionsBackgroundEstimator\n",
    "from gammapy.detect import TSMapEstimator, find_peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure the logger, so that the spectral analysis\n",
    "# isn't so chatty about what it's doing.\n",
    "import logging\n",
    "\n",
    "logging.basicConfig()\n",
    "log = logging.getLogger(\"gammapy.spectrum\")\n",
    "log.setLevel(logging.ERROR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select observations\n",
    "\n",
    "Like explained in [cta_1dc_introduction.ipynb](cta_1dc_introduction.ipynb), a Gammapy analysis usually starts by creating a `DataStore` and selecting observations.\n",
    "\n",
    "This is shown in detail in the other notebook, here we just pick three observations near the galactic center."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/cta-1dc/index/gps\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just as a reminder: this is how to select observations\n",
    "# from astropy.coordinates import SkyCoord\n",
    "# table = data_store.obs_table\n",
    "# pos_obs = SkyCoord(table['GLON_PNT'], table['GLAT_PNT'], frame='galactic', unit='deg')\n",
    "# pos_target = SkyCoord(0, 0, frame='galactic', unit='deg')\n",
    "# offset = pos_target.separation(pos_obs).deg\n",
    "# mask = (1 < offset) & (offset < 2)\n",
    "# table = table[mask]\n",
    "# table.show_in_browser(jsviewer=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_id = [110380, 111140, 111159]\n",
    "observations = data_store.get_observations(obs_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_cols = [\"OBS_ID\", \"GLON_PNT\", \"GLAT_PNT\", \"LIVETIME\"]\n",
    "data_store.obs_table.select_obs_id(obs_id)[obs_cols]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make sky images\n",
    "\n",
    "### Define map geometry\n",
    "\n",
    "Select the target position and define an ON region for the spectral analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axis = MapAxis.from_edges(\n",
    "    np.logspace(-1.0, 1.0, 10), unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
    ")\n",
    "geom = WcsGeom.create(\n",
    "    skydir=(0, 0), npix=(500, 400), binsz=0.02, coordsys=\"GAL\", axes=[axis]\n",
    ")\n",
    "geom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute images\n",
    "\n",
    "Exclusion mask currently unused. Remove here or move to later in the tutorial?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_position = SkyCoord(0, 0, unit=\"deg\", frame=\"galactic\")\n",
    "on_radius = 0.2 * u.deg\n",
    "on_region = CircleSkyRegion(center=target_position, radius=on_radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exclusion_mask = geom.to_image().region_mask([on_region], inside=False)\n",
    "exclusion_mask = WcsNDMap(geom.to_image(), exclusion_mask)\n",
    "exclusion_mask.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "maker = MapMaker(geom, offset_max=\"2 deg\")\n",
    "maps = maker.run(observations)\n",
    "print(maps.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The maps are cubes, with an energy axis.\n",
    "# Let's also make some images:\n",
    "images = maker.run_images()\n",
    "\n",
    "excess = images[\"counts\"].copy()\n",
    "excess.data -= images[\"background\"].data\n",
    "images[\"excess\"] = excess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Show images\n",
    "\n",
    "Let's have a quick look at the images we computed ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images[\"counts\"].smooth(2).plot(vmax=5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images[\"background\"].plot(vmax=5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images[\"excess\"].smooth(3).plot(vmax=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Source Detection\n",
    "\n",
    "Use the class [gammapy.detect.TSMapEstimator](https://docs.gammapy.org/dev/api/gammapy.detect.TSMapEstimator.html) and [gammapy.detect.find_peaks](https://docs.gammapy.org/dev/api/gammapy.detect.find_peaks.html) to detect sources on the images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = Gaussian2DKernel(1, mode=\"oversample\").array\n",
    "plt.imshow(kernel);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "ts_image_estimator = TSMapEstimator()\n",
    "images_ts = ts_image_estimator.run(images, kernel)\n",
    "print(images_ts.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sources = find_peaks(images_ts[\"sqrt_ts\"], threshold=8)\n",
    "sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_pos = SkyCoord(sources[\"ra\"], sources[\"dec\"])\n",
    "source_pos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot sources on top of significance sky image\n",
    "images_ts[\"sqrt_ts\"].plot(add_cbar=True)\n",
    "\n",
    "plt.gca().scatter(\n",
    "    source_pos.ra.deg,\n",
    "    source_pos.dec.deg,\n",
    "    transform=plt.gca().get_transform(\"icrs\"),\n",
    "    color=\"none\",\n",
    "    edgecolor=\"white\",\n",
    "    marker=\"o\",\n",
    "    s=200,\n",
    "    lw=1.5,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial analysis\n",
    "\n",
    "See other notebooks for how to run a 3D cube or 2D image based analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectrum\n",
    "\n",
    "We'll run a spectral analysis using the classical reflected regions background estimation method,\n",
    "and using the on-off (often called WSTAT) likelihood function.\n",
    "\n",
    "### Extraction\n",
    "\n",
    "The first step is to \"extract\" the spectrum, i.e. 1-dimensional counts and exposure and background vectors, as well as an energy dispersion matrix from the data and IRFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "bkg_estimator = ReflectedRegionsBackgroundEstimator(\n",
    "    observations=observations,\n",
    "    on_region=on_region,\n",
    "    exclusion_mask=exclusion_mask,\n",
    ")\n",
    "bkg_estimator.run()\n",
    "bkg_estimate = bkg_estimator.result\n",
    "bkg_estimator.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "extract = SpectrumExtraction(\n",
    "    observations=observations, bkg_estimate=bkg_estimate\n",
    ")\n",
    "extract.run()\n",
    "observations = extract.spectrum_observations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model fit\n",
    "\n",
    "The next step is to fit a spectral model, using all data (i.e. a \"global\" fit, using all energies)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "model = models.PowerLaw(\n",
    "    index=2, amplitude=1e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n",
    ")\n",
    "fit = SpectrumFit(observations, model)\n",
    "fit.run()\n",
    "print(fit.result[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral points\n",
    "\n",
    "Finally, let's compute spectral points. The method used is to first choose an energy binning, and then to do a 1-dim likelihood fit / profile to compute the flux and flux error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flux points are computed on stacked observation\n",
    "stacked_obs = extract.spectrum_observations.stack()\n",
    "print(stacked_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ebounds = EnergyBounds.equal_log_spacing(1, 40, 4, unit=u.TeV)\n",
    "\n",
    "seg = SpectrumEnergyGroupMaker(obs=stacked_obs)\n",
    "seg.compute_groups_fixed(ebounds=ebounds)\n",
    "\n",
    "fpe = FluxPointEstimator(\n",
    "    obs=stacked_obs, groups=seg.groups, model=fit.result[0].model\n",
    ")\n",
    "flux_points = fpe.run()\n",
    "flux_points.table_formatted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot\n",
    "\n",
    "Let's plot the spectral model and points. You could do it directly, but there is a helper class.\n",
    "Note that a spectral uncertainty band, a \"butterfly\" is drawn, but it is very thin, i.e. barely visible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_result = SpectrumResult(model=fit.result[0].model, points=flux_points)\n",
    "\n",
    "total_result.plot(\n",
    "    energy_range=[1, 40] * u.TeV,\n",
    "    fig_kwargs=dict(figsize=(8, 8)),\n",
    "    point_kwargs=dict(color=\"green\"),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "* Re-run the analysis above, varying some analysis parameters, e.g.\n",
    "    * Select a few other observations\n",
    "    * Change the energy band for the map\n",
    "    * Change the spectral model for the fit\n",
    "    * Change the energy binning for the spectral points\n",
    "* Change the target. Make a sky image and spectrum for your favourite source.\n",
    "    * If you don't know any, the Crab nebula is the \"hello world!\" analysis of gamma-ray astronomy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print('hello world')\n",
    "# SkyCoord.from_name('crab')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What next?\n",
    "\n",
    "* This notebook showed an example of a first CTA analysis with Gammapy, using simulated 1DC data.\n",
    "* This was part 2 for CTA 1DC turorial, the first part was here: [cta_1dc_introduction.ipynb](cta_1dc_introduction.ipynb)\n",
    "* More tutorials (not 1DC or CTA specific) with Gammapy are [here](../index.ipynb)\n",
    "* Let us know if you have any question or issues!"
   ]
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}