DP02_08_Truth_Tables.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e9a2fea9",
   "metadata": {},
   "source": [
    "# Comparing Object and Truth Tables\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: Jeff Carlin <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: beginner <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "849eee07-7215-47bd-b6c4-e9579d46981c",
   "metadata": {},
   "source": [
    "**Description:** An introduction to using the truth data for the Dark Energy Science Collaboration's DC2 data set, which formed the basis for the DP0 data products."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3ca9f7d-1627-451a-bea5-6b524f568606",
   "metadata": {},
   "source": [
    "**Skills:** Use the TAP service with table joins to retrieve truth data matched to the Object catalog."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f60def38-6678-4c7e-9099-b7f803ca8d01",
   "metadata": {},
   "source": [
    "**LSST Data Products:** TAP dp02_dc2_catalogs.Object, .MatchesTruth, and .TruthSummary tables. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "207b4c9c-ac56-4750-a4b6-cba3b284e73b",
   "metadata": {},
   "source": [
    "**Packages:** lsst.rsp.get_tap_service, lsst.rsp.retrieve_query"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93f24e59-8555-4157-aef8-e90e3f822474",
   "metadata": {},
   "source": [
    "**Credit:** Originally developed by Jeff Carlin and the Rubin Community Science Team in the context of the Rubin DP0, with some help from Melissa Graham in the update from DP0.1 to DP0.2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c413b992-86c5-4581-89c7-b01bcd238f95",
   "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": "1b181bce",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1.0. Introduction\n",
    "\n",
    "This tutorial demonstrates how to use the TAP service to query and retrieve data from the two truth tables and the Object table for DP0.2.\n",
    "Joining these three tables enables users to compare the recovered (measured) properties (e.g., fluxes, positions, magnitudes, etc. from the Object table, which is comprised of SNR>5 detections in the deepCoadds) to the simulated values that were assigned to each object when creating the DC2 simulations.\n",
    "\n",
    "> **Thanks to the DESC!** The DC2 simulations which make up DP0 were generated by the <a href=\"https://lsstdesc.org/\">Dark Energy Science Collaboration</a> for their Data Challenge 2 (DC2).\n",
    "A full description of the simulated truth data can be found in the ApJS paper <a href=\"https://ui.adsabs.harvard.edu/abs/2021ApJS..253...31L/abstract\">The LSST DESC DC2 Simulated Sky Survey</a>, with more information in the DESC's [DC2 Data Release Note](https://ui.adsabs.harvard.edu/abs/2021arXiv210104855L/abstract).\n",
    "\n",
    "### 1.1. Package imports\n",
    "\n",
    "The [`matplotlib`](https://matplotlib.org/), [`numpy`](http://www.numpy.org/), and [`pandas`](https://pandas.pydata.org/docs/) libraries are widely used Python libraries for plotting, scientific computing, and astronomical data analysis. We will use these packages below, including the `matplotlib.pyplot` plotting sublibrary.\n",
    "\n",
    "We also use the `lsst.rsp` package to access the TAP service and query the DP0 catalogs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a75a320a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas\n",
    "from lsst.rsp import get_tap_service, retrieve_query"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "379c34ed",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and parameters\n",
    "\n",
    "Set the pandas parameter for the maximum number of rows to display to 200."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38e9a72f-ef35-4904-b4a5-e84a1039acb9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "pandas.set_option('display.max_rows', 200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c82043b-8712-4867-83c9-77c560f46ba4",
   "metadata": {},
   "source": [
    "Set `matplotlib` to show plots inline, within the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72b7fc50-681b-427d-94f5-792249424076",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "308cb6ea-be77-4129-92d5-344ebbb431cb",
   "metadata": {},
   "source": [
    "Set up colors and plot symbols corresponding to the _ugrizy_ bands. These colors are the same as those used for _ugrizy_ bands in Dark Energy Survey (DES) publications, and are defined in <a href=\"https://github.com/DarkEnergySurvey/descolors\">this github repository</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4beb776-9bc1-4485-bb63-f39a16d23792",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "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": "261349f3-12d5-4332-bd6d-f4d1fe78c66e",
   "metadata": {},
   "source": [
    "To access tables, we will use the TAP service in a similar manner to what we showed in the \n",
    "DP0.2 introductory tutorial notebook,\n",
    "and explored further in the DP0.2 tutorial notebook 02 catalog queries with TAP.\n",
    "See those notebooks for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5803048d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6cc2670",
   "metadata": {},
   "source": [
    "## 2.0. Discover truth data\n",
    "\n",
    "The <a href=\"https://dp0-2.lsst.io\">DP0.2 Documentation</a> contains a <a href=\"https://dp0-2.lsst.io/data-products-dp0-2/index.html#catalogs\">list of all DP0.2 catalogs</a>, and also a link to the <a href=\"https://dm.lsst.org/sdm_schemas/browser/dp02.html\">DP0.2 Schema Browser</a> where users can read about the available tables and their contents.\n",
    "\n",
    "Alternatively, the Portal Aspect of the Rubin Science Platform can be used to browse catalog data.\n",
    "\n",
    "Below, we show how to browse catalog data from a Notebook using the TAP service.\n",
    "\n",
    "### 2.1. Print the names of all available tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc7d10ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT description, table_name FROM TAP_SCHEMA.tables\")\n",
    "results_tab = results.to_table()\n",
    "\n",
    "for tablename in results_tab['table_name']:\n",
    "    print(tablename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79e632c8-092f-479e-9686-dd6f04a20c7e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del results, results_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca9b3dc9",
   "metadata": {},
   "source": [
    "### 2.2. Print the table schema for MatchesTruth\n",
    "\n",
    "Use the `.to_pandas()` method, and not just `.to_table()` (astropy table), so that all rows of the second cell display."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63c516fd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT column_name, datatype, description,\\\n",
    "                          unit from TAP_SCHEMA.columns\\\n",
    "                          WHERE table_name = 'dp02_dc2_catalogs.MatchesTruth'\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b88d1e0c-6878-4039-8312-86f5f337c6b3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f58edc05-fd41-46dd-9e9e-94b484401aaa",
   "metadata": {},
   "source": [
    "The above is fine if the full description is not needed, but there are some important details that are being hidden by the line truncation above.\n",
    "\n",
    "Try this instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2420334-1db2-4b47-9667-47ab1687b3cd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for c, columnname in enumerate(results['column_name']):\n",
    "    print('%-25s %-200s' % (columnname, results['description'][c]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fc304a7-dc1a-4da1-a310-2f9a88906bf4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83d7cf85-30e4-48fd-bf7a-eb4fb49af9b9",
   "metadata": {},
   "source": [
    "### 2.2. Print the table schema for TruthSummary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31d870ea",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT column_name, datatype, description,\\\n",
    "                          unit from TAP_SCHEMA.columns\\\n",
    "                          WHERE table_name = 'dp02_dc2_catalogs.TruthSummary'\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48c2b52d-6b66-4237-80f0-e8e1773ebfb2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "074e3ce7-1a34-4bff-b317-3fe643ec751b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for c, columnname in enumerate(results['column_name']):\n",
    "    print('%-25s %-200s' % (columnname, results['description'][c]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a454949d-d0cc-4b3c-a2c7-aa8a637682ee",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e2c51b5-e046-4772-ad0c-6633bc7254bc",
   "metadata": {},
   "source": [
    "## 3.0. Retrieve truth data\n",
    "\n",
    "### 3.1. Join MatchesTruth and TruthSummary\n",
    "\n",
    "As described in the column description, the column `id_truth_type` should be used to join the `MatchesTruth` and `TruthSummary` tables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58227ef0-b168-49e8-8888-4de08d318b09",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT mt.id_truth_type, mt.match_objectId, ts.ra, ts.dec, ts.truth_type \"\\\n",
    "        \"FROM dp02_dc2_catalogs.MatchesTruth AS mt \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.TruthSummary AS ts ON mt.id_truth_type = ts.id_truth_type \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', ts.ra, ts.dec), CIRCLE('ICRS', 62.0, -37.0, 0.10)) = 1 \"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fcf4839-102f-4aed-9a4a-51751a9e699b",
   "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": "0e6cf56f-3259-4daf-8489-f3e96fc41e96",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f26a422-e3bc-4e7d-b304-a37e02fb5db1",
   "metadata": {},
   "source": [
    "Notice that not all objects from the truth table have matches in the dp02_dc2_catalogs.Object table of detections in the deep coadded images (i.e., their \"match_objectId\" is blank, meaning there was no match). Print the fraction of retrieved truth objects that are matched to the dp02_dc2_catalogs.Object table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df8dfec5-31be-440d-982b-bd7674f7ee3a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tx = np.where(results['match_objectId'] > 1)[0]\n",
    "print('Number: ', len(tx))\n",
    "print('Fraction: ', np.round(len(tx)/len(results),2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad4591ce-dd10-42b8-a750-885bed08d5fa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24a99585-2c73-4a8b-aa4b-dfaada2b227a",
   "metadata": {},
   "source": [
    "### 3.2. Triple-join MatchesTruth, TruthSummary, and Objects\n",
    "\n",
    "The `MatchesTruth` table provides identifying information to pick out objects that have matches to truth objects. In order to compare _measured_ quantities (e.g., fluxes) from the `Object` table to the simulated truth values, we also need data from the `TruthSummary` table. This requires joining all three tables.\n",
    "\n",
    "> **Notice:** When restricting a JOIN query to just a portion of the sky (e.g., a cone search), restricting on coordinates in `Object`, `DiaObject`, or `Source` where possible can result in significantly improved performance.\n",
    "\n",
    "Note that the table length is now equal to the number of matched objects printed above (14850)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25bb2f05-9f85-4290-a999-1186bb3b5b5f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT mt.id_truth_type, mt.match_objectId, ts.ra, ts.dec, ts.truth_type, \"\\\n",
    "        \"obj.coord_ra, obj.coord_dec \"\\\n",
    "        \"FROM dp02_dc2_catalogs.MatchesTruth AS mt \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.TruthSummary AS ts ON mt.id_truth_type = ts.id_truth_type \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.Object AS obj ON mt.match_objectId = obj.objectId \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', obj.coord_ra, obj.coord_dec), CIRCLE('ICRS', 62.0, -37.0, 0.10)) = 1 \"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dce2bf26-2188-4396-b22a-a1bad5f3e5d4",
   "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": "377af263-8eca-437a-b864-fe78a4710c4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37264a3c-0da4-4ee6-b188-d1bddd2a0f9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af9db6c1-8a1b-49cf-8c63-c2e377f8c0f7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9524c446-e669-40b4-9140-16893d0503af",
   "metadata": {},
   "source": [
    "### 3.3. Efficiently return truth matched data for a single Object\n",
    "\n",
    "It may be a common query to check whether truth matched data exists for a given Object that was detected in a deepCoadd.\n",
    "\n",
    "> **Warning:** Please note that the restriction for the given Object is written in the query below specifically as `WHERE obj.objectId=1486698050427598336`. If we were to write `WHERE mt.match_objectId=1486698050427598336` instead, **the query could take orders of magnitude longer to execute**.\n",
    "\n",
    "This potential stumbling block stems from how the DP0 tables are stored in <a href=\"https://qserv.lsst.io/\">Qserv</a> (the distributed database): the `TruthSummary` and `Object` tables are stored in Qserv as what are known as \"director\" tables, while the `MatchesTruth` table used to join between them is stored as a somewhat more restricted \"ref match\" table.  Qserv has special mechanics to optimize queries with `WHERE` restrictions expressed in terms of director tables, and can often dispatch these queries to just a few involved data shards.  These same mechanics, however, cannot be applied in general to ref match tables, so the seemingly same restriction, if expressed in terms of the _ref match_ table, would necessitate a full scan of the entire catalog which could be quite time-consuming.\n",
    "\n",
    "While we would like Qserv to be able to notice and rewrite such queries on its own so our users would not have to be made aware of such details, this is the situation as it exists today, for DP0.  In general:\n",
    "\n",
    "> **Advisory:** If ever a simple query that seems it should take seconds takes an hour (or even ten minutes!) instead of seconds, please submit a <a href=\"https://github.com/rubin-dp0/Support\">GitHub Issue</a> or post in the <a href=\"https://community.lsst.org/c/support/dp0-rsp-issues/52\">DP0 RSP Service Issues</a> category of the Rubin Community Forum.  Staff there with specific familiarity with Qserv are happy to engage to investigate and help tweak queries for optimal execution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5300d283-6b2a-4b60-8f94-c4acd6b341da",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT mt.id_truth_type AS mt_id_truth_type, \"\\\n",
    "        \"mt.match_objectId AS mt_match_objectId, \"\\\n",
    "        \"obj.objectId AS obj_objectId, \"\\\n",
    "        \"ts.redshift AS ts_redshift \"\\\n",
    "        \"FROM dp02_dc2_catalogs.MatchesTruth AS mt \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.TruthSummary AS ts \"\\\n",
    "        \"ON mt.id_truth_type=ts.id_truth_type \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.Object AS obj \"\\\n",
    "        \"ON mt.match_objectId=obj.objectId \"\\\n",
    "        \"WHERE obj.objectId=1486698050427598336 \"\\\n",
    "        \"AND ts.truth_type=1 \"\\\n",
    "        \"AND obj.detect_isPrimary=1 \"\\\n",
    "        \"ORDER BY obj_objectId DESC\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c85a1180-0101-4bec-ad8a-d2b2aacec858",
   "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": "5c12473b-ebc0-46e4-88d5-09a25e4e94bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b06cd9b-e82b-41d9-b877-dd6b84895776",
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "415dd1e3-b2c7-4660-8696-3da4e56dff15",
   "metadata": {},
   "outputs": [],
   "source": [
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc36e8cf-9294-4904-8c3b-5714abb9d3ec",
   "metadata": {},
   "source": [
    "### 3.4. Retrieve additional data for true galaxies that are matched to detected objects\n",
    "\n",
    "With the query below we retrieve much more data for true galaxies, such as their true and measured fluxes and extendedness. Since we are only retrieving measurement data from the Object table for true galaxies, we retrieve the `cModelFlux` instead of the `psfFlux`, the latter being appropriate for point-sources (e.g., stars).\n",
    "\n",
    "> **Notice:** Above, the column names in the retrieved results have no provenance: once the data are in the results table, it is unclear from which table it originated (i.e., whether it is a true coordinate or a measured coordinate). Below, we use the `AS` statement to rename columns to start with their origin table, in order to keep track of what is from a truth table (mt_ and ts_) and what is from the object table (obj_).\n",
    "\n",
    "> **Notice:** Below, we use `truth_type = 1` to only retrieve truth and measurement data for \"true galaxies.\"\n",
    "\n",
    "The following query will retrieve 14501 results, 349 fewer than the query above, due to the specification of `truth_type = 1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "990789e8-0581-4c43-ac62-bf9cc00424b8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT mt.id_truth_type AS mt_id_truth_type, \"\\\n",
    "        \"mt.match_objectId AS mt_match_objectId, \"\\\n",
    "        \"ts.ra AS ts_ra, \"\\\n",
    "        \"ts.dec AS ts_dec, \"\\\n",
    "        \"ts.truth_type AS ts_truth_type, \"\\\n",
    "        \"ts.mag_r AS ts_mag_r, \"\\\n",
    "        \"ts.is_pointsource AS ts_is_pointsource, \"\\\n",
    "        \"ts.redshift AS ts_redshift, \"\\\n",
    "        \"ts.flux_u AS ts_flux_u, \"\\\n",
    "        \"ts.flux_g AS ts_flux_g, \"\\\n",
    "        \"ts.flux_r AS ts_flux_r, \"\\\n",
    "        \"ts.flux_i AS ts_flux_i, \"\\\n",
    "        \"ts.flux_z AS ts_flux_z, \"\\\n",
    "        \"ts.flux_y AS ts_flux_y, \"\\\n",
    "        \"obj.coord_ra AS obj_coord_ra, \"\\\n",
    "        \"obj.coord_dec AS obj_coord_dec, \"\\\n",
    "        \"obj.refExtendedness AS obj_refExtendedness, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(obj.r_cModelFlux) AS obj_cModelMag_r, \"\\\n",
    "        \"obj.u_cModelFlux AS obj_u_cModelFlux, \"\\\n",
    "        \"obj.g_cModelFlux AS obj_g_cModelFlux, \"\\\n",
    "        \"obj.r_cModelFlux AS obj_r_cModelFlux, \"\\\n",
    "        \"obj.i_cModelFlux AS obj_i_cModelFlux, \"\\\n",
    "        \"obj.z_cModelFlux AS obj_z_cModelFlux, \"\\\n",
    "        \"obj.y_cModelFlux AS obj_y_cModelFlux \"\\\n",
    "        \"FROM dp02_dc2_catalogs.MatchesTruth AS mt \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.TruthSummary AS ts ON mt.id_truth_type = ts.id_truth_type \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.Object AS obj ON mt.match_objectId = obj.objectId \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', obj.coord_ra, obj.coord_dec), CIRCLE('ICRS', 62.0, -37.0, 0.10)) = 1 \"\\\n",
    "        \"AND ts.truth_type = 1 \"\\\n",
    "        \"AND obj.detect_isPrimary = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1ee3fa8-4ed0-48fc-8e79-150680c7b51a",
   "metadata": {},
   "source": [
    "This query might take a couple of minutes to execute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06f7ad35-8acd-4687-b34d-c16b1fbeeb63",
   "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": "2e182996-da71-48db-aa2c-12694db8b84d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0797e739-a009-4679-9b2b-f3e783caba72",
   "metadata": {},
   "source": [
    "Notice that there is no `del results` statement here.\n",
    "\n",
    "Keep these results and use them below, in Section 4, to explore the retrieved data for true galaxies.\n",
    "\n",
    "## 4.0. Compare true and measured properties for true galaxies\n",
    "\n",
    "### 4.1. Plot coordinate offsets for true galaxies\n",
    "\n",
    "Below, plot the difference between the true and measured declination versus the difference between the true and measured right ascension. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "104605c3-a495-44a3-9ee3-720e87816fa4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 4))\n",
    "plt.plot(3600*(results['ts_ra']-results['obj_coord_ra']), \\\n",
    "         3600*(results['ts_dec']-results['obj_coord_dec']), \\\n",
    "         'o', ms=2, alpha=0.2, mew=0)\n",
    "plt.xlabel('Right Ascension (true-measured; [\"])', fontsize=12)\n",
    "plt.ylabel('Declination (true-measured; [\"])', fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e56bdbe6-d845-4f04-975c-9c6c43875f77",
   "metadata": {},
   "source": [
    "> Figure 1: A scatter plot of the difference between the true and measured coordinates of right ascension (RA) versus declination (Dec), in arcseconds.\n",
    "\n",
    "We see that the scatter is less than about 0.5 arcseconds. For the (DC2-simulated) LSST Science Camera's platescale of 0.2 arcsec per pixel, that's a measurement accuracy of 2.5 pixels. Note also that most galaxies' positions are measured to sub-pixel accuracy.\n",
    "\n",
    "### 4.2. How many true galaxies are measured as point sources?\n",
    "\n",
    "The number of true galaxies that are measured as point sources, or in other words, have a measured `refExtendedness` equal to zero (see the [schema for the DP0.2 Object table](https://dm.lsst.org/sdm_schemas/browser/dp02.html#Object) for more info about the `refExtendedness` column).\n",
    "\n",
    "In this example, the following cell will report that 19% of true galaxies appear as point sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98171043-8a71-4a84-a784-45f208c0566e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "x = np.where(results['obj_refExtendedness'] == 0)[0]\n",
    "print('Number: ', len(x))\n",
    "print('Fraction: ', np.round(len(x)/len(results['ts_is_pointsource']),2))\n",
    "del x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d25d690-c3c0-4a7b-9462-5617806dc390",
   "metadata": {},
   "source": [
    "The fact that 19% of the true galaxies retrieved from the catalog appear point-like does not necessarily indicate an error in the measurement pipelines. For example, very small or very distant galaxies can appear point-like, even if they were simulated as extended objects. It is left as an exercise for the learner to explore what types of galaxies are measured to be point-like."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f18cf63-e43c-4d77-94bc-95e0d9605424",
   "metadata": {},
   "source": [
    "### 4.3. Compare true and measured r-band magnitudes for true galaxies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "397e8863-e238-422b-bb6f-54ec5426cee7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 4))\n",
    "plt.plot([18,32], [18,32], ls='solid', color='black', alpha=0.5)\n",
    "x = np.where(results['obj_refExtendedness'] == 1)[0]\n",
    "plt.plot(results['ts_mag_r'][x], results['obj_cModelMag_r'][x], \\\n",
    "         'o', ms=4, alpha=0.2, mew=0, color=plot_filter_colors['r'],\\\n",
    "         label='measured as extended')\n",
    "del x\n",
    "x = np.where(results['obj_refExtendedness'] == 0)[0]\n",
    "plt.plot(results['ts_mag_r'][x], results['obj_cModelMag_r'][x], \\\n",
    "         'o', ms=2, alpha=0.5, mew=0, color='black',\\\n",
    "         label='measured as point-like')\n",
    "del x\n",
    "plt.xlabel('true r-band magnitude', fontsize=12)\n",
    "plt.ylabel('measured cModel r-band magnitude', fontsize=12)\n",
    "plt.legend(loc='lower right')\n",
    "plt.xlim([18,30])\n",
    "plt.ylim([18,30])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81d6b37c-8937-474a-9c93-5991b2d78a60",
   "metadata": {},
   "source": [
    "> Figure 2: The true $r$-band magnitude versus the measured `cModel` $r$-band magnitude for extended (dark red) and point-like (black) objects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59c13f2d-3b60-4601-9d3c-89508a437ae0",
   "metadata": {},
   "source": [
    "### 4.4. Compare true and measured fluxes in all filters for true galaxies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5240a127-72d4-4453-b9a8-66b2d7c52fbf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 3, figsize=(10, 7))\n",
    "i=0\n",
    "j=0\n",
    "for f,filt in enumerate(plot_filter_labels):\n",
    "    ax[i,j].plot([0.1,1e6], [0.1,1e6], ls='solid', color='black', alpha=0.5)\n",
    "    ax[i,j].plot(results['ts_flux_'+filt], results['obj_'+filt+'_cModelFlux'], \\\n",
    "                 plot_filter_symbols[filt], color=plot_filter_colors[filt], \\\n",
    "                 alpha=0.1, mew=0, label=filt)\n",
    "    ax[i,j].loglog()\n",
    "    ax[i,j].text(0.1, 0.9, filt, horizontalalignment='center', verticalalignment='center',\n",
    "                 transform = ax[i,j].transAxes, color=plot_filter_colors[filt], fontsize=14)\n",
    "    ax[i,j].set_xlim([0.1,1e6])\n",
    "    ax[i,j].set_ylim([0.1,1e6])\n",
    "    j += 1\n",
    "    if j == 3:\n",
    "        i += 1\n",
    "        j = 0\n",
    "ax[0,0].set_ylabel('measured cModelFlux', fontsize=12)\n",
    "ax[1,0].set_ylabel('measured cModelFlux', fontsize=12)\n",
    "ax[1,0].set_xlabel('true flux', fontsize=12)\n",
    "ax[1,1].set_xlabel('true flux', fontsize=12)\n",
    "ax[1,2].set_xlabel('true flux', fontsize=12)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca1a9c06-4ad2-4ed3-ae25-56f3eceaa94a",
   "metadata": {},
   "source": [
    "> Figure 3: The true flux versus the measured `cModel` flux for extended objects (galaxies) by filter."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25af1fce-40e8-4f4d-b88e-9adeae02dbaa",
   "metadata": {},
   "source": [
    "### 4.5. Compare color-magnitude diagrams (CMDs) for true and measured properties of true galaxies\n",
    "\n",
    "The following cells plot the true CMD at left in black, and the measured CMD at right in grey.\n",
    "\n",
    "The first pair of plots uses the r-band for magnitude, and the g-r color. The second pair of plots uses the r-band for magnitude, and i-z for color.\n",
    "\n",
    "In the first set of plots, the effects of measurement uncertainties are correlated between the _x_ and _y_ axes because the r-band data are included in both axes. In the second set of plots, the i-band and the z-band are instead used for color. Notice how the effect of measurement uncertainties changes.\n",
    "\n",
    "Recall that these plots do not contain data for stars, as only true galaxies were retrieved from the truth tables.\n",
    "\n",
    "> **Warning:** Pink \"RuntimeWarning\" errors will appear due to a few of the measured fluxes in the denominator being zero. It is OK to ignore these warnings for the context of this tutorial, which focuses on retrieving truth data, but for scientific analyses users should follow up and understand such warnings (e.g., use flags to reject poor flux measurements from their samples)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffe668ed-bfe5-4ed0-9cfb-e3a4d198c397",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax[0].plot(-2.5*np.log10(results['ts_flux_g']/results['ts_flux_r']), results['ts_mag_r'], \\\n",
    "           'o', ms=2, alpha=0.2, mew=0, color='black')\n",
    "\n",
    "ax[1].plot(-2.5*np.log10(results['obj_g_cModelFlux']/results['obj_r_cModelFlux']), results['obj_cModelMag_r'], \\\n",
    "           'o', ms=2, alpha=0.2, mew=0, color='grey')\n",
    "ax[0].set_xlabel('true color (g-r)', fontsize=12)\n",
    "ax[0].set_ylabel('true magnitude (r-band)', fontsize=12)\n",
    "ax[0].set_xlim([-2, 4])\n",
    "ax[0].set_ylim([30, 18])\n",
    "ax[1].set_xlabel('measured color (g-r)', fontsize=12)\n",
    "ax[1].set_ylabel('measured magnitude (r-band)', fontsize=12)\n",
    "ax[1].set_xlim([-2, 4])\n",
    "ax[1].set_ylim([30, 18])\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f509166-2832-481a-991c-977fe1f6455e",
   "metadata": {},
   "source": [
    "> Figure 4: At left, the true $g-r$ color versus true $r$-band magnitude. At right, the measured $g-r$ color versus measured $r$-band magnitude. This plot demonstrates the impact of measurement errors on galaxy colors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95bc8957-80ce-4448-8ed0-d2ecf1769a0b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax[0].plot(-2.5*np.log10(results['ts_flux_i']/results['ts_flux_z']), results['ts_mag_r'], \\\n",
    "           'o', ms=2, alpha=0.2, mew=0, color='black')\n",
    "\n",
    "ax[1].plot(-2.5*np.log10(results['obj_i_cModelFlux']/results['obj_z_cModelFlux']), results['obj_cModelMag_r'], \\\n",
    "           'o', ms=2, alpha=0.2, mew=0, color='grey')\n",
    "ax[0].set_xlabel('true color (i-z)', fontsize=12)\n",
    "ax[0].set_ylabel('true magnitude (r-band)', fontsize=12)\n",
    "ax[0].set_xlim([-2, 4])\n",
    "ax[0].set_ylim([30, 18])\n",
    "ax[1].set_xlabel('measured color (i-z)', fontsize=12)\n",
    "ax[1].set_ylabel('measured magnitude (r-band)', fontsize=12)\n",
    "ax[1].set_xlim([-2, 4])\n",
    "ax[1].set_ylim([30, 18])\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60041aa8-2cce-4ea5-a368-ff03929ffc42",
   "metadata": {},
   "source": [
    "> Figure 5: Similar to Figure 4, but with $i-z$ color instead of $g-r$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c28ee04d-07b6-46a3-b7ab-50e7edd4eaa2",
   "metadata": {},
   "source": [
    "## 5.0 Exercises for the learner\n",
    "\n",
    "1. Repeat the query in Section 3.3, but instead of only retrieving true galaxies (`ts.truth_type = 1`), include stars, which have a `truth_type` of 2. Since stars are point sources, instead of only retrieving `cModelFlux` measurements, also retrieve `psfFlux` from the Object catalog, because PSF-fit fluxes are more appropriate for point sources.\n",
    "\n",
    "2. As mentioned in Section 4.2, it is left as an exercise for the learner to explore what types of galaxies are measured to be point-like.\n",
    "\n",
    "3. Explore the truth data for Type Ia supernovae (`truth_type = 3`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5d23342-9904-4b6f-ad99-a228547dde57",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}