{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Accessing variable and transient objects in the truth catalog\n",
    "\n",
    "**Notebook owner**: Yao-Yuan Mao [@yymao](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@yymao).\n",
    "\n",
    "**Last Run:** 2019-02-23\n",
    "\n",
    "In addition to transients being included in images, information on the magnitudes of the transients are included in truth catalogs (input to the image simulation), which can eventually be compared to e.g., the fluxes determined from a difference image. This tutorial shows how to access the truth tables for these transients, and how to search for transients based on properties like magnitude or RA, Dec position.\n",
    "\n",
    "**Learning Objectives:** Use GCR to load summary table and light curves of variable and transient objects in the truth catalogs\n",
    "\n",
    "**Notes:** \n",
    "- Variable and transient objects are **not** available in truth catalog 1.1\n",
    "- To run this notebook, follow the instructions to setup Jupyter-dev at NERSC: https://confluence.slac.stanford.edu/x/1_ubDQ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages and methods that will be used in this notebook\n",
    "\n",
    "import healpy as hp\n",
    "import numpy as np\n",
    "import GCRCatalogs\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary catalog\n",
    "\n",
    "First, we'll load the **summary catalog** (`sc`). The summary catalog is just a simple table that contains the coordinates, types, and id of the objects. \n",
    "\n",
    "The `filters` and `native_filters` options work just like how they work for the static object catalog (see the `truth_gcr_intro` tutorial notebook). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc = GCRCatalogs.load_catalog('dc2_truth_run1.2_variable_summary')\n",
    "sc.list_all_quantities(include_native=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now plot the locations of supernova and AGNs on the sky"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(np.ones(hp.nside2npix(2)), title='', cmap='Greys', min=0.9, max=1.5, cbar=None) # to get an empty sky map\n",
    "\n",
    "data = sc.get_quantities(['ra', 'dec'], native_filters=['sn == 1'])\n",
    "hp.projscatter(data['ra'], data['dec'], lonlat=True, s=1)\n",
    "\n",
    "data = sc.get_quantities(['ra', 'dec'], native_filters=['agn == 1'])\n",
    "hp.projscatter(data['ra'], data['dec'], lonlat=True, s=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## light curve catalog\n",
    "\n",
    "The way to access light curve catalog (`lc`) is a bit different.\n",
    "\n",
    "The data returned by the reader are the light curve tables, which have 5 columns: `uniqueId`, `obshistid`, `mjd`, `filter`, `mag`; and the `filters` option would work on any of these columns.\n",
    "\n",
    "However, the `native_filters` option works differently. They would be imposed on the columns of the \"summary catalog\" we mentioned above, so that one can for example select SNs from a certain RA ranges. \n",
    "\n",
    "Another difference is that, when `return_iterator` is turned on, the reader will iterator over the objects (i.e., returns one light curve table per objects). \n",
    "\n",
    "It's probably easier to illustrate this by an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc = GCRCatalogs.load_catalog('dc2_truth_run1.2_variable_lightcurve')\n",
    "lc.list_all_quantities(include_native=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the i-band light curves for 3 supernovae. \n",
    "\n",
    "Note that here `filter` can be an integer of 0 to 5, and they correspond to `ugrizy` filters respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_plotted = 0\n",
    "for q in lc.get_quantities(['mjd', 'mag'], filters=['filter == 3'], native_filters=['sn == 1'], return_iterator=True):\n",
    "    if len(q['mjd']) >= 20:\n",
    "        plt.plot(q['mjd'], q['mag'], '.-')\n",
    "        n_plotted += 1\n",
    "    if n_plotted >= 3:\n",
    "        break\n",
    "plt.ylim(29, 24)\n",
    "plt.xlabel('MJD')\n",
    "plt.ylabel('$i$ [mag]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1: getting RA and Dec of all SNe whose brightest magnitude is brighter than a certain value\n",
    "\n",
    "In this example, you will need to first iterate over the objects in the light curve table. \n",
    "You can do that by setting `return_iterator=True`. \n",
    "For each object, you then find the brightest r-band magnitude in the light curve, and see if it satisfies your criteria. \n",
    "If so, then you record the object ID (`uniqueId`), which you then use to get RA, Dec from summery table later. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc_data = []\n",
    "for q in lc.get_quantities(['uniqueId', 'mag'], filters=['filter == 2'], native_filters=['sn == 1'], return_iterator=True):\n",
    "    if len(q['mag']) and q['mag'].min() < 25.5:\n",
    "        lc_data.append((q['uniqueId'][0], q['mag'].min()))\n",
    "lc_data = np.array(lc_data, dtype=np.dtype([('uniqueId', int), ('mag_brightest', float)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sc.get_quantities(['uniqueId', 'ra', 'dec'], filters=[(lambda x: np.in1d(x, lc_data['uniqueId'], True), 'uniqueId')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's check if all the IDs actually match\n",
    "# Food for thought: is this always true? What if the order is not preserved between the two tables?\n",
    "\n",
    "(data['uniqueId'] == lc_data['uniqueId']).all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(data['ra'], data['dec'], c=lc_data['mag_brightest'], s=4);\n",
    "plt.xlabel('RA');\n",
    "plt.ylabel('Dec');\n",
    "plt.colorbar(label='mag');\n",
    "\n",
    "# Food for thought: there seem to be different densities in different areas. Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: a histogram of brighest mag for all SNe in a particular area\n",
    "\n",
    "This is an \"opposite\" task of Example 1. Here, we want to plot a histogram of brighest mag for all SNe in a particular area. \n",
    "\n",
    "When querying the light curve table, you can use `native_filters` to filter on the columns of the summary table. \n",
    "Hence, for this task, you **only** need to use the light curve table (surprise!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "area_cut = ['ra > 58', 'dec > -30', 'dec < -28']\n",
    "sn_cut = ['sn == 1']\n",
    "all_cuts = area_cut + sn_cut\n",
    "\n",
    "mags = []\n",
    "for q in lc.get_quantities(['mag'], filters=['filter == 1'], native_filters=all_cuts, return_iterator=True):\n",
    "    if len(q['mag']):\n",
    "        mags.append(q['mag'].min())\n",
    "\n",
    "mags = np.array(mags)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(mags, 20);\n",
    "plt.xlabel('$g$ [mag]');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "desc-python",
   "language": "python",
   "name": "desc-python"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}