tutorials/detect_ts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Source detection with Gammapy\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This notebook show how to do source detection with Gammapy using one of the methods available in [gammapy.detect](https://docs.gammapy.org/dev/detect/index.html).\n",
    "\n",
    "We will do this:\n",
    "\n",
    "* produce 2-dimensional test-statistics (TS) images using Fermi-LAT 3FHL high-energy Galactic center dataset\n",
    "* run a peak finder to make a source catalog\n",
    "* do some simple measurements on each source\n",
    "* compare to the 3FHL catalog\n",
    "\n",
    "Note that what we do here is a quick-look analysis, the production of real source catalogs use more elaborate procedures.\n",
    "\n",
    "We will work with the following functions and classes:\n",
    "\n",
    "* [gammapy.maps.WcsNDMap](https://docs.gammapy.org/dev/api/gammapy.maps.WcsNDMap.html)\n",
    "* [gammapy.detect.TSMapEstimator](https://docs.gammapy.org/dev/api/gammapy.detect.TSMapEstimator.html)\n",
    "* [gammapy.detect.find_peaks](https://docs.gammapy.org/dev/api/gammapy.detect.find_peaks.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "As always, let's get started 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": [
    "import numpy as np\n",
    "from astropy import units as u\n",
    "from gammapy.maps import Map\n",
    "from gammapy.detect import TSMapEstimator, find_peaks\n",
    "from gammapy.catalog import source_catalogs\n",
    "from gammapy.cube import PSFKernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in input images\n",
    "\n",
    "We first read in the counts cube and sum over the energy axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts = Map.read(\"$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz\")\n",
    "background = Map.read(\n",
    "    \"$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-background.fits.gz\"\n",
    ")\n",
    "exposure = Map.read(\n",
    "    \"$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-exposure.fits.gz\"\n",
    ")\n",
    "\n",
    "maps = {\"counts\": counts, \"background\": background, \"exposure\": exposure}\n",
    "\n",
    "kernel = PSFKernel.read(\n",
    "    \"$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-psf.fits.gz\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "estimator = TSMapEstimator()\n",
    "images = estimator.run(maps, kernel.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "images[\"sqrt_ts\"].plot(add_cbar=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "images[\"flux\"].plot(add_cbar=True, stretch=\"sqrt\", vmin=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "images[\"niter\"].plot(add_cbar=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Source catalog\n",
    "\n",
    "Let's run a peak finder on the `sqrt_ts` image to get a list of sources (positions and peak `sqrt_ts` values)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sources = find_peaks(images[\"sqrt_ts\"], threshold=8)\n",
    "sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot sources on top of significance sky image\n",
    "plt.figure(figsize=(15, 5))\n",
    "\n",
    "_, ax, _ = images[\"sqrt_ts\"].plot(add_cbar=True)\n",
    "\n",
    "ax.scatter(\n",
    "    sources[\"ra\"],\n",
    "    sources[\"dec\"],\n",
    "    transform=plt.gca().get_transform(\"icrs\"),\n",
    "    color=\"none\",\n",
    "    edgecolor=\"w\",\n",
    "    marker=\"o\",\n",
    "    s=600,\n",
    "    lw=1.5,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measurements\n",
    "\n",
    "* TODO: show cutout for a few sources and some aperture photometry measurements (e.g. energy distribution, significance, flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare to 3FHL\n",
    "\n",
    "TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fermi_3fhl = source_catalogs[\"3fhl\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selection = counts.geom.contains(fermi_3fhl.positions)\n",
    "fermi_3fhl.table = fermi_3fhl.table[selection]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fermi_3fhl.table[[\"Source_Name\", \"GLON\", \"GLAT\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "TODO: put one or more exercises"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start exercises here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What next?\n",
    "\n",
    "In this notebook, we have seen how to work with images and compute TS images from counts data, if a background estimate is already available.\n",
    "\n",
    "Here's some suggestions what to do next:\n",
    "\n",
    "- TODO: point to background estimation examples\n",
    "- TODO: point to other docs ..."
   ]
  }
 ],
 "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
}