DP02_07a_DiaObject_Samples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "# DiaObject Sample Identification\n",
    "\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\">\n",
    "<br>\n",
    "Contact author: Melissa Graham <br>\n",
    "Last verified to run: 2024-12-17 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_50 <br>\n",
    "Container Size: medium <br>\n",
    "Targeted learning level: intermediate <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** To use the DiaObject table parameters to identify a sample of time-variable objects of interest, and to investigate the forced photometry lightcurves of a specific object of interest with the ForcedSourceOnDiaObject table."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Use the TAP service and the DP0.2 DiaObject and DiaSource tables."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** TAP tables dp02_dc2_catalogs.DiaObject, DiaSource, and ForcedSourceOnDiaObject."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c67fab9-136a-4adc-bb42-142b91ab69dd",
   "metadata": {},
   "source": [
    "**Packages:** lsst.rsp, astropy.cosmology, astropy.stats, numpy, matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "metadata": {},
   "source": [
    "**Credit:** Originally developed by Melissa Graham and the Rubin Community Science Team for Data Preview 0. Updated for ForcedSourceOnDiaObject by Ryan Lau."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28e91cbf-ab7f-4e26-9276-b00299d6065e",
   "metadata": {},
   "source": [
    "**Get Support:**\n",
    "Find DP0-related documentation and resources at <a href=\"https://dp0.lsst.io\">dp0.lsst.io</a>. Questions are welcome as new topics in the <a href=\"https://community.lsst.org/c/support/dp0\">Support - Data Preview 0 Category</a> of the Rubin Community Forum. Rubin staff will respond to all questions posted there."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfc73be0",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This notebook demonstrates how to use the DiaObject table's lightcurve summary statistics to identify samples of objects for further study, and how to retrieve and plot lightcurves for DiaObjects of interest from the DiaSource table.\n",
    "\n",
    "As an example, a sample of potential low-redshift Type Ia supernovae are identified and retrieved from the DiaObject table.\n",
    "\n",
    "> **Notice:** This notebook does not produce a pure or complete sample of SNIa, just a sample of *potential well-sampled low-redshift* SNIa. This is decent enough to use as a demonstration, or a starting point for a more rigorous classification process.\n",
    "Different science goals will have different requirements on sample identification.\n",
    "\n",
    "There are three different types of catalogs built from the difference images:\n",
    "1. `DiaSource`: astrometric and photometric measurements for sources detected with signal-to-noise ratio > 5 in the difference images.\n",
    "2. `DiaObject`: derived summary parameters for `DiaSources` associated by sky location, including lightcurve statistics.\n",
    "3. `ForcedSourceOnDiaObject`: point-source forced-photometry measurements on individual single-epoch visit images and difference images, based on and linked to the entries in the `DiaObject` table.\n",
    "\n",
    "Learn more about the contents of the DiaObject, DiaSource, ForcedSourceOnDiaObjects tables in the <a href=\"https://dp0-2.lsst.io\">DP0.2 Documentation</a>.\n",
    "\n",
    "This notebook uses the Table Access Protocol (TAP) service to access the catalogs, but they are also available via the Butler."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc36f107",
   "metadata": {},
   "source": [
    "### 1.1. Package Imports\n",
    "\n",
    "**lsst.rsp:** The LSST Science Pipelines package for RSP functionality such as the TAP service (<a href=\"http://pipelines.lsst.io\">pipelines.lsst.io</a>).\n",
    "\n",
    "**astropy.cosmology:** An open-source package of cosmology tools (<a href=\"https://docs.astropy.org/en/stable/cosmology/index.html\">the astropy cosmology documentation</a>).\n",
    "\n",
    "**astropy.stats:** An open-source package of astronomy-related statistical tools, used in this notebook to look for outliers in the forced-source lightcurves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cddc1458",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "from astropy.stats import sigma_clipped_stats\n",
    "\n",
    "from lsst.rsp import get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c217adff-25ed-4fce-95e7-8aa04630f6cc",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cc807cf-3937-4a92-9fc8-9efbf8827447",
   "metadata": {},
   "source": [
    "Start the TAP service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e32778c-067e-41c3-a777-dcc4128d7de4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "633e872a-137c-4622-8bc7-a32d7f020658",
   "metadata": {},
   "source": [
    "Set the cosmology to use with the astropy.cosmology package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6753ded-fabd-4536-9d3b-66cfee8ebb1f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cosmo = FlatLambdaCDM(H0=70, Om0=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6a8c049-2ea5-4370-b80e-8d2df94464ee",
   "metadata": {},
   "source": [
    "Set a few parameters to use later, when plotting lightcurves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a72c70e-e682-458c-8aeb-56ecfb4567f2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.style.use('tableau-colorblind10')\n",
    "\n",
    "plot_filter_labels = ['u', 'g', 'r', 'i', 'z', 'y']\n",
    "plot_filter_colors = {'u': '#0c71ff', 'g': '#49be61', 'r': '#c61c00',\n",
    "                      'i': '#ffc200', 'z': '#f341a2', 'y': '#5d0000'}\n",
    "plot_filter_symbols = {'u': 'o', 'g': '^', 'r': 'v', 'i': 's', 'z': '*', 'y': 'p'}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b7a23cf-6498-4d36-9a62-a04831371a48",
   "metadata": {},
   "source": [
    "## 2. Understand the DiaObject table's contents\n",
    "\n",
    "**DiaSource Table:** Measurements for sources detected with a signal-to-noise ratio > 5 in the difference images (i.e., lightcurves).\n",
    "Note that these measurements are not forced photometry at the locations of all DiaObjects, but detections only.\n",
    "\n",
    "**DiaObject Table:** Summary parameters for DiaSources associated by coordinate (including lightcurve summary statistics).\n",
    "\n",
    "Option to print all of the available column names for the DiaObject table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab763af8-1dbf-43d2-82c4-a89ae31b1b2f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT column_name from TAP_SCHEMA.columns \"\n",
    "#                          \"WHERE table_name = 'dp02_dc2_catalogs.DiaObject'\")\n",
    "# for column_name in results['column_name']:\n",
    "#     print(column_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1286bda-77bd-44f1-b46a-d991515212ec",
   "metadata": {},
   "source": [
    "### 2.1. Review these descriptions for selected lightcurve summary statistics\n",
    "\n",
    "The DiaObjects table includes the following lightcurve summary statistics, which are derived from the contents of the DiaSource table.\n",
    "\n",
    "**nDiaSources:** The number of difference-image detections in any filter (i.e., number of DiaSources associated with a given DiaObject).\n",
    " \n",
    "The following statistics are all based on _difference-image point source (PS) flux values_ for each filter [f].<br>\n",
    "**[f]PSFluxMin:** The faintest flux. <br>\n",
    "**[f]PSFluxMax:** The brightest flux. <br>\n",
    "**[f]PSFluxMean:** The average flux. <br>\n",
    "**[f]PSFluxSigma:** The standard deviation of the fluxes. <br>\n",
    "**[f]PSFluxMAD:** The mean absolute deviation of the fluxes (i.e., the average distance from the mean). <br>\n",
    "**[f]PSFluxChi2:** The Chi2 statistic for the scatter of the fluxes around the mean. <br>\n",
    "**[f]PSFluxNdata:** The number of data points used to compute [f]PSFluxChi2. <br>\n",
    "**[f]PSFluxSkew:** A measure of asymmentry in the distribution of fluxes about the mean (where 0 means symmetric). <br>\n",
    "**[f]PSFluxStetsonJ:** A variability index developed for Cepheids (defined in <a href=\"https://ui.adsabs.harvard.edu/abs/1996PASP..108..851S/abstract\">Stetson 1996</a>). <br>\n",
    "**[f]PSFluxPercentile05, 25, 50, 75, 95:** Derived from the cumulative distribution of flux values. <br>\n",
    "\n",
    "The following statistics are all based on the _direct-image total (TOT) flux values_ for each filter [f]. <br>\n",
    "**[f]TOTFluxMean:** The average flux.  <br>\n",
    "**[f]TOTFluxSigma:** The standard deviation of the fluxes.  <br>\n",
    "\n",
    "> **Notice:** The DP0.2 DiaObject table is missing some variability characterization parameters (<a href=\"https://dmtn-118.lsst.io\">DMTN-118</a>) and host association parameters (<a href=\"https://dmtn-151.lsst.io\">DMTN-151</a>) which will exist for future data releases."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e56a765a-29c7-49cc-bc30-dcabf9d29514",
   "metadata": {},
   "source": [
    "### 2.2. Plot summary-statistic histograms for a random sample of DiaObjects\n",
    "\n",
    "In order to learn a bit more about these lightcurve summary statistics, histograms of their values for a random selection of DiaObjects are plotted.\n",
    "\n",
    "First, retrieve a random sample of DiaObjects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58c773c7-890f-420d-a6d7-daa7fd9460ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT TOP 100000 \"\\\n",
    "        \"ra, decl, diaObjectId, nDiaSources, \"\\\n",
    "        \"rPSFluxMin, rPSFluxMax, rPSFluxMean, rPSFluxSigma, \"\\\n",
    "        \"rPSFluxMAD, rPSFluxChi2, rPSFluxNdata, rPSFluxSkew, \"\\\n",
    "        \"rPSFluxStetsonJ, rPSFluxPercentile05, rPSFluxPercentile25, \"\\\n",
    "        \"rPSFluxPercentile50, rPSFluxPercentile75, rPSFluxPercentile95, \"\\\n",
    "        \"rTOTFluxMean, rTOTFluxSigma \"\\\n",
    "        \"FROM dp02_dc2_catalogs.DiaObject\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a368e1b5-ba8f-415b-8b6f-c224ff9fbd1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ace0d37-925b-465c-83b0-903d419b63ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "DiaObjs = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73d454f2-ad25-422f-b090-173f0ea7b866",
   "metadata": {},
   "source": [
    "**Plot the number of DiaSources per DiaObject**\n",
    "\n",
    "Plot the distribution of the number of DiaSources per DiaObject (i.e., the total number of difference-image detections in any filter; at left), and the distribution of the number of r-band DiaSources per DiaObject (at right).\n",
    "\n",
    "Notice how the distribution is peaked at small numbers of DiaSources (or r-band detections) -- these are time-variable sources that were only detected in a few difference images, and are most likely to have a faint time-variable flux (e.g., high-_z_ supernovae that were only detectable for a short duration)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5b24c43-caae-49aa-af04-be5bbec4210d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)\n",
    "\n",
    "ax[0].hist(DiaObjs['nDiaSources'], bins=50, log=True, color='gray')\n",
    "ax[0].set_xlabel('nDiaSources')\n",
    "ax[0].set_ylabel('log(Number of DiaObjects)')\n",
    "\n",
    "ax[1].hist(DiaObjs['rPSFluxNdata'], bins=50, log=True, color=plot_filter_colors['r'])\n",
    "ax[1].set_xlabel('rPSFluxNdata')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1253ecb-7149-44ff-b3a7-91cf7212473e",
   "metadata": {},
   "source": [
    "> Figure 1: At left, a histogram (grey) showing the distribution of the number of `diaSource`s per `diaObject`, across all filters. At right, a histogram (dark red) of the number of r-band PSF flux detections per `diaObject`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "360c31db-1ac4-45c7-b0b6-499d5d99292c",
   "metadata": {},
   "source": [
    "**Plot distributions of min/max PS fluxes**\n",
    "\n",
    "Plot the distribution of minimum and maximum r-band PS flux from the difference-image detections.\n",
    "\n",
    "Notice how the PS fluxes can be negative because they are measured on the difference images.\n",
    "Any source that is fainter in the direct image than in the template image will be detected as a negative-flux source in the difference image. \n",
    "\n",
    "> **Warning:** Take care when converting difference-image fluxes into magnitudes, because they can be negative.\n",
    "\n",
    "Notice how these distributions are peaked at low time-variable flux components -- likely the same types of faint sources which were only detectable for a short duration (as shown above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b119d72-6caa-4c7c-b122-d834d58b1314",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)\n",
    "\n",
    "ax[0].hist(DiaObjs['rPSFluxMin'], bins=100, log=True, color=plot_filter_colors['r'])\n",
    "ax[0].set_xlabel('rPSFluxMin')\n",
    "ax[0].set_ylabel('log(Number of DiaObjects)')\n",
    "\n",
    "ax[1].hist(DiaObjs['rPSFluxMax'], bins=100, log=True, color=plot_filter_colors['r'])\n",
    "ax[1].set_xlabel('rPSFluxMax')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66ba1709-02f5-4564-bdc1-763c8f719a73",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-08-19T20:26:21.334487Z",
     "iopub.status.busy": "2024-08-19T20:26:21.333755Z",
     "iopub.status.idle": "2024-08-19T20:26:21.338682Z",
     "shell.execute_reply": "2024-08-19T20:26:21.337941Z",
     "shell.execute_reply.started": "2024-08-19T20:26:21.334459Z"
    }
   },
   "source": [
    "> Figure 2: The distribution of the minimum (left) and maximum (right) r-band PSF flux measurements for `diaObject`s."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c40eb33-9331-4365-9dde-85bfa1e178e2",
   "metadata": {},
   "source": [
    "**Confirm that DiaObjects with few DiaSources have min/max PS fluxes close to zero**\n",
    "\n",
    "Plot the correlation between the number of r-band DiaSources per DiaObject, and the minimum (left) and maximum (right) difference-image PS flux values. \n",
    "\n",
    "Notice that it is, indeed, the fainter DiaObjects that have fewer DiaSource detections, as hypothesized above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "862a4ba2-cff1-42a8-9676-a5625a90d74a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)\n",
    "\n",
    "ax[0].plot(numpy.log10(numpy.absolute(DiaObjs['rPSFluxMin'])),\n",
    "           DiaObjs['rPSFluxNdata'],\n",
    "           'o', ms=5, alpha=0.2, mew=0, color=plot_filter_colors['r'])\n",
    "ax[0].set_xlabel('log |rPSFluxMin|')\n",
    "ax[0].set_ylabel('rPSFluxNdata')\n",
    "\n",
    "ax[1].plot(numpy.log10(numpy.absolute(DiaObjs['rPSFluxMax'])),\n",
    "           DiaObjs['rPSFluxNdata'],\n",
    "           'o', ms=5, alpha=0.2, mew=0, color=plot_filter_colors['r'])\n",
    "ax[1].set_xlabel('log |rPSFluxMax|')\n",
    "ax[1].set_ylabel('rPSFluxNdata')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aa6afb9-7ca8-482b-8a19-d837c70bfb00",
   "metadata": {},
   "source": [
    "> Figure 3: The r-band PSF flux minimum (left) and maximum (right) for `diaObject`s, versus the number of r-band PSF flux detections. This plot demonstrates that `diaObject`s with fewer detections (smaller y-axis values) are fainter (have lower PSF flux values)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cc0a428-ae54-4ff0-a3cf-d81ca60b1f55",
   "metadata": {},
   "source": [
    "**Other PS Flux statistics**\n",
    "\n",
    "It is left as an exercise for the learner to investigate typical values of the other PS Flux summary statistics, such as rPSFluxMean, rPSFluxSigma, rPSFluxMAD, rPSFluxChi2, rPSFluxSkew, and/or rPSFluxStetsonJ."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a63b9b17-e4d8-4a03-91a4-5c59f18f5668",
   "metadata": {},
   "source": [
    "**Plot distributions of TOT flux mean and sigma**\n",
    "\n",
    "Plot the distributions of the DiaObject parameters (mean and sigma) derived from the total fluxes from the direct-images. \n",
    "\n",
    "Note that the TOT fluxes _cannot be_ negative because they are measured on the direct images.\n",
    "The TOT fluxes are safe to convert to magnitudes.\n",
    "\n",
    "Notice how the distribution is not so sharply peaked at low flux, compared to the PS Flux plots above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8970a042-5155-4652-a161-617a2bd3acab",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)\n",
    "\n",
    "ax[0].hist(DiaObjs['rTOTFluxMean'], bins=50, log=True, color=plot_filter_colors['r'])\n",
    "ax[0].set_xlabel('rTOTFluxMean')\n",
    "ax[0].set_ylabel('log(Number of DiaObjects)')\n",
    "\n",
    "ax[1].hist(DiaObjs['rTOTFluxSigma'], bins=50, log=True, color=plot_filter_colors['r'])\n",
    "ax[1].set_xlabel('rTOTFluxSigma')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f5eaf49-0a61-43ec-b333-e4dc7ec16d13",
   "metadata": {},
   "source": [
    "> Figure 4: The distribution of the mean (left) and standard deviation (right) of the total flux (from the 'direct' or science image) measurements for `diaObject`s."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb77e8c5-79a0-485d-b1e1-d07e4971efe2",
   "metadata": {},
   "source": [
    "### 2.3. Investigate the DiaSource data for one random, bright, well-sampled DiaObject\n",
    "\n",
    "Choose a DiaObject that was detected >40 times in an r-band difference image, had an average total (direct-image) flux > 1e6, and an average PS difference-image flux > 5e5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c45cc504-1e97-492d-80e3-69d8ce913a2d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tx = numpy.where((DiaObjs['rPSFluxNdata'] > 40)\n",
    "                 & (DiaObjs['rTOTFluxMean'] > 1000000)\n",
    "                 & (DiaObjs['rPSFluxMean'] > 500000))[0]\n",
    "\n",
    "use_index = tx[0]\n",
    "use_diaObjectId = DiaObjs['diaObjectId'][tx[0]]\n",
    "del tx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b0aec16-2c06-4533-8168-6ca66f2239c0",
   "metadata": {},
   "source": [
    "Retrieve the DiaSource data for this DiaObject."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bfb39b9-6cc1-4b2e-afa5-2873bf790186",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT ra, decl, diaObjectId, diaSourceId, ccdVisitId,\"\\\n",
    "        \"filterName, midPointTai, psFlux, totFlux, totFluxErr, \"\\\n",
    "        \"apFlux, apFluxErr,psFluxErr, snr  \"\\\n",
    "        \"FROM dp02_dc2_catalogs.DiaSource \"\\\n",
    "        \"WHERE diaObjectId = \"+str(use_diaObjectId)\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9378bdce-ec0f-4764-8672-93d447ad1976",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa715e76-1087-47f8-842e-216eb66b6eb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "DiaSrcs = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f00ba02e-5ab8-4cab-80e7-657cdeb4d6e1",
   "metadata": {},
   "source": [
    "Sort the DiaSrcs table based on the mid-point time of the exposure, and create an array `fx` that indexes all observations obtained with the r-band filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b08bfb4-e9b3-4808-9468-49ced614bd7d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "DiaSrcs.sort('midPointTai')\n",
    "fx = numpy.where(DiaSrcs['filterName'] == 'r')[0]\n",
    "print('Number of r-band observations: ', len(fx))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de36bdd9-555b-4a14-90ef-015ddf1a0738",
   "metadata": {},
   "source": [
    "**Plot the lightcurves for one DiaObject**\n",
    "\n",
    "Plot the difference-image (PSFlux) and direct-image (TOTFlux) lightcurves. \n",
    "Mark the DiaObject minimum and maximum flux with solid lines; average flux with dashed lines; and the standard deviation in flux with dotted lines.\n",
    "\n",
    "The purpose of this plot is to serve as a visual aid for understanding the lightcurve statistic summary parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57166c37-81f7-478d-a2db-5aad76820427",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, figsize=(10, 8), sharey=False, sharex=False)\n",
    "\n",
    "ax[0].plot(DiaSrcs['midPointTai'][fx], DiaSrcs['psFlux'][fx],\n",
    "           plot_filter_symbols['r'], ms=15, mew=0, alpha=0.5, color=plot_filter_colors['r'])\n",
    "ax[0].axhline(DiaObjs['rPSFluxMin'][use_index])\n",
    "ax[0].axhline(DiaObjs['rPSFluxMax'][use_index])\n",
    "ax[0].axhline(DiaObjs['rPSFluxMean'][use_index], ls='dashed')\n",
    "ax[0].axhline(DiaObjs['rPSFluxMean'][use_index] - DiaObjs['rPSFluxSigma'][use_index], ls='dotted')\n",
    "ax[0].axhline(DiaObjs['rPSFluxMean'][use_index] + DiaObjs['rPSFluxSigma'][use_index], ls='dotted')\n",
    "ax[0].set_xlabel('Modified Julian Date')\n",
    "ax[0].set_ylabel('PS Flux')\n",
    "ax[0].set_title('Difference-Image PSFlux r-band Lightcurve')\n",
    "\n",
    "ax[1].plot(DiaSrcs['midPointTai'][fx], DiaSrcs['totFlux'][fx],\n",
    "           plot_filter_symbols['r'], ms=15, mew=0, alpha=0.5, color=plot_filter_colors['r'])\n",
    "ax[1].axhline(DiaObjs['rTOTFluxMean'][use_index], ls='dashed')\n",
    "ax[1].axhline(DiaObjs['rTOTFluxMean'][use_index] - DiaObjs['rTOTFluxSigma'][use_index], ls='dotted')\n",
    "ax[1].axhline(DiaObjs['rTOTFluxMean'][use_index] + DiaObjs['rTOTFluxSigma'][use_index], ls='dotted')\n",
    "ax[1].set_xlabel('Modified Julian Date')\n",
    "ax[1].set_ylabel('Total Flux')\n",
    "ax[1].set_title('Direct-Image TOTFlux r-band Lightcurve')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "178246d3-daae-4193-86f2-c805b44ce3c4",
   "metadata": {},
   "source": [
    "> Figure 5: The r-band light curve for a random, well-sampled `diaObject`. The top panel uses the difference-image flux whereas the bottom panel uses the direct-image flux. Horizontal lines mark the minimum and maximum flux (solid), average flux (dashed), and the standard deviation (dotted)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "592a9baa-74a0-41cc-810e-3c8b3ac4df07",
   "metadata": {},
   "source": [
    "**Plot the distribution of flux measurements for one DiaObject**\n",
    "\n",
    "Plot the distribution (left) and normalized cumulative distribution (right) of difference-image fluxes (PSFlux), along with the relevant DiaObject characterization parameters (e.g., mean, sigma, and percentiles).\n",
    "Mark the average flux with dashed lines and the standard deviation in flux with dotted lines.\n",
    "\n",
    "The purpose of this plot is to serve as a visual aid for understanding the lightcurve statistic summary parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "560af696-3468-4c77-aa47-71027f2c131d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)\n",
    "\n",
    "ax[0].hist(DiaSrcs['psFlux'][fx], bins=20, color=plot_filter_colors['r'],\n",
    "           label='skew=%5.2f' % DiaObjs['rPSFluxSkew'][use_index])\n",
    "ax[0].axvline(DiaObjs['rPSFluxMean'][use_index], ls='dashed')\n",
    "ax[0].axvline(DiaObjs['rPSFluxMean'][use_index] - DiaObjs['rPSFluxSigma'][use_index], ls='dotted')\n",
    "ax[0].axvline(DiaObjs['rPSFluxMean'][use_index] + DiaObjs['rPSFluxSigma'][use_index], ls='dotted')\n",
    "ax[0].set_xlabel('r-band PS Flux')\n",
    "ax[0].set_ylabel('Number of DiaSources')\n",
    "ax[0].legend(loc='upper left')\n",
    "\n",
    "ax[1].hist(DiaSrcs['psFlux'][fx], bins=len(fx), color=plot_filter_colors['r'],\n",
    "           cumulative=True, density=True, histtype='step')\n",
    "ax[1].plot(DiaObjs['rPSFluxPercentile05'][use_index], 0.05, '*', ms=10, color='black',\n",
    "           label='percentiles: 0.05, 0.25, 0.50, 0.75, and 0.95')\n",
    "ax[1].plot(DiaObjs['rPSFluxPercentile25'][use_index], 0.25, '*', ms=10, color='black')\n",
    "ax[1].plot(DiaObjs['rPSFluxPercentile50'][use_index], 0.50, '*', ms=10, color='black')\n",
    "ax[1].plot(DiaObjs['rPSFluxPercentile75'][use_index], 0.75, '*', ms=10, color='black')\n",
    "ax[1].plot(DiaObjs['rPSFluxPercentile95'][use_index], 0.95, '*', ms=10, color='black')\n",
    "ax[1].set_xlabel('r-band PS Flux')\n",
    "ax[1].set_ylabel('Cumulative Fraction')\n",
    "ax[1].legend(loc='upper left')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55aae84e-f25f-42c1-a846-f182df75493a",
   "metadata": {},
   "source": [
    "> Figure 6: A visual demonstration of `diaObject` lightcurve parameters. At left, the distribution of r-band PSF difference-image flux measurements for the selected `diaObject`, the skew value of which is noted in the legend. Vertical lines mark the average flux (dashed) and standard deviation (dotted). Note the influence of temporal sampling, as the histogram shows many detections near the maximum PSF flux, yet in Figure 5 it is evident these are temporally clustered observations. At right, the normalized cumulative distribution of flux measurements with percentile values marked as five-point stars."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f28be602-fdbb-40cb-b4b1-56a2259a39ae",
   "metadata": {},
   "source": [
    "Clean up by deleting a few arrays that are no longer needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "656f5166-84ce-415a-a9e4-eadcb90ccc13",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del DiaObjs, DiaSrcs, fx\n",
    "del use_index, use_diaObjectId"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88e53a69-f6dd-4a04-8d00-93d54b325668",
   "metadata": {},
   "source": [
    "## 3. Identify DiaObjects that are potential Type Ia Supernovae\n",
    "\n",
    "For this example, a sample of ***potential*** low-redshift, well-sampled Type Ia supernovae (SNIa) are identified.\n",
    "\n",
    "Compared to other types of supernovae, SNIa have homogenous lightcurves with very similar peak absolute brightnesses (about -19 mag in B-band), and similar rise and decline times (i.e., similar observed lightcurve amplitudes and durations, for a survey with a given limiting magnitude).\n",
    "\n",
    "In LSST-like data sets with an r-band limiting magnitude of about 24.5 mag, such as the DC2 simulation, low-redshift SNIa (0.1 < _z_ < 0.3) lightcurves would look approximately like those plotted in the image below.\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/rubin-dp0/tutorial-notebooks/main/data/template_lowz_Vband.png\" width=400>\n",
    "\n",
    "> Figure 7: Example light curves of SNe Ia at redshift 0.1 (green) and 0.3 (orange), from Nugent et al. (2002; [templates](https://c3.lbl.gov/nugent/nugent_templates.html)).\n",
    "\n",
    "In this example, the expected peak brightnesses, amplitude, and duration for the lightcurves of low-redshift SNIa are used to construct a TAP query on the DiaObjects table to identify a sample of ***potential*** SNIa.\n",
    "\n",
    "> **Warning:** The following is a very rough, back-of-the-envelope way to identify a bunch of DiaObjects with parameters that suggest they are SNIa-like. The resulting sample is not scientifically useful (not pure, not complete) – it is just a few simple cuts that will get you a bunch of DiaObjects for which we can demonstrate how to plot lightcurves using the DiaSource table. Real SN-finding algorithms are photometric classifiers that fit the lightcurves with templates, rigorously, and provide probabilistic classifications. The text in Section 1 of the notebook describes this."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec51ac0b",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.1. Establish lightcurve parameter constraints to identify potential SNIa\n",
    "\n",
    "Define the desired redshift boundaries, the known approximate peak absolute magnitude for SNeIa (-19 mag), and the desired peak range to use to create the TAP query. This range roughly includes the intrinsic diversity in SNIa brightness and color, and host-galaxy reddening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23eb3470-bd49-45a1-a4a2-ab36697c63c6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "redshift_min = 0.1\n",
    "redshift_max = 0.3\n",
    "\n",
    "snia_peak_mag = -19.0\n",
    "snia_peak_mag_range = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d8e4f0a-4b8e-4847-9989-b7fac84f66f5",
   "metadata": {},
   "source": [
    "Use the astropy.cosmology package to convert redshift to distance modulus. \n",
    "Define the minimum and maximum peak apparent r-band magnitudes -- allowing for the intrinsic diversity range specified above -- to use in the \n",
    "TAP query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "862cec58-1dc6-44b7-afed-34a3bbeb12c6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "snia_peak_mr_min = cosmo.distmod(redshift_min).value + snia_peak_mag - snia_peak_mag_range\n",
    "snia_peak_mr_max = cosmo.distmod(redshift_max).value + snia_peak_mag + snia_peak_mag_range\n",
    "print('The minimum and maximum apparent r-band magnitudes '\n",
    "      'to use in the TAP query are %5.2f and %5.2f mag.' %\n",
    "      (snia_peak_mr_min, snia_peak_mr_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab83fd6f-8215-4af7-be4a-d25b07d6fea7",
   "metadata": {},
   "source": [
    "Define maximum magnitudes in the g- and i-bands to use in the TAP query. The point of this is to simply enforce detection in at least the three filters g, r, and i. With knowledge of SNIa colors, this could be made more constraining, but for the purposes of this example these values are fine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1df6cb93-53ee-4f09-9581-455bdd0f5466",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "snia_peak_mg_max = 24.0\n",
    "snia_peak_mi_max = 24.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cadc9f1-0d97-4724-a97a-8a98e0733712",
   "metadata": {},
   "source": [
    "Define the r-band minimum and maximum lightcurve amplitudes to use in the TAP query (i.e., the difference between the brightest and faintest detections in the difference image, in magnitudes). Well-sampled, low-redshift SNIa should be observed to change by at least 1.5 mag (_z_=0.3), and up to 5.5 mag (_z_=0.1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61f81c3b-e7ab-4d66-b534-97a0b477d3ce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "snia_ampl_mr_min = 1.5\n",
    "snia_ampl_mr_max = 5.5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a336a81-cceb-4488-bad7-6b6278c7720a",
   "metadata": {},
   "source": [
    "Define the minimum and maximum number of DiaSources (i.e., difference-image detections) to use in the TAP query.\n",
    "The goal of this example is to identify potential _well-sampled_ Type Ia supernovae, and here a minimum of 15 detections is used.\n",
    "Since the DC2 dataset was simulated using a baseline observing strategy (and does not include deep drilling fields), there are no more than 100 visits per year per field.\n",
    "Any DiaObject with >100 DiaSources had a duration >1 year, and is not a SNIa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d83e88cb-2032-4641-a152-092670595fa6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "snia_nDiaSources_min = 15\n",
    "snia_nDiaSources_max = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f2c49cb-bcb1-4deb-aad9-42d6530c24be",
   "metadata": {},
   "source": [
    "Define the minimum and maximum lightcurve duration, in days. \n",
    "The duration is the time between the first and last difference-image detection in any filter.\n",
    "As seen in the template lightcurve plot above, SNIa at redshifts 0.1 < _z_ < 0.3 will have durations of 50 to 300 days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a77ce538-57cf-4bb8-97d8-a2a5c787b0b9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "snia_duration_min = 50\n",
    "snia_duration_max = 300"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d41d9032-612e-4134-8a2e-4b82470b6268",
   "metadata": {},
   "source": [
    "> **Notice:** Of the above parameters defined to identify potential SNIa, only the lightcurve duration is _not_ represented in the DiaObjects table, and cannot be included as a constraint in the TAP query used below.\n",
    "Instead, the lightcurve durations are calculated and used in Section 3.3."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c2c05d1-7f62-458c-b5e9-cb6a9f399a6d",
   "metadata": {},
   "source": [
    "### 3.2. Retrieve a sample of potentially SNIa-like DiaObjects\n",
    "\n",
    "Only retrieve 1000 DiaObjects for this example.\n",
    "When the TAP query completes, transfer the results to an astropy table.\n",
    "\n",
    "This TAP query takes about a minute.\n",
    "\n",
    "Note the use of the function `scisql_nanojanskyToAbMag` to retrieve magnitudes directly rather than having to convert them from fluxes after the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81b8cd59-1ba3-4eaa-846f-6478ed0c3cf5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT TOP 1000 \"\\\n",
    "        \"ra, decl, diaObjectId, nDiaSources, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(rPSFluxMin) AS rMagMax, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(rPSFluxMax) AS rMagMin, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(gPSFluxMax) AS gMagMin, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(iPSFluxMax) AS iMagMin \"\\\n",
    "        \"FROM dp02_dc2_catalogs.DiaObject \"\\\n",
    "        \"WHERE nDiaSources > \"+str(snia_nDiaSources_min)+\" \"\\\n",
    "        \"AND nDiaSources < \"+str(snia_nDiaSources_max)+\" \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(rPSFluxMax) > \"+str(snia_peak_mr_min)+\" \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(rPSFluxMax) < \"+str(snia_peak_mr_max)+\" \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(gPSFluxMax) < \"+str(snia_peak_mg_max)+\" \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(iPSFluxMax) < \"+str(snia_peak_mi_max)+\" \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(rPSFluxMin)\"\\\n",
    "        \" - scisql_nanojanskyToAbMag(rPSFluxMax) < \"+str(snia_ampl_mr_max)+\" \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(rPSFluxMin)\"\\\n",
    "        \" - scisql_nanojanskyToAbMag(rPSFluxMax) > \"+str(snia_ampl_mr_min)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "749297d6-19c4-4446-9ba8-302865986300",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e920287-d639-4c8f-b749-9b7deeb9f025",
   "metadata": {},
   "outputs": [],
   "source": [
    "DiaObjs = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16f32806-c831-466f-b1fb-1997e3c6b7a3",
   "metadata": {},
   "source": [
    "Calculate the r-band magnitude amplitude as the difference between the brightest and faintests difference-image magnitude, and add it to the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fb8c860-0d1d-408a-9c01-16de5b7a9800",
   "metadata": {},
   "outputs": [],
   "source": [
    "DiaObjs['rMagAmp'] = DiaObjs['rMagMin'] - DiaObjs['rMagMax']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1768203c-de36-4aeb-a447-f3f7e40cdb9c",
   "metadata": {},
   "source": [
    "**Option**: display the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85db3b35-5317-4bb5-888a-56bb56290e2b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# DiaObjs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ff735a9-9a1a-4b5f-8a8a-dfa02f6e2593",
   "metadata": {},
   "source": [
    "**Plot histograms of brightest magnitude and lightcurve amplitude**\n",
    "\n",
    "Plot histograms to characterize the brightest detected r-band magnitudes (left) and the r-band amplitudes (right) for the DiaObjects in this initial sample of potentially SNIa-like time-domain objects.\n",
    "\n",
    "Notice that the brightest r-band magnitude distribution (left) is bimodal, and that the distribution of r-band amplitudes declines until about 3.5 mag and then extends as a tail of brighter time-variable objects to amplitudes of 5.5 mag.\n",
    "The underlying cause of both of these effects is likely the fact that the only time-variable sources in DC2 are stars (mostly brighter) and Type Ia supernovae (mostly fainter) -- but confirming this is left as an exercise for the learner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89c95b1b-159f-461d-9057-c4ea95255ca6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)\n",
    "\n",
    "ax[0].hist(DiaObjs['rMagMin'], bins=20, color=plot_filter_colors['r'])\n",
    "ax[0].set_xlabel('Brightest Detected r-band Magnitude')\n",
    "ax[0].set_ylabel('Number of Potential SNIa')\n",
    "\n",
    "ax[1].hist(DiaObjs['rMagAmp'], bins=20, color=plot_filter_colors['r'])\n",
    "ax[1].set_xlabel('Amplitude in r-band Magnitude')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7059e4ae-664c-467e-b39d-8d4654dc02b6",
   "metadata": {},
   "source": [
    "> Figure 8: At left, the distribution of the brightest difference-image magnitude for all the potential SNIa. At right, the distribution of the difference-image magnitude amplitude for all the potential SNIa."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "424d4c96-66ca-4b0d-a4c8-9aa684981737",
   "metadata": {},
   "source": [
    "### 3.3. Calculate lightcurve duration and identify potential SNIa\n",
    "\n",
    "The lightcurve duration -- time between the first and last detected DiaSource in any filter -- is not included in the DiaObject table.\n",
    "It is calculated below, using all of the DiaSources for each DiaObject.\n",
    "\n",
    "Time is reported in the DiaSource table as `midPointTai`, which is in the SI unit of \"TAI\" (<a href=\"https://en.wikipedia.org/wiki/International_Atomic_Time\">International Atomic Time</a>), and is presented in days (in particular, as \"<a href=\"https://en.wikipedia.org/wiki/Julian_day\">Modified Julian Days</a>\").\n",
    "\n",
    "First, create a string list of the `diaObjectId` that is formatted like \"(1249546586455802795, 1248684569339625641, 1248684569339626125)\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dae27d2-353d-4c6d-b623-b2c0bb0500f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_ids = DiaObjs['diaObjectId']\n",
    "temp_list = '('\n",
    "for i, id in enumerate(all_ids):\n",
    "    temp_list = temp_list + str(id)\n",
    "    if i < len(all_ids)-1:\n",
    "        temp_list = temp_list + ', '\n",
    "    elif i == len(all_ids)-1:\n",
    "        temp_list = temp_list + ')'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3df6731f-5ef8-4c89-9d4f-f4fbbf5578f6",
   "metadata": {},
   "source": [
    "Option to print the 1000-element long list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de12d5ba-cafd-463d-bcbe-20bf2ba673a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(temp_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df1997fd-472b-45e4-a23c-b475b5521e04",
   "metadata": {},
   "source": [
    "Retrieve the `midPointTai` for all observations (`diaSource`s) for all the `diaObject`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5056ddb-88ff-49c6-b1ae-56234d9c0f96",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT diaObjectId, midPointTai \"\\\n",
    "        \"FROM dp02_dc2_catalogs.DiaSource \"\\\n",
    "        \"WHERE diaObjectId IN \" + temp_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4e95b20-40d3-44d3-bdbd-3bc8a3494f5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c9b3a7d-3505-4f23-b845-1f6025f2bbb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03d9f7c7-00b1-44c3-b2a7-aede4247f08e",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49286229-60c9-4c01-9326-4d49d7b3f3d8",
   "metadata": {},
   "source": [
    "Calculate the durations for each `diaObject` and add them to the `DiaObjs` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "910d465f-94ca-4b91-a7ac-9717a4415541",
   "metadata": {},
   "outputs": [],
   "source": [
    "DiaObjs['duration'] = numpy.zeros(len(DiaObjs), dtype='float')\n",
    "\n",
    "for j, DiaObjId in enumerate(DiaObjs['diaObjectId']):\n",
    "    tx = numpy.where(results['diaObjectId'] == DiaObjId)[0]\n",
    "    DiaObjs['duration'][j] = numpy.max(results['midPointTai'][tx]) - numpy.min(results['midPointTai'][tx])\n",
    "    del tx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2353c54e-de16-41d2-af88-be7bc2dcb143",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-01-31T20:39:35.064380Z",
     "iopub.status.busy": "2024-01-31T20:39:35.063650Z",
     "iopub.status.idle": "2024-01-31T20:39:35.068845Z",
     "shell.execute_reply": "2024-01-31T20:39:35.068162Z",
     "shell.execute_reply.started": "2024-01-31T20:39:35.064356Z"
    }
   },
   "source": [
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ef761c3-b610-4393-b5a5-9067afe7c616",
   "metadata": {},
   "outputs": [],
   "source": [
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16aafdfa-e78d-4345-982e-c4ee37b44489",
   "metadata": {},
   "source": [
    "Select only DiaObjects that have lightcurve durations within the specified range for SNIa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aed21a5f-101c-44bb-868c-a3b18e97fe27",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tx = numpy.where((DiaObjs['duration'] > snia_duration_min)\n",
    "                 & (DiaObjs['duration'] < snia_duration_max))[0]\n",
    "print(len(tx))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23cdc39f-3fbb-4bfb-87d3-0e07f5c34096",
   "metadata": {},
   "source": [
    "**Plot a histogram of lightcurve durations, and a scatter plot of duration vs. amplitude, for the potential SNIa**\n",
    "\n",
    "Plot a histogram of lightcurve durations (left) and a plot of duration versus r-band amplitude (right) to further characterize this sample of potential SNIa. \n",
    "\n",
    "Notice the spread in duration, and how it is correlated with lightcurve amplitude, as expected for SNIa. \n",
    "\n",
    "Tracking down the origin of the two outliers is left as an exercise for the learner: if they are truly SNIa observed for >100 days, then they should have a brighter amplitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "508ae771-c983-4299-89a4-97b497490b95",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)\n",
    "\n",
    "ax[0].hist(DiaObjs['duration'][tx], bins=20, color='gray')\n",
    "ax[0].set_xlabel('Lightcurve Duration [days]')\n",
    "ax[0].set_ylabel('Number of Potential SNIa')\n",
    "\n",
    "ax[1].plot(DiaObjs['duration'][tx], DiaObjs['rMagAmp'][tx], 'o', color='gray')\n",
    "ax[1].set_xlabel('Lightcurve Duration [days]')\n",
    "ax[1].set_ylabel('Amplitude in r-band [mag]')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7571a60a-b0fe-4d72-bcd7-ae4280dd21db",
   "metadata": {},
   "source": [
    "> Figure 9: At left, the distribution of the lightcurve duration, in days, for potential SNIa. At right, the lightcurve duration versus r-band amplitude."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f6216ee-2188-43e0-b188-a70f42d2f11e",
   "metadata": {},
   "source": [
    "### 3.4. Plot multi-band lightcurves for the first 20 potential SNIa\n",
    "\n",
    "Notice that all of these lightcurves do resemble Type Ia supernovae!\n",
    "That is a promising sign that the methodology of using the DiaObject summary parameters, along with the lightcurve duration (which must be calculate separately), can identify an initial sample of potential SNIa.\n",
    "\n",
    "> **Keep in mind that this is not necessarily a pure or complete sample of SNIa,** just a sample of *potential well-sampled low-redshift* SNIa - decent enough to use as a demonstration, or a starting point for a more rigorous classification process.\n",
    "Different science goals will have different requirements on sample identification.\n",
    "\n",
    "Below, the lightcurves plotted are based on the signal-to-noise ratio > 5 detections from the DiaSource table only.\n",
    "In Section 4, lightcurves based on forced photometry are used instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41bb4666-7c90-4f35-b9e8-b3dbd59b0185",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(5, 4, figsize=(10, 10), sharey=False, sharex=False)\n",
    "\n",
    "x = 0\n",
    "\n",
    "for i in range(5):\n",
    "    for j in range(4):\n",
    "        results = service.search(\"SELECT ra, decl, diaObjectId, diaSourceId, \"\n",
    "                                 \"filterName, midPointTai, \"\n",
    "                                 \"scisql_nanojanskyToAbMag(psFlux) AS psAbMag \"\n",
    "                                 \"FROM dp02_dc2_catalogs.DiaSource \"\n",
    "                                 \"WHERE diaObjectId = \"+str(DiaObjs['diaObjectId'][tx[x]]))\n",
    "        results = results.to_table()\n",
    "\n",
    "        for f, filt in enumerate(plot_filter_labels):\n",
    "            fx = numpy.where(results['filterName'] == filt)[0]\n",
    "            ax[i, j].plot(results['midPointTai'][fx], results['psAbMag'][fx],\n",
    "                          plot_filter_symbols[filt], ms=10, mew=0, alpha=0.5,\n",
    "                          color=plot_filter_colors[filt])\n",
    "            del fx\n",
    "\n",
    "        ax[i, j].invert_yaxis()\n",
    "        ax[i, j].set_title(DiaObjs['diaObjectId'][tx[x]])\n",
    "\n",
    "        if i == 4:\n",
    "            ax[i, j].xaxis.set_label_text('MJD (days)')\n",
    "        if j == 0:\n",
    "            ax[i, j].yaxis.set_label_text('mag')\n",
    "\n",
    "        x += 1\n",
    "        del results\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8592eaa7-7480-4153-bcc2-d6af68ac6a39",
   "metadata": {},
   "source": [
    "> Figure 10: Multi-band lightcurves for 20 potential SNIa. Legend: $u$, blue circle; $g$, green triangle; $r$, dark red inverted triangle; $i$, golden square; $z$, pink star; $y$, brown pentagram."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "251cbca7-8b24-437f-9442-e140133ef96a",
   "metadata": {},
   "source": [
    "## 4. Forced photometry lightcurves\n",
    "\n",
    "Above, the DiaSource table is used to plot lightcurves, but it only includes SNR>5 detections in the difference images.\n",
    "If your science goal requires lower-SNR measurements, use the forced photometry in the ForcedSourceOnDiaObjects table instead, as demonstrated below.\n",
    "\n",
    "Let's say there is a specific object whose lightcurve looks particularly interesting. In order to explore the flux of the object in all sets of differences images (not just those with detected sources), the ForcedSourceOnDiaObject catalog can be utilized. \n",
    "\n",
    "**ForcedSourceOnDiaObject Table**: point-source forced-photometry measurements on individual single-epoch visit images and difference images, based on and linked to the entries in the `DiaObject` table.\n",
    "\n",
    "For more information, see <a href=\"https://dm.lsst.org/sdm_schemas/browser/dp02.html#ForcedSourceOnDiaObject\">the schema of the ForcedSourceOnDiaObject catalog</a>.\n",
    "\n",
    "**Option**: print all of the available column names and descriptions for the ForcedSourceOnDiaObject table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7a14ffd-9999-47d7-8e5a-3f71c3f329a1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT column_name, description from TAP_SCHEMA.columns \"\n",
    "#                          \"WHERE table_name = 'dp02_dc2_catalogs.ForcedSourceOnDiaObject'\")\n",
    "# for column_name, description in zip(results['column_name'], results['description']):\n",
    "#     print(column_name, '-', description)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d0bae65-b864-4938-9a7d-ee2509db6bda",
   "metadata": {},
   "source": [
    "### 4.1. Retrieve photometry for a DiaObject\n",
    "\n",
    "As an example, let's look at the forced photometry on the difference images of the object 1250953961339360185, which has an interesting lightcurve. Note that it is unlikely this object is shown in the 20 SNIa candidates above since they will change every time the search is run. That object should however be in the list of 1000 SNIa candidates that we started with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "298f86a6-d019-4a93-81c4-7bf0ccdc014b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "DiaObjID = 1250953961339360185\n",
    "print(DiaObjID)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3259766-b214-4564-8739-2efec82fdaa4",
   "metadata": {},
   "source": [
    "Here, we will define the query to extract the information from the ForcedSourceOnDiaObject table on our DiaObject of interest in order to construct a lightcurve that includes forced photometry (i.e., sub-threshold measurements with SNR<5).\n",
    "\n",
    "Below, a JOIN to the CcdVisit table is used to obtain the exposure time mid-point in modified julian date (MJD) format (expMidptMJD), which is necessary to construct the lightcurve.\n",
    "\n",
    "The ForcedSourceOnDiaObject contains forced photometry on both the difference image (psfDiffFlux, psfDiffFluxErr) and the processed visit image (PVI), also called the \"direct\" image (psfFlux, psfFluxErr).\n",
    "Both types of fluxes are retrieved below for our DiaObject of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9abfbc2b-3269-4748-80ef-0e7d851f40f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT fsodo.coord_ra, fsodo.coord_dec, \"\\\n",
    "        \"fsodo.diaObjectId, fsodo.ccdVisitId, fsodo.band, \"\\\n",
    "        \"fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, \"\\\n",
    "        \"fsodo.psfFlux, fsodo.psfFluxErr, \"\\\n",
    "        \"cv.expMidptMJD \"\\\n",
    "        \"FROM dp02_dc2_catalogs.ForcedSourceOnDiaObject as fsodo \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.CcdVisit as cv ON cv.ccdVisitId = fsodo.ccdVisitId \"\\\n",
    "        \"WHERE fsodo.diaObjectId = \"+str(DiaObjID)\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a986514f-9e93-4c44-b73c-865651dd984e",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "282cbca1-d5e3-4499-837c-d3d7d25b2d3c",
   "metadata": {},
   "source": [
    "Save the results in a table \"FrcdSrc\" (short for forced source, to distinguish it from the results from DiaSource saved as DiaSrc, below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1af5d4d-624b-48fe-9efd-88d243561e2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "FrcdSrc = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91724b8a-a56d-4e8c-9a49-0fc9d2c057a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FrcdSrc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2609a6c3-c9fa-49c4-a70a-48a4dc4723af",
   "metadata": {},
   "source": [
    "Also retrieve the DiaSource measurements for the DiaObject of interest.\n",
    "Recall that this will only return flux measurements from the difference images with SNR>5 detection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cc6be77-26f0-455b-b231-0e997fb0aae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT ra, decl, diaObjectId, diaSourceId, psFlux, psFluxErr, \"\\\n",
    "        \"filterName, midPointTai, SNR, ccdVisitId, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(psFlux) AS psAbMag \"\\\n",
    "        \"FROM dp02_dc2_catalogs.DiaSource \"\\\n",
    "        \"WHERE diaObjectId = \"+str(DiaObjID)\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "322391e5-82e4-402f-bbba-8486539bcd89",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63d1cf14-9d69-4d03-9278-41b8658080d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "DiaSrc = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce60a0dd-b89f-4d96-a2c5-9134661e1540",
   "metadata": {},
   "source": [
    "Unlike the DiaSource catalog, the ForcedSourceOnDiaObject catalog does not provide the SNR for the flux measurements. Since it provides a useful diagnostic, we define 'diffSNR' for our object by the ratio of the absolute value of 'psfDiffFlux' and 'psfDiffFluxErr.'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a701ae4-874f-4558-a1b1-f63270ef64e7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "FrcdSrc['DiffSNR'] = abs(FrcdSrc['psfDiffFlux']/FrcdSrc['psfDiffFluxErr'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9a5d68f-1aaa-4067-8218-4722805b6803",
   "metadata": {},
   "source": [
    "Option to show that the DiffSNR column has been added to our \"FrcdSrc\" table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ff11d98-a048-418d-a9f4-179fed460b79",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# FrcdSrc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1b6da40-19e7-4191-a363-6c8beb0f3714",
   "metadata": {},
   "source": [
    "### 4.2. Plot the lightcurves\n",
    "\n",
    "For our DiaObject of interest, compare the lightcurves that include forced photometry from the ForcedSourceOnDiaObject table (open markers, left plots) with lightcurves that include only SNR>5 detections from the DiaSource table (filled markers, right plots).\n",
    "\n",
    "> **Notice**: The forced photometry lightcurves (left) include flux measurements from before and after the source was identified in the DIAsource catalogs. The shaded gray regions indicate the dates spanned by the detections only (DiaSources only).\n",
    "\n",
    "Define the MJD ranges for the \"full\" lightcurve from the forced photometry, and for the \"detected\" part of the lightcurve from the DiaSource table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fe9873e-2da6-46f4-9677-9213847c9585",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mjd_full_start = numpy.min(FrcdSrc['expMidptMJD'])\n",
    "mjd_full_end = numpy.max(FrcdSrc['expMidptMJD'])\n",
    "\n",
    "mjd_det_start = numpy.min(DiaSrc['midPointTai'])\n",
    "mjd_det_end = numpy.max(DiaSrc['midPointTai'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5d82867-4ade-4c86-98fe-fbeea609d833",
   "metadata": {},
   "source": [
    "#### 4.2.1. Lightcurve for the full survey"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dfb35a4-d68c-4fd0-a482-0402142d0e67",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(12, 4), sharey=True, sharex=False)\n",
    "\n",
    "for f, filt in enumerate(plot_filter_labels):\n",
    "    fx = numpy.where(FrcdSrc['band'] == filt)[0]\n",
    "    ax[0].plot(FrcdSrc['expMidptMJD'][fx], FrcdSrc['psfDiffFlux'][fx],\n",
    "               plot_filter_symbols[filt], ms=10, mew=2, mec=plot_filter_colors[filt],\n",
    "               alpha=0.5, color='none', label=filt)\n",
    "    del fx\n",
    "\n",
    "ax[0].set_xlabel('Modified Julian Date')\n",
    "ax[0].set_ylabel('Difference-Image Flux')\n",
    "ax[0].set_title('Forced Photometry (ForcedSourceOnDiaObject)')\n",
    "ax[0].set_xlim(mjd_full_start, mjd_full_end)\n",
    "ax[0].legend(loc='upper left')\n",
    "ax[0].axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')\n",
    "\n",
    "for f, filt in enumerate(plot_filter_labels):\n",
    "    fx = numpy.where(DiaSrc['filterName'] == filt)[0]\n",
    "    ax[1].plot(DiaSrc['midPointTai'][fx], DiaSrc['psFlux'][fx],\n",
    "               plot_filter_symbols[filt], ms=10, mew=0, color=plot_filter_colors[filt],\n",
    "               alpha=0.5, label=filt)\n",
    "    del fx\n",
    "\n",
    "ax[1].set_xlabel('Modified Julian Date')\n",
    "ax[1].set_title('SNR>5 Detections (DiaSource)')\n",
    "ax[1].set_xlim(mjd_full_start, mjd_full_end)\n",
    "ax[1].legend(loc='upper left')\n",
    "ax[1].axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')\n",
    "plt.subplots_adjust(wspace=.0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89410825-2ae5-4ea9-9ec6-15266158ab4f",
   "metadata": {},
   "source": [
    "> Figure 11: The multi-band lightcurve for the chosen `diaObject` based on forced photometry in all difference images (left) and measurement photometry for signal-to-noise ratio $>$ 5 difference-image detections only (right)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b28bf1e6-6362-4a9a-83c3-48107b923b94",
   "metadata": {},
   "source": [
    "#### 4.2.2. Lightcurve for epochs with SNR>5 detections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64eb4993-4684-41f8-a5a5-72197c41f066",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(12, 4), sharey=True, sharex=False)\n",
    "\n",
    "for f, filt in enumerate(plot_filter_labels):\n",
    "    fx = numpy.where(FrcdSrc['band'] == filt)[0]\n",
    "    ax[0].plot(FrcdSrc['expMidptMJD'][fx], FrcdSrc['psfDiffFlux'][fx],\n",
    "               plot_filter_symbols[filt], ms=10, mew=2, mec=plot_filter_colors[filt],\n",
    "               alpha=0.5, color='none', label=filt)\n",
    "    del fx\n",
    "ax[0].set_xlabel('Modified Julian Date')\n",
    "ax[0].set_ylabel('Difference-Image Flux')\n",
    "ax[0].set_title('Forced Photometry (ForcedSourceOnDiaObject)')\n",
    "ax[0].set_xlim(mjd_det_start-50, mjd_det_end+50)\n",
    "ax[0].legend(loc='upper left')\n",
    "ax[0].axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')\n",
    "\n",
    "for f, filt in enumerate(plot_filter_labels):\n",
    "    fx = numpy.where(DiaSrc['filterName'] == filt)[0]\n",
    "    ax[1].plot(DiaSrc['midPointTai'][fx], DiaSrc['psFlux'][fx],\n",
    "               plot_filter_symbols[filt], ms=10, mew=0, color=plot_filter_colors[filt],\n",
    "               alpha=0.5, label=filt)\n",
    "    del fx\n",
    "\n",
    "ax[1].set_xlabel('Modified Julian Date')\n",
    "ax[1].set_title('SNR>5 Detections (DiaSource)')\n",
    "ax[1].set_xlim(mjd_det_start-50, mjd_det_end+50)\n",
    "ax[1].legend(loc='upper left')\n",
    "ax[1].axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')\n",
    "plt.subplots_adjust(wspace=.0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b3203ea-1512-4046-bb41-934a20d91699",
   "metadata": {},
   "source": [
    "> Figure 12: As in Figure 11, but zoomed-in on the dates for which the `diaObject` was detected."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "052bd8aa-a6eb-4e93-9110-48670550078c",
   "metadata": {},
   "source": [
    "> **Notice**: In the above plots, the forced photometry (psfDiffFlux from ForcedSourceOnDiaObject, open symbols, left panels) is nearly identical to the detection photometry (psFlux from DiaSource, filled symbols, right panels), but the ForcedSourceOnDiaObject contains additional measurements from epochs in which the source in the difference image had a SNR<5. No such cutoff is applied to the ForcedSourceOnDiaObject catalog, in which photometry from all visits is included.\n",
    "\n",
    "### 4.3. Identify \"missing measurements\" from SNR<5 epochs\n",
    "\n",
    "Within the \"gray zone\" above (defined by epochs in which the source was detected with SNR>5 in the difference images), identify how many measurements there are in the DiaSource and ForcedSourceOnDiaObject tables, and compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c4681de-b970-4c00-8f50-55a68907d611",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cov = numpy.where((FrcdSrc['expMidptMJD'] >= mjd_det_start)\n",
    "                  & (FrcdSrc['expMidptMJD'] <= mjd_det_end))\n",
    "\n",
    "print('Within the date range of SNR>5 detections, there were:')\n",
    "print(len(DiaSrc), 'measurements in DiaSource catalog of %.i' % (DiaObjID))\n",
    "print(len(FrcdSrc[cov]), 'measurements in ForcedSourceOnDiaObject catalog of %.i' % (DiaObjID))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0fbe561-56c3-4fa2-b233-b93cdbfc1b94",
   "metadata": {},
   "source": [
    "Identify only the SNR<5 measurements from the ForcedSourceOnDiaObject table by looking for unique \"ccdVisitId\" values that are not in the DiaSource results.\n",
    "Note that we are still applying the condition where the measurements are taken within the date range of the DiaSource detections.\n",
    "\n",
    "First, identify the unique visit IDs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88694b90-e64b-4f3b-bd8e-9748ac9b6cff",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "UniqueVID = set(FrcdSrc[cov]['ccdVisitId']) - set(DiaSrc['ccdVisitId'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bee4c7bd-b501-4f82-91b4-0def50993753",
   "metadata": {},
   "source": [
    "Next, construct a table of ForcedSourceOnDiaObject \"missing measurements\" that are not in the DiaSource catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5bcb312-3e8c-449c-bdb6-8422b976857e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "UniqueObj = []\n",
    "for x in UniqueVID:\n",
    "    UniqueObj += numpy.where(FrcdSrc['ccdVisitId'] == x)\n",
    "UniqueObj = numpy.concatenate(UniqueObj)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc3c0e76-dd5f-4888-8f11-21e507a852cd",
   "metadata": {},
   "source": [
    "Inspect the table of \"missing measurements\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80f12a04-b2bc-4f9b-8409-6909381fb8cd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "FrcdSrc[UniqueObj]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0721e50-1c9c-46fd-b1ad-8434502e5adf",
   "metadata": {},
   "source": [
    "Finally, verify why these \"missing measurements\" are not in the DiaSource catalog (their SNR is <5) by printing the maximum SNR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cd157fe-4d8a-42d2-8be7-2f16265bf704",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "numpy.max(FrcdSrc[UniqueObj]['DiffSNR'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73bca2c0-4872-4173-8ee3-5fb44befdfde",
   "metadata": {},
   "source": [
    "### 4.4. Investigate potential SNR<5 variability\n",
    "\n",
    "Going back to the forced photometry lightcurve of our object of interest, it looked like there may have been some past variability in the z band (purple star) at around MJD 60250. Let's investigate this point. **This is for demonstration purposes only since there are no simulated precursor outbursts for the supernovae in DP0.2**. \n",
    "\n",
    "Let's check if it is at least greater than 3$\\sigma$ of the baseline (median) measurements using \"sigma_clipped_stats\"."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ac4a87b-7d8d-4def-b764-9db39353ee73",
   "metadata": {},
   "source": [
    "Define the filter in which flaring (precursor activity) is investigated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5e75997-954e-443a-b846-fa65051474e3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "flare_filt = 'z'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30219e89-92d3-4ee1-a602-c011e8471561",
   "metadata": {},
   "source": [
    "Specify the median and standard deviation in the \"flare\" filter band for the flux measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f9b613e-b48e-4585-862a-dca1c939a154",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "_, med, std = sigma_clipped_stats(FrcdSrc['psfDiffFlux'][numpy.where(FrcdSrc['band'] == flare_filt)[0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d66f2b38-54d1-42b1-a340-1ba0fa8c7013",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print('Median psfDiffFlux for all z-band forced sources: ', med, FrcdSrc['psfDiffFlux'].unit)\n",
    "print('Standard deviation in psfDiffFlux for all z-band forced sources: ',std, FrcdSrc['psfDiffFlux'].unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef2dc868-9912-41d2-b732-05273003eeef",
   "metadata": {},
   "source": [
    "Define where any outlier variability (above 3$\\sigma$) might be in the \"flare\" band, BEFORE or AFTER the \"gray zone\" defined by the epochs in which our DiaObject of interest was detected with SNR>5 in any filter.\n",
    "\n",
    "First, index the flare-band epochs \"before\" and \"after\" within FrcdSrc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfee369c-8866-4e67-afd8-14adda14ab94",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ex = numpy.where(((FrcdSrc['expMidptMJD'] < mjd_det_start)\n",
    "                  | (FrcdSrc['expMidptMJD'] > mjd_det_end))\n",
    "                 & (FrcdSrc['band'] == flare_filt))[0]\n",
    "print(len(ex))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03377608-7760-4170-88a4-f44e83954202",
   "metadata": {},
   "source": [
    "Next, index which of _these_ is >3$\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfab7f63-20a4-41c4-9cde-bbf7505d3645",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ix = numpy.where(FrcdSrc['psfDiffFlux'][ex] > med+3*std)[0]\n",
    "print('Number of >3𝜎 z-band measurements (outside SNR>5 detections):',len(ix))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b41250d3-4058-4fdb-8c1f-214747f99ff3",
   "metadata": {},
   "source": [
    "Inspect the single >3$\\sigma$ measurement in the forced photometry that occurred outside of the \"gray zone\".\n",
    "\n",
    "**As a reminder, 3$\\sigma$ is defined from the distribution of forced photometry measurements and is _not_ the photometric signal-to-noise ratio (i.e., DiffSNR)**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4542760b-ea27-46ea-9138-b9bcfc9a60b3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "FrcdSrc[ex[ix]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd561917-30cd-461a-a585-cb3e8f8be2d6",
   "metadata": {},
   "source": [
    "**Statistics check:** recall that $\\pm 3\\sigma$ encompasses 99.7% of values in a normally-distributed random sample. So out of the 77 z-band forced photometry measurements we could have reasonably expected 0.23 at the $\\sigma$ > 3 level, and the fact that we found one is consistent with small-number statistics, i.e., consistent with noise.\n",
    "\n",
    "**Again, since this is from DP0.2, we _know_ this \"flaring activity\" is just noise but we use it to demonstrate how the ForcedSourceOnDIAObject catalog can be applied to search for pre- or post-transient variability.**\n",
    "\n",
    "As a final demonstration, plot the point(s) on the lightcurve that exceed 3$\\sigma$ and are outside the DiaSource date coverage.\n",
    "Below, we plot the 3$\\sigma$ flux range in red and the apparent variability outlier(s) as a yellow-filled plot marker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "255bf00e-0331-4c15-8126-43efbf23fd48",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(7, 4))\n",
    "\n",
    "for f, filt in enumerate(plot_filter_labels):\n",
    "    fx = numpy.where(FrcdSrc['band'] == filt)[0]\n",
    "    ax.plot(FrcdSrc['expMidptMJD'][fx], FrcdSrc['psfDiffFlux'][fx],\n",
    "            plot_filter_symbols[filt], ms=10, mew=2, mec=plot_filter_colors[filt],\n",
    "            alpha=0.5, color='none', label=filt)\n",
    "    del fx\n",
    "\n",
    "for i in ix:\n",
    "    ax.plot(FrcdSrc['expMidptMJD'][ex[i]], FrcdSrc['psfDiffFlux'][ex[i]],\n",
    "            plot_filter_symbols['z'], ms=14, mew=2, mec='black', alpha=1.0,\n",
    "            color='yellow', label='z-band >3$\\sigma$ FrcdSrc psfDiffFlux')\n",
    "\n",
    "ax.hlines(y=med, xmin=numpy.min(FrcdSrc['expMidptMJD']), xmax=numpy.max(FrcdSrc['expMidptMJD']),\n",
    "          linewidth=2, color='black', ls='--', label='median z-band FrcdSrc psfDiffFlux')\n",
    "ax.axhspan(med, med+3*std, alpha=0.3, color='tab:red', label='stddev of z-band FrcdSrc psfDiffFlux')\n",
    "\n",
    "ax.set_xlabel('Modified Julian Date')\n",
    "ax.set_ylabel('Difference-Image Flux')\n",
    "ax.set_title('Forced Photometry (ForcedSourceOnDiaObject)')\n",
    "ax.legend(loc='upper left', ncol=2)\n",
    "ax.set_xlim(mjd_full_start, mjd_full_end)\n",
    "ax.axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7908d9de-27fd-44ac-b2a1-e57a970da125",
   "metadata": {},
   "source": [
    "> Figure 13: The same lightcurve as in the left panel of Figure 11, but with statistics for the $z$-band forced flux: median (dashed line) and standard deviation (pink)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a79e66e-9e2b-4c13-b6bc-18260e203d18",
   "metadata": {},
   "source": [
    "## 5. Exercises for the learner\n",
    "\n",
    "From here, there are a variety of exercises that a learner could undertake.\n",
    "\n",
    "1. Change the query constraints and identify a sample of high-redshift SNIa.\n",
    "2. Tighten the query contraints to identify only SNIa with pre-peak observations.\n",
    "3. Add error bars to the lightcurves. Magnitude errors can be retrieved during a TAP survey with, e.g., `scisql_nanojanskyToAbMagSigma(psFlux, psFluxErr) as psAbMagErr`.\n",
    "4. Apply a lightcurve-template fitter, or try a photometric classifier on the sample.\n",
    "5. Investigate the forced source photometry lightcurves of other SNIa candidates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2916ed47-1f7d-421d-a90c-27f70ef7b581",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "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.11.9"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}