DP02_12b_PSF_Science_Demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "# Point Spread Function (PSF) Analysis for Cosmology\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(s): Andrés A. Plazas Malagón<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: advanced <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** Demonstration of the use of measured PSF properties in weak lensing analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Use of the catalog data products for PSF analysis. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** DP0.2 `Object` and `Source` catalogs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c67fab9-136a-4adc-bb42-142b91ab69dd",
   "metadata": {},
   "source": [
    "**Packages:** lsst.analysis.tools.atools"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "metadata": {},
   "source": [
    "**Credits:** Developed by Andrés A. Plazas Malagón in collaboration with Melissa Graham, Jeff Carlin, Douglas Tucker, Ryan Lau, and the Rubin Community Science Team for DP0.2. The section on PSF size and its correlation function using the `treecorr` software is based on a Jupyter notebook authored by Claire-Alice Hebert. The section on PSF size histograms and their dependency on the central filter wavelength is derived from a Jupyter Notebook by Tianqing Zhang. The section on `rho` statistics using `lsst.analysis.tools` draws from a Jupyter Notebook by Arun Kannawadi.\n",
    "This notebook incorporates suggestions from the [Accessible Authoring Checklist](https://iota-school.github.io/accessibility_hackathon/hack#checklist) and utilizes NASA's Astrophysics Data System Bibliographic Services."
   ]
  },
  {
   "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": {
    "tags": []
   },
   "source": [
    "## 1. Introduction\n",
    "\n",
    "In the context of astronomical surveys, such as the Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST), \n",
    "a significant objective is to extract cosmological information by scrutinizing subtle deviations in the \n",
    "shapes and orientations of galaxies. \n",
    "These deviations are caused by gravitational lensing from large-scale structures in the foreground. \n",
    "However, the observed shapes of galaxies are not solely influenced by the gravitational lensing effects \n",
    "of the universe's large-scale structure (referred to as \"cosmic shear\"), \n",
    "but also by a combined blurring effect known as the point spread function (PSF). \n",
    "The PSF arises from various factors such as atmospheric conditions (for ground-based instruments), \n",
    "telescope optics, and the image sensor.\n",
    "PSF analysis is thus an integral part of cosmological studies (and is important in other astronomical fields as well).\n",
    "\n",
    "Whereas the introductory DP0.2 PSF tutorial notebook (DP02_12a_PSF_Data_Products) demonstrated how to access\n",
    "PSF information for `calexp` and `deepCoadd` images, this tutorial focuses on the measured PSF information\n",
    "that is available in the `Source` and `Object` catalogs for sources detected in the `calexp` and `deepCoadd`\n",
    "images, respectively.\n",
    "\n",
    "**Section 2** introduces and explores the PSF measurements that are available in the `Object` and `Source` catalogs.\n",
    "Two TAP queries are executed: one to retrieve a set of `Objects` detected in the `deepCoadd` images,\n",
    "and one to retrieve a set of `Sources` detected in a `calexp` (single visit) image at the same location.\n",
    "\n",
    "**Section 3** calculates and displays the size of the PSF in a particular photometric band (`i`) as a function of sky location. \n",
    "It also computes the size-size two-point correlation function (autocorrelation) using the `treecorr` software. \n",
    "This software is widely used in state-of-the-art weak lensing analyses to derive cosmological parameters. \n",
    "These calculations are performed using both `calexp` and `deepCoadd` images.\n",
    "\n",
    "**Section 4** illustrates one chromatic dependence of the PSF: the wavelength dependence of the PSF size.\n",
    "\n",
    "**Section 5** calculates the rho statistics, a series of two-point correlation functions used to diagnose \n",
    "the quality of a PSF model for weak lensing analyses. This calculation is done using the `lsst.analysis.tools` framework."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de0dbefc-e2f3-4ec2-b49c-318ad5132417",
   "metadata": {},
   "source": [
    "**Additional resources.** \n",
    "\n",
    "It is not required to review and understand the following resources before working through this tutorial.\n",
    "\n",
    "For an overall review on PSF modelling and a look towards the future, see\n",
    "\"Point spread function modelling for astronomical telescopes: a review focused on weak gravitational lensing studies\" by\n",
    "[Liaudat et al. 2023](https://ui.adsabs.harvard.edu/abs/2023arXiv230607996L/abstract).\n",
    "For papers on the impact of chromatic effects on the PSF on weak lensing cosmology see \n",
    "\"Atmospheric Dispersion Effects in Weak Lensing Measurements\" by \n",
    "[Plazas and Bernstein 2012](https://ui.adsabs.harvard.edu/abs/2012PASP..124.1113P/abstract),\n",
    "or \"Impact of Atmospheric Chromatic Effects on Weak Lensing Measurements\" by\n",
    "[Meyers and Burchat 2015](https://ui.adsabs.harvard.edu/abs/2015ApJ...807..182M/abstract).\n",
    "For a paper on the impact of differential chromatic refraction (DCR) on supernova cosmlogy, see\n",
    "\"The Dark Energy Survey Supernova Program: Corrections on Photometry Due to Wavelength-dependent Atmospheric Effects\" by\n",
    "[Lee et al. 2023](https://ui.adsabs.harvard.edu/abs/2023AJ....165..222L/abstract).\n",
    "\n",
    "For more details on PSF modeling and null tests for galaxy surveys, including `rho` statistics, see\n",
    "\"Dark Energy Survey year 3 results: point spread function modelling\" by \n",
    "[Jarvis et al. 2021](https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.1282J/abstract),\n",
    "or\n",
    "\"The three-year shear catalog of the Subaru Hyper Suprime-Cam SSP Survey\" by\n",
    "[Li et al. 2022](https://ui.adsabs.harvard.edu/abs/2022PASJ...74..421L/abstract).\n",
    "And finally, for an approach to a mathematical framework that aids in producing well-defined coadded PSFs for static-sky sources,\n",
    "see \"PSFs of coadded images\" by\n",
    "[Mandelbaum et al. 2022](https://ui.adsabs.harvard.edu/abs/2023OJAp....6E...5M/abstract)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc36f107",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1.1. Import packages\n",
    "\n",
    "The [matplotlib](https://matplotlib.org/) (and especially sublibrary `matplotlib.pyplot`), [numpy](http://www.numpy.org/),\n",
    "and [scipy](https://scipy.org/) libraries are widely used Python libraries for plotting and scientific computing and analysis.\n",
    "The [treecorr](https://rmjarvis.github.io/TreeCorr/_build/html/index.html) package by\n",
    "[Jarvis et al. 2004](https://ui.adsabs.harvard.edu/abs/2004MNRAS.352..338J/abstract) \n",
    "is a widely used tool for fast correlations measurements based on a ball tree method \n",
    "(similar to a [k-d tree](https://en.wikipedia.org/wiki/K-d_tree)). \n",
    "In particular, it is widely used in state-of-the-art weak lensing cosmological analyses. \n",
    "\n",
    "The `lsst.rsp` package provides access to the Table Access Protocol (TAP) service for queries to the DP0 catalogs.\n",
    "The `lsst.analysis.tools` package is in-developement by Rubin staff for the \n",
    "creation of consistent and repeatable plots and metrics.\n",
    "Intermediate and advanced users are welcome to use it too; more information is available in the \n",
    "[analysis tools documentation](https://pipelines.lsst.io/modules/lsst.analysis.tools/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d185358a-4c15-4b04-8e3d-13083191ee72",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "import treecorr\n",
    "\n",
    "from lsst.rsp import get_tap_service\n",
    "from lsst.analysis.tools.atools import RhoStatistics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa240aaf-d675-4d72-bb68-d736200f4980",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78fc9bc8-29d1-48e7-a9a5-d6b0fa1a1cfa",
   "metadata": {},
   "source": [
    "Instantiate the TAP service for DP0.2 catalog access."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c493f00-7a7d-4c8d-85c5-46e45bd4ae30",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1221e437-0e2c-4034-8af1-04fb8e9b6307",
   "metadata": {},
   "source": [
    "Set a standard color style."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "400c0d16-523f-49cb-8369-3f3d19bb8fc0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.style.use('tableau-colorblind10')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "955b0a82-9fb4-46fb-bc9b-2db4c04497ba",
   "metadata": {},
   "source": [
    "**Function**: `powerLaw`\n",
    "\n",
    "Defines a power-law function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d654417-a7d2-4a2c-9846-de3ee6bcbafb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def powerLaw(lam, a, b):\n",
    "    \"\"\"Function to return a power-law.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    lam : `array`\n",
    "       array of values\n",
    "    a : `float`\n",
    "       constant factor\n",
    "    b : `float`\n",
    "       power\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    result : `float`\n",
    "        a*lam**(-b)\n",
    "    \"\"\"\n",
    "\n",
    "    result = a*lam**(-b)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5285320-1e9f-4618-857e-dbdd50ec2f12",
   "metadata": {},
   "source": [
    "## 2. The measured PSF properties in the catalogs\n",
    "\n",
    "Recall that the `Objects` and `Sources` catalogs contains measurements for astrophysical objects \n",
    "detected in the `deepCoadd` and `calexp` (single-visit) images, respectively.\n",
    "\n",
    "These measurements include a variety of parameters that characterize the PSF as measured for\n",
    "that `Object` or `Source`.\n",
    "While these parameters are calculated for all detections, they\n",
    "are only relevant for point-sources (stars) and not extended sources (galaxies).\n",
    "\n",
    "### 2.1. PSF properties for `Objects`\n",
    "\n",
    "To explore what PSF-related measurements exist, one option is to page search for the term \"psf\" in the \n",
    "[DP0.2 Catalogs Schema Browser](https://dm.lsst.org/sdm_schemas/browser/dp02.html) for the `Object` table.\n",
    "\n",
    "**Option:** To retrieve the column names and descriptions for the `Object` table and search for the term \"psf\",\n",
    "uncomment the lines in the following cell and execute the cell.\n",
    "The result is quite a long list, so this cell is commented-out by default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6afee619-2c4a-4ce5-8b1e-550f11ff7f2c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT column_name, description \"\n",
    "#                          \"from TAP_SCHEMA.columns \"\n",
    "#                          \"WHERE table_name = 'dp02_dc2_catalogs.Object'\"\n",
    "#                          ).to_table().to_pandas()\n",
    "# for cname, descr in zip(results['column_name'], results['description']):\n",
    "#     if (cname.find('psf') > -1) | (cname.find('PSF') > -1):\n",
    "#         print(cname + ' --- ' + descr)\n",
    "# del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54a159d4-7ed6-41ce-a6e8-14e45a7924d6",
   "metadata": {
    "tags": []
   },
   "source": [
    "The `Object` catalog contains measurements for sources detected in the `deepCoadd` images for *all* six filters, *ugrizy*.\n",
    "Thus, for all the PSF-related measurements in the `Object` catalog, there is one per band.\n",
    "\n",
    "Some of the most useful PSF-related parameters in the `Object` catalog include the following.\n",
    "\n",
    "`{band}_psfFlux` and `{band}_psfFluxErr`:\n",
    "The flux (and its error) derived from linear least-squares forced-fit of PSF model.\n",
    "\n",
    "`{band}_psfFlux_flag`: A general failure flag for the linear least-squares forced-fit of PSF model (`True` or `1` if failed).\n",
    "\n",
    "`{band}_ixxPSF`, `{band}_ixyPSF`, and `{band}_iyyPSF`: Second-moments of the shape of the fit PSF model ($I_{xx}$, $I_{yy}$, and $I_{xy}$).\n",
    "\n",
    "`{band}_calib_psf_used`: A flag that is `True` (or `1`) if the `Object` is a star that was used to establish the PSF model\n",
    "for the image (i.e., bright isolated stars with well-measured PSF properties).\n",
    "\n",
    "`{band}_pixelFlags_inexact_psfCenter`: A flag that is `True` (or `1`) if the `Object` is close to image pixels which do not\n",
    "have an exact PSF, typically because they are near an image edge. \n",
    "Objects near these edges can have a discontinuous PSF."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd891c38-bca7-4413-bc90-eb88f3bdeff6",
   "metadata": {},
   "source": [
    "#### 2.1.1. Retrieve a sample of `Objects`\n",
    "\n",
    "Use the same central coordinates as were used to retrieve the images in the introductory PSF tutorial notebook\n",
    "(DP02_12a_PSF_Data_Products.ipynb).\n",
    "\n",
    "Use the same default filter, *i*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8396cdb0-cc76-4223-88bd-447561217225",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "use_center_coords = \"53.0, -34.1\"\n",
    "band = 'i'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25322224-42db-45d3-a968-4c2f4afd1c04",
   "metadata": {},
   "source": [
    "The following query retrieves the PSF-related columns mentioned above, plus a few more columns\n",
    "which will be used later in this notebook, when calculating rho statistics with the analysis tools pacakge.\n",
    "\n",
    "The `WHERE` statements for this query require object to:\n",
    " - be within a 2 degree radius around the center coordinates;\n",
    " - have a `detect_isPrimary` flag that is `True` (`1`), which is always recommended as a way to avoid duplicates and unblended objects;\n",
    " - have a flux $> 360$ nJy (i.e., be bright);\n",
    " - have an extendedness parameter $= 0$ (stars only, no extended objects);\n",
    " - have a `{band}_calib_psf_used` flag that is `True` (`1`), such that only stars with well-measured PSF properties are returned;\n",
    " - and have a `{band}_pixelFlags_inexact_psfCenter` flag that is `False` (`0`) to avoid stars near the image edge.\n",
    " \n",
    "The choice to requre that `i_pixelFlags_inexact_psfCenter = 0` is influenced by the fact that the `EDGE` \n",
    "flag in a `calexp` corresponds to the `SENSOR_EDGE` flag in the `coadd`. \n",
    "Objects near these edges have a discontinuous PSF. \n",
    "Additionally, the `inexact_psfCenter` results from an `OR` operation that includes `SENSOR_EDGE`. \n",
    "The other properties retrieved through the query will play a role in calculating the `rho` statistics, \n",
    "which will be the subject of Section 5 in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "602d9d74-e0df-4a6e-b1a8-e150998eb392",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT coord_ra, coord_dec, \"\\\n",
    "        f\"{band}_psfFlux, \"\\\n",
    "        f\"{band}_psfFluxErr, \"\\\n",
    "        f\"{band}_psfFlux_flag, \"\\\n",
    "        f\"{band}_pixelFlags_saturatedCenter, \"\\\n",
    "        f\"{band}_extendedness, \"\\\n",
    "        f\"{band}_extendedness_flag, \"\\\n",
    "        f\"{band}_ixy, {band}_ixyPSF, \"\\\n",
    "        \"u_ixx, u_iyy, u_ixxPSF, u_iyyPSF, \"\\\n",
    "        \"g_ixx, g_iyy, g_ixxPSF, g_iyyPSF, \"\\\n",
    "        \"r_ixx, r_iyy, r_ixxPSF, r_iyyPSF, \"\\\n",
    "        \"i_ixx, i_iyy, i_ixxPSF, i_iyyPSF, \"\\\n",
    "        \"xy_flag, sky_object, \"\\\n",
    "        \"detect_isPatchInner, \"\\\n",
    "        \"detect_isDeblendedSource \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \"+use_center_coords+\", 2.0)) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1 \"\\\n",
    "        f\"AND {band}_calibFlux > 360 \"\\\n",
    "        f\"AND {band}_extendedness = 0.0 \"\\\n",
    "        f\"AND {band}_calib_psf_used = 1 \"\\\n",
    "        f\"AND {band}_pixelFlags_inexact_psfCenter = 0 \"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce881c22-bc0a-46f7-b25e-bb0ebf505f15",
   "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": "0a7beba7-896b-4d4b-a299-918a36b2ec79",
   "metadata": {},
   "outputs": [],
   "source": [
    "objects = job.fetch_result().to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76f78f93-5857-48ab-b3de-30b5f1d08066",
   "metadata": {},
   "source": [
    "**Option:** View the table. This query returns 8426 rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fe9bf20-49e4-49d0-bc7b-92cc06aded58",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# objects"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "801681d1-de2b-4946-9aea-7623178dbafb",
   "metadata": {},
   "source": [
    "Calculate the trace radius (PSF size, $r_t$) for the *u*, *g*, *r*, and *i*-band \n",
    "from the second moments of the PSF shape as $r_t = \\sqrt{\\frac{I_{xx} + I_{yy}}{2}}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c895d01b-17f2-454d-9249-b327f95cd122",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for band in ['u', 'g', 'r', 'i']:\n",
    "    objects['size_'+band+'_PSF'] = np.sqrt((objects[band + '_ixxPSF']\n",
    "                                            + objects[band + '_iyyPSF']) / 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a75a4896-b20f-4e77-81dd-f478006ff3e3",
   "metadata": {},
   "source": [
    "The PSF sizes in different filters will be used below to explore the wavelength dependence of seeing."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60d3b3d1-5b9f-45a5-9206-c7c3023c91dd",
   "metadata": {},
   "source": [
    "### 2.2. PSF properties for `Sources`\n",
    "\n",
    "**Option:** To retrieve the column names and descriptions for the `Source` table and search for the term \"psf\",\n",
    "uncomment the lines in the following cell and execute the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dc04273-4311-4046-b4fc-d4ec608af297",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT column_name, description \"\n",
    "#                          \"from TAP_SCHEMA.columns \"\n",
    "#                          \"WHERE table_name = 'dp02_dc2_catalogs.Source'\"\n",
    "#                          ).to_table().to_pandas()\n",
    "# for cname, descr in zip(results['column_name'], results['description']):\n",
    "#     if (cname.find('psf') > -1) | (cname.find('PSF') > -1):\n",
    "#         print(cname + ' --- ' + descr)\n",
    "# del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55724e21-ce31-4430-bd44-e9fae424f80a",
   "metadata": {},
   "source": [
    "The `Source` catalog contains mostly the same PSF-related columns as the `Object` table.\n",
    "\n",
    "However, since `Sources` are measured only in one `calexp` (single visit) image, columns for each of the\n",
    "six filters are not needed.\n",
    "The `Sources` catalog has a column named `band` that stores the filter of the `calexp`, which is the filter\n",
    "in which all the PSF-related measurements for the `Source` have been made."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fd58914-51b9-464a-a06c-31abb7a27ae7",
   "metadata": {},
   "source": [
    "#### 2.2.1. Retrieve a sample of `Sources`\n",
    "\n",
    "The following query retrieves PSF-related columns for `Sources` detected in visit number 192350,\n",
    "which is the same visit used in the introductory PSF notebook and which is already known to be in the *i*-band\n",
    "and to overlap the sky coordinates used for the `Object` query.\n",
    "\n",
    "Similar `WHERE` constraints are used as above, but since the `Source` table does not have an \n",
    "equivalent to `{band}_pixelFlags_inexact_psfCenter` in the `Object` table,\n",
    "the constraint that the `pixelFlags_suspectCenter` flag is `False` (`0`) is used to exclude sources near\n",
    "\n",
    "The `Source` table is much larger than the `Object` table and so this query can take a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9227633d-5e0b-4169-bedf-0235c54b77de",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT coord_ra, coord_dec, \"\\\n",
    "        \"ixx, iyy, ixy, ixxPSF, iyyPSF, ixyPSF \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Source \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \"+use_center_coords+\", 2.0)) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1 \"\\\n",
    "        \"AND calibFlux > 360 \"\\\n",
    "        \"AND extendedness = 0.0 \"\\\n",
    "        \"AND calib_psf_used = 1 \"\\\n",
    "        \"AND pixelFlags_suspectCenter = 0 \"\\\n",
    "        \"AND visit = 192350\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa9942a6-e1e0-4d6b-b576-bb5efb2befa1",
   "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": "85a616a0-2487-4508-9515-990954916047",
   "metadata": {},
   "outputs": [],
   "source": [
    "sources = job.fetch_result().to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce5d2987-156b-4ba8-b772-7aa26d6e9fdb",
   "metadata": {},
   "source": [
    "**Option:** View the table. This query returns 8049 rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed2e594e-dedc-4c47-bac3-1574156e9fb5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# sources"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2bf59f9-82ee-47c2-a145-b914d757b1e5",
   "metadata": {},
   "source": [
    "Calculate the trace radius."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfa1dea8-b6e4-4846-8716-61eafdb22c4e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sources['size_PSF'] = np.sqrt((sources['ixxPSF'] + sources['iyyPSF']) / 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf7aa3d2-560e-48f5-855f-3be97a65c012",
   "metadata": {},
   "source": [
    "## 3. Visualize the PSF size and its correlation function with `treecorr`\n",
    "\n",
    "The two-point correlation functions hold significant importance in weak lensing cosmological analysis, \n",
    "as they allow for the inference of crucial cosmological parameters. \n",
    "\n",
    "Notably, `treecorr` serves as the standard tool employed by various collaborations, including the \n",
    "Dark Energy Survey and the Hyper-Suprime Camera Survey. \n",
    "It's also a prominent choice utilized by the LSST Science Pipelines code.\n",
    "\n",
    "In a broader sense, the two-point correlation function serves as a statistical metric that aids in comprehending the distribution and potential correlations of random variables within space or time.\n",
    "In cosmology, this function was introduced as a means to quantify the excess probability of locating two galaxies at a specific distance apart, as opposed to an expected uniform distribution ([Pebbles 1980](https://ui.adsabs.harvard.edu/abs/1980lssu.book.....P/abstract)). \n",
    "In modern galaxy surveys such as the Dark Energy Survey (DES), this tool is applied for cross-correlating key quantities like galaxy number counts and weak lensing shear in a \"3x2\" point analysis. \n",
    "This analysis methodology is employed to deduce constraints on cosmological parameters, and it's an approach that will also be adopted by the Legacy Survey of Space and Time (LSST) (refer to [DES Collaboration 2022](https://ui.adsabs.harvard.edu/abs/2022PhRvD.105b3520A/abstract)).\n",
    "\n",
    "Below, the `treecorr` software ([Jarvis et al. 2015](https://ui.adsabs.harvard.edu/abs/2015ascl.soft08007J/abstract)) \n",
    "is used to compute the two-point correlation function of the PSF model size as a function of angular separation.\n",
    "For more information about `treecorr`, check its [documentation](https://rmjarvis.github.io/TreeCorr/_build/html/index.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e280ddb-afe6-4336-bdbc-eb8e9206720c",
   "metadata": {},
   "source": [
    "### 3.1. For `Objects`\n",
    "\n",
    "Temporarily store the `object` right ascention and declination coordinates in the `obj_ra `and `obj_dec` variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5219e352-7f96-41c3-ab29-78835a8544cf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "obj_ra = objects['coord_ra']\n",
    "obj_dec = objects['coord_dec']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4de61d34-d0ed-44aa-a380-98a90f19ec65",
   "metadata": {},
   "source": [
    "Define `treecorr`'s catalog and adjust a few of its parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0321df4-9cbd-4d54-9830-28bf0fb1a7ab",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cat = treecorr.Catalog(ra=obj_ra, dec=obj_dec,\n",
    "                       k=objects[f'size_{band}_PSF']\n",
    "                       - np.mean(objects[f'size_{band}_PSF']),\n",
    "                       ra_units='deg', dec_units='deg')\n",
    "\n",
    "kk_config = {'max_sep': .06, 'min_sep': .0001, 'nbins': 12}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a18542a2-c4b1-4042-9b4d-750d78881c99",
   "metadata": {},
   "source": [
    "Calculate the two-point correlation function of the size. In this case, this is an auto-correlation (size-size two point correlation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0ed6dee-6240-4875-a444-89f3eb8a7b90",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "kk = treecorr.KKCorrelation(kk_config)\n",
    "kk.process(cat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b099ebcb-9640-4f59-931e-3f9ed5a76fae",
   "metadata": {},
   "source": [
    "Retrieve the correlation function, denoted by `xi`, and the angular separation arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f12e1e1-b2b1-4893-8308-655b9695943e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xi = kk.xi\n",
    "bins = kk.rnom"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13ad3d23-1a9a-47bd-aaea-003c5ad1de0d",
   "metadata": {},
   "source": [
    "Plot the PSF size as a function of sky coordinates, and the correlation function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e58bd5ff-0213-4465-a22b-06241c254d37",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "_, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.5),\n",
    "                             gridspec_kw={'wspace': .3})\n",
    "\n",
    "ax1.set_title(f'PSF size. Band : {band}')\n",
    "scatter_plot = ax1.scatter(obj_ra, obj_dec, c=objects[f'size_{band}_PSF'],\n",
    "                           s=1, cmap='cividis')\n",
    "ax1.set_xlabel('RA [deg]')\n",
    "ax1.set_ylabel('DEC [deg]')\n",
    "plt.colorbar(scatter_plot, ax=ax1, label='PSF size [pixels]')\n",
    "\n",
    "ax2.set_title(f'PSF size correlation. Band : {band}')\n",
    "ax2.plot(np.degrees(bins), xi*1e4, 'o-', color='darkblue')\n",
    "ax2.axhline(0, linestyle='--', color='lightgrey')\n",
    "ax2.set_xscale('log')\n",
    "ax2.set_ylabel(r'$\\xi \\times 10^{4}$', labelpad=10)\n",
    "ax2.set_xlabel(r'$\\theta$ [degree]')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8e07914-4f7a-48f5-8aa4-78be1833dc4b",
   "metadata": {},
   "source": [
    "> **Figure 1:** *The figure above consists of two panels. The left panel illustrates the PSF size within the investigated circle and band, as a function of declination (ranging from approximately negative 36 to negative 32 degrees) and right ascension (ranging from around 51 to 55 degrees). The PSF size or trace radius is depicted using a color scale, with values ranging from about 1.68 to 1.78 pixels. The right panel of the figure showcases the two-point autocorrelation function of the PSF size, plotted against angular separations in degrees. The correlation function values are expressed as multiples of 1e4, spanning from negative 1 to 4. Angular separations are marked on the abscissa axis, ranging from 0.01 to 2 degrees, and are represented on a logarithmic scale.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02b22d24-1ff7-4688-8938-d38092af4a96",
   "metadata": {},
   "source": [
    "### 3.2. For `Sources`\n",
    "\n",
    "Temporarily store the `sources` right ascention and declination coordinates in the `src_ra `and `src_dec` variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db013cb9-af0f-44b7-b5d8-cdbe55b189dc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "src_ra = sources['coord_ra']\n",
    "src_dec = sources['coord_dec']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "189b93a3-aacd-47cc-996f-28ead69d47d1",
   "metadata": {},
   "source": [
    "Define `treecorr`'s catalog and adjust a few of its parameters, calculate the two-point correlation function of the size, and retrieve it, along with the angular seprations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c85aeb14-e189-496d-8aa1-3cc4982983a4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cat = treecorr.Catalog(ra=src_ra, dec=src_dec,\n",
    "                       k=sources['size_PSF']\n",
    "                       - np.mean(sources['size_PSF']),\n",
    "                       ra_units='deg', dec_units='deg')\n",
    "kk_config = {'max_sep': .06, 'min_sep': .0001, 'nbins': 12}\n",
    "kk = treecorr.KKCorrelation(kk_config)\n",
    "kk.process(cat)\n",
    "xi_source = kk.xi\n",
    "bins_source = kk.rnom"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cdf2448-223e-4f65-85fc-609e0346e85d",
   "metadata": {},
   "source": [
    "Plot the PSF size as a function of sky coordinates, and the correlation function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfcdcbc8-f964-43ad-8fa3-bdc86ab1a796",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "_, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.5),\n",
    "                             gridspec_kw={'wspace': .3})\n",
    "\n",
    "ax1.set_title('PSF size. Visit : 192350, band : i')\n",
    "scatter_plot = ax1.scatter(src_ra, src_dec,\n",
    "                           c=sources['size_PSF'],\n",
    "                           s=1, cmap='cividis')\n",
    "ax1.set_xlabel('RA [deg]')\n",
    "ax1.set_ylabel('DEC [deg]')\n",
    "plt.colorbar(scatter_plot, ax=ax1, label='[pixels]')\n",
    "\n",
    "ax2.set_title('PSF size correlation. Visit : 192350, band : i')\n",
    "ax2.plot(np.degrees(bins_source), xi_source*1e4, 'o-', color='darkblue')\n",
    "ax2.axhline(0, linestyle='--', color='lightgrey')\n",
    "ax2.set_xscale('log')\n",
    "ax2.set_ylabel(r'$\\xi \\times 10^{4}$', labelpad=-12)\n",
    "ax2.set_xlabel(r'$\\theta$ [degree]')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dce500ff-0dee-46fd-962f-f6e2f399b399",
   "metadata": {},
   "source": [
    "> **Figure 2:** *The figure above consists of two panels. The left panel illustrates the PSF size for sources detected in a single visit as a function of sky coordinate, and which are within two degrees of the specified central coordinates. These coordinates are clearly near the edge of the field of view, as the outline of the CCD array edge is seen. The PSF size or trace radius is depicted using a color scale, with values ranging from about 1.59 to 1.63 pixels. The right panel of the figure showcases the two-point autocorrelation function of the PSF size, plotted against angular separations in degrees. The correlation function values are expressed as multiples of 1e4, spanning from negative 0.2 to 0.6. Angular separations are marked on the abscissa axis, ranging from 0.01 to 2 degrees, and are represented on a logarithmic scale.*\n",
    "\n",
    "\n",
    "### 3.3. Discussion \n",
    "\n",
    "By comparing Figure 1 and Figure 2, it is clear that the PSF size exhibits variations in both the `deepCoadds` (`Objects`) and in individual `calexps` (`Sources`). For the latter, the gradient observed is dominated by Kolmogorov turbulence, demonstrating the turbulent seeing variation across the focal plane. These variations are averaged down in the case of the `deepCoadds`, and therefore the correlation function decreases faster with distance. Both are in the `i` band in this case, and as a consequence, there are no chromatic variations in size coming from effects such as Differential Chromatic Refraction (DCR) or the seeing dependence with wavelength, which will be explored in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c700951-55c1-4875-ac01-78e10f115488",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 4. The wavelength dependence of seeing  \n",
    "\n",
    "To correct for the effects of the Point Spread Function (PSF) and to deduce the weak lensing signal (referred to as \"shear\") from galaxy shapes, the PSF is typically modeled using stars, often bright ones. Subsequently, the PSF of galaxy images is deconvolved by utilizing the estimated PSF convolution kernel. Implicit in this methodology is the assumption that the convolution kernel for galaxies remains the same as that for stars. However, this assumption becomes invalid if the PSF varies with wavelength due to the distinctive Spectral Energy Distributions (SEDs) of stars and galaxies, resulting in disparate PSFs.\n",
    "\n",
    "In this section, use the point source sample derived from the `Object` catalog, which was obtained in the preceding sections. Aftwerward, compute one of the wavelength-dependent components contributing to the PSF owing to atmospheric effects: the wavelength dependence on seeing. As per the standard theory of atmospheric turbulence, the linear dimension `θ` of the atmospheric convolution kernel (i.e., seeing) is inversely proportional to the wavelength, expressed as $\\theta \\propto \\lambda^{-\\frac{1}{5}}$.\n",
    "\n",
    "Additional wavelength-dependent contributions to the PSF encompass Differential Chromatic Refraction (DCR), telescope diffraction, variations in the absorption length of camera CCDs, and optical chromatic aberrations. In general, each of these factors exhibits distinct functional wavelength dependencies and may exert dominance in different spectral bands. Nevertheless, the wavelength-dependent seeing that we explore in this section tends to be the primary contributor to the PSF. These chromatic PSF effects introduce biases in observables such as photometry, astrometry, measurements of galaxy shape and size, and higher-order moments of the PSF. Addressing these effects is imperative for precise cosmological analysis employing weak gravitational lensing. For more detailed information, refer to the works of [Plazas and Bernstein 2012](https://ui.adsabs.harvard.edu/abs/2012PASP..124.1113P/abstract), [Lee et al. 2023](https://ui.adsabs.harvard.edu/abs/2023AJ....165..222L/abstract), [Meyers and Burchat 2015](https://ui.adsabs.harvard.edu/abs/2015ApJ...807..182M/abstract), [Tianqing et al. 2023](https://ui.adsabs.harvard.edu/abs/2023MNRAS.520.2328Z/abstract), and [Meyers and Burchat 2015b](https://ui.adsabs.harvard.edu/abs/2015JInst..10C6004M/abstract)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "665398c4-00a4-44f6-abd3-e3122e71937d",
   "metadata": {},
   "source": [
    "To replicate this wavelength relationship using the date retreived above, compute first the central wavelengths of the broadband filters employed by LSST, alongside determining the mean size and its standard deviation for each respective band. The boundaries of the filters are available in Table 2.1 of the [LSST Science Book](https://ui.adsabs.harvard.edu/abs/2009arXiv0912.0201L/abstract)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7371d2cc-a861-43e3-89a1-2d986ba34f3f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wavelength = 0.5*np.array([320+400, 400+552, 552+691, 691+818])\n",
    "\n",
    "mean_size = [np.mean(objects['size_u_PSF']),\n",
    "             np.mean(objects['size_g_PSF']),\n",
    "             np.mean(objects['size_r_PSF']),\n",
    "             np.mean(objects['size_i_PSF'])]\n",
    "\n",
    "std_size = [np.std(objects['size_u_PSF']),\n",
    "            np.std(objects['size_g_PSF']),\n",
    "            np.std(objects['size_r_PSF']),\n",
    "            np.std(objects['size_i_PSF'])]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f1dc1a3-11f5-4eb0-993e-b27e878c3b06",
   "metadata": {},
   "source": [
    "Fit now a power-law to the mean PSF size in relation to wavelength by employing [`scipy.optimize.curvefit`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html). \n",
    "\n",
    "This function takes the functional form to be fitted as its first argument, followed by two arrays containing the data. The output is a tuple consisting of the optimized parameters that minimize the squared residuals of the function with respect to the measured data, labeled as `popt`. Additionally, it provides the estimated 2D covariance matrix of the fitted parameters, denoted as `pcov`, wherein the diagonal elements offer an estimation of their respective one-standard-deviation errors (for more comprehensive information, refer to [the documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce4e6234-d4da-48e0-874a-b0c64f793bce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "popt, pcov = curve_fit(powerLaw, wavelength, mean_size)\n",
    "perr = np.sqrt(np.diag(pcov))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9afcec9-316a-4299-82f0-248ace3b9df4",
   "metadata": {},
   "source": [
    "Use the fitted parameters to define fitted values, and print the parameters of the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "450779bf-77bf-42c5-9805-297985790b4f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "lamb_grid = np.arange(300, 800, 1)\n",
    "fit_size = powerLaw(lamb_grid, *popt)\n",
    "amplitude, exponent = popt\n",
    "amplitude_error, exponent_error = perr\n",
    "print(f\"Amplitude: {amplitude:.3} +/- {amplitude_error:.3}\")\n",
    "print(f\"Exponent: {exponent:.3} +/- {exponent_error:.3}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea4b9ee6-ee39-4ff3-8ae7-f8675fc129cb",
   "metadata": {},
   "source": [
    "The fitted exponent, `-0.196 +/- 0.0235` is within one-sigma of the expected value of `-1/5=-0.2`. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5841e9fa-406c-4b17-bce9-1f9d2d9025f3",
   "metadata": {},
   "source": [
    "Plot the mean PSF size per band as a function of wavelength."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bb842f4-1a74-4d21-bb83-06236a0ee812",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6, 4))\n",
    "plt.errorbar(wavelength, mean_size, std_size/np.sqrt(len(objects)),\n",
    "             fmt='o--', label='data')\n",
    "plt.plot(lamb_grid, fit_size,\n",
    "         label=r'$\\theta = %.3f \\lambda^{-%.3f}$' % (amplitude, exponent))\n",
    "plt.title(\"Seeing dependence on wavelength.\")\n",
    "plt.xlabel(r'$\\lambda$ (nm)')\n",
    "plt.ylabel(r'$\\sqrt{(Ixx+Iyy)/2}$ (pixel)')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea79d584-117b-4ae6-ae3f-0526053512d0",
   "metadata": {},
   "source": [
    "> **Figure 3**: *The plot above depicts the mean PSF size (expressed as the trace radius in pixels) in relation to the central wavelength for the previously specified LSST filters. Additionally, a power-law fit is presented on the same plot, with the associated amplitude and fitting parameters indicated within the legend.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2906b27f-96fc-4b78-bbee-9f413856960e",
   "metadata": {},
   "source": [
    "Superimpose the size hisograms per band and the plot in Figure 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6da64b0-7505-43d6-adcb-12ff8d1a393b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(figsize=(6, 6))\n",
    "\n",
    "color = 'tab:red'\n",
    "ax1.set_xlabel(r'$\\lambda$ (nm)')\n",
    "ax1.set_ylabel(r'$\\sqrt{(Ixx+Iyy)/2}$ (pixel)')\n",
    "ax1.errorbar(wavelength, mean_size, std_size/np.sqrt(len(objects)),\n",
    "             fmt='o--', label='data')\n",
    "ax1.plot(lamb_grid, fit_size,\n",
    "         label=r'$\\theta = %.3f \\lambda^{-%.3f}$' % (amplitude, exponent))\n",
    "ax1.legend(loc='lower right')\n",
    "ax1.tick_params(axis='y')\n",
    "\n",
    "ax2 = ax1.twiny()\n",
    "\n",
    "plot_filter_colors = {'u': '#56b4e9', 'g': '#008060', 'r': '#ff4000',\n",
    "                      'i': '#850000', 'z': '#6600cc', 'y': '#000000'}\n",
    "\n",
    "ax2.set_xlabel('counts')\n",
    "plt.hist(objects['size_u_PSF'], bins=100, range=[1.6, 2.2],\n",
    "         color=plot_filter_colors['u'], alpha=0.3, label='u',\n",
    "         orientation=\"horizontal\")\n",
    "plt.hist(objects['size_g_PSF'], bins=100, range=[1.6, 2.2],\n",
    "         color=plot_filter_colors['g'], alpha=0.3, label='g',\n",
    "         orientation=\"horizontal\")\n",
    "plt.hist(objects['size_r_PSF'], bins=100, range=[1.6, 2.2],\n",
    "         color=plot_filter_colors['r'], alpha=0.3, label='r',\n",
    "         orientation=\"horizontal\")\n",
    "plt.hist(objects['size_i_PSF'], bins=100, range=[1.6, 2.2],\n",
    "         color=plot_filter_colors['i'], alpha=0.3, label='i',\n",
    "         orientation=\"horizontal\")\n",
    "plt.legend()\n",
    "ax2.tick_params(axis='y')\n",
    "\n",
    "fig.suptitle(\"PSF size histogram per band superimposed with \"\n",
    "             \"mean PSF size as a funtion of wavelength\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0aa29a7-6ddd-497c-b70f-5b147f836ec6",
   "metadata": {},
   "source": [
    "> **Figure 4**: *The figure above superimposes the size distributions per band onto the previously generated plot in Figure 3. The figure includes histograms of the size distributions corresponding to each band (`u`, `g`, `r`, and `i`), rotated 270 degrees to align their bases with the left vertical axis. The upper horizontal axis displays the histogram counts, spanning from 0 to approximately 1200 counts. Notably, this plot serves a dual role, as it incorporates both the data for the PSF size as a function of wavelength and the power-law fit from the preceding plot.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5fecd7f-018f-48df-ba9d-45e96a4b4801",
   "metadata": {},
   "source": [
    "## 5. Rho statistics with the `analysis.tools` package\n",
    "\n",
    "Utilize [`lsst.analysis.tools`](https://github.com/lsst/analysis_tools) to compute the `rho` statistics, which are a collection of two-point correlation functions used to assess the spatial relationships between errors in PSF models. The presence of non-zero values in these statistics signifies the existence of systematic errors in the weak lensing shear correlation function. To understand their definitions, refer to [the analysis tools documentation](https://pipelines.lsst.io/v/daily/py-api/lsst.analysis.tools.actions.vector.CalcRhoStatistics.html#lsst.analysis.tools.actions.vector.CalcRhoStatistics). \n",
    "\n",
    "`lsst.analysis.tools` is a package designed to generate quality assurance (QA) plots and metrics from the outputs of the LSST Science Pipelines. It employs `treecorr` to compute the `rho` statistics. For an introductory demonstration of `analysis_tools`, you can refer to [this Jupyter Notebook created for a session during Rubin's 2022 Project and Community Workshop](https://github.com/lsst-dm/analysis_tools_examples/blob/main/analysis_tools_demo_pcw2022.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d639a047-d86b-4540-b897-54564189bf3f",
   "metadata": {},
   "source": [
    "Create the `RhoStatistics` analysis structure using the `AnalysisTool`. As outlined in the [2022 Rubin PCW `analysis_tools` demo notebook](https://github.com/lsst-dm/analysis_tools_examples/blob/main/analysis_tools_demo_pcw2022.ipynb), the `AnalysisTool` is one of the three analysis structures provided by `analysis_tools`. It encompasses several `AnalysisAction` components, which collectively conduct specific types of analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fa13319-c6b5-4c42-a6f4-21eacfb7ad23",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "atool = RhoStatistics()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62fd3266-2af3-4712-89f8-17795abf216d",
   "metadata": {},
   "source": [
    "There are three types of `AnalysisAction`s or stages: \n",
    "\n",
    "- `prep` (for initial filtering of data), \n",
    "- `process` (for modyfing the data), and \n",
    "- `produce` (for generating final plots of metrics). \n",
    "\n",
    "Adjust a few configuration parameters from the `prep` and `process` stages."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb6ad2b6-b872-45ff-9bba-1c06ed30d59f",
   "metadata": {},
   "source": [
    "Select stars with with signal-to-noise ratio (SNR) > 20."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40b60d5c-0529-4e76-bf0e-4be4b81ac22a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "atool.prep.selectors.snSelector.threshold = 20.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a295e266-f0a9-4c7b-ab3a-23b144e8cfc0",
   "metadata": {},
   "source": [
    "As in a previous section, we can adjust the `treecorr` parameters before calculations, but this time it will be done via `lsst.analysis.tools`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b47bb00f-b4e5-4525-a259-1b1831440f30",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "atool.process.calculateActions.rho.treecorr.nbins = 21\n",
    "atool.process.calculateActions.rho.treecorr.min_sep = 0.1\n",
    "atool.process.calculateActions.rho.treecorr.max_sep = 100.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36b4b77d-4ed0-4354-8057-6ec0867df185",
   "metadata": {},
   "source": [
    "Call `atool.finalize()` after the `atool.prep` and `atool.process` stages are configured."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00ff3bbe-bfde-4756-b149-1e63c3c0e247",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "atool.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcd8f06d-be02-45c1-ba24-76ddf93de7a0",
   "metadata": {},
   "source": [
    "The TAP query in section 3.1 included fields that not used in that section, but that will be used in this section. \n",
    "\n",
    "In order to know which entries will be needed by `analysys_tools` to calculate the `rho` statistics, look at the input schema of the analysis tool `RhoStatistics`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ce8b391-dbf1-451f-b6da-e5c25adfacd5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "input_schema = atool.getInputSchema()\n",
    "needed_catalog_fields = [name[0] for name in list(atool.getInputSchema())]\n",
    "print(needed_catalog_fields)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19ea9678-00bd-4b5c-8b5e-3f7394dd0670",
   "metadata": {},
   "source": [
    "Compute the `rho` statistics in the particular band that we have chosen previously (`i`)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d711b878-1e1c-427e-8405-34f56c663e5b",
   "metadata": {},
   "source": [
    "Start with the `prep` stage in the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3876cc98-a846-4bb7-9adb-0f1aaff0e9d0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "prepResults = atool.prep(objects, band=f\"{band}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6201731f-83eb-4205-8765-a832f5c57a65",
   "metadata": {},
   "source": [
    "Execute the `process` stage, utilizing the outcomes obtained from the `prep` stage. Ignore the `INFO` lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d92dd649-59f2-40b4-861c-4b94b9b04be0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "processResults = atool.process(prepResults, band=f\"{band}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ebdc877-1f65-4c90-b8d7-29aa052c60b5",
   "metadata": {},
   "source": [
    "Access the results from `treecorr` for each `rho` statistic in order to examine the values or perform any other operation on them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ca7d5d0-43fe-48e2-a56a-ec844a7a82d6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(\"Mean angular separation:\\n\", processResults['rho1'].meanr, \"\\n\")\n",
    "print(\"Correlation function:\\n\", processResults['rho1'].xip, \"\\n\")\n",
    "print(\"Error in the correlation function:\\n\", processResults['rho1'].varxip)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f52f3aae-5c3d-4d28-89a7-088de01b97e4",
   "metadata": {},
   "source": [
    "Plot the `rho` statistics using `analysis_tools` "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29d0cbac-2701-4604-8b07-586dc7e1455a",
   "metadata": {},
   "source": [
    "By default, the rho statistics are plotted on a symmetric log scale.\n",
    "The threshold where the switch happens from linear to log scale can be set as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "400821f2-f4ce-4012-8886-071a5cdd7644",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "atool.produce.plot.rhoPlots[\"rho3\"].yLinThresh = 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed77ff9a-ecbb-46ad-b392-2f56bd0df271",
   "metadata": {},
   "source": [
    "At the moment, `atool.produce` does nto provide a direct handle to the plotting keyword arguments. \n",
    "\n",
    "Overwrite the global plotting parameters in order to make figures with a size suitable for a Jupyter notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cee4fab9-399e-4e24-8245-24597ec5c2de",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (2, 2)\n",
    "plt.rcParams['lines.markersize'] = 2\n",
    "plt.rcParams['xtick.labelsize']=5\n",
    "plt.rcParams['ytick.labelsize']=5\n",
    "plt.rcParams['font.size']=5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24ea3e87-d27e-4a2f-8717-33a5d961a41b",
   "metadata": {},
   "source": [
    "Use `analysis_tools` to produce the `rho` statistics plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee19ba4b-a4a0-4da8-b986-e496e978fd5f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "produceResults = atool.produce(processResults, band=f\"{band}\", skymap=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3453cd0-de55-4469-add6-344e82922f54",
   "metadata": {},
   "source": [
    "> **Figure 5** *This figure displays six plots of the `rho` statistics using `analysis_tools`. \n",
    "For other possible configuration options, you can refer to \n",
    "[the documentation](https://pipelines.lsst.io/v/daily/py-api/lsst.analysis.tools.actions.plot.XYPlot.html#lsst.analysis.tools.actions.plot.XYPlot). \n",
    "The vertical axis of each plot displays the specific `rho` statistic, and the horizontal axis represents the separation in arcminutes, ranging from approximately 0.1 to 100 arcminutes on a logarithmic scale.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2229b46-2a92-4596-8021-837afb6f844d",
   "metadata": {},
   "source": [
    "As explained in [Jarvis et al. 2021](https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.1282J/abstract), `rho1` (ρ1), `rho3` (ρ3), and `rho4` (ρ4) represent direct systematic errors in `xi+` (ξ+) with leading coefficients of order unity. `rho2` (ρ2) and `rho5` (ρ5) describe the extent of \"PSF leakage\" that occurs during the shear measurement process.\n",
    "\n",
    "The DC2 area is relatively small and does not exhibit significant PSF variations. The `rho` statistics obtained in the previous examples (for the `i` band, considering that `rho` statistics can vary with different bands) are relatively small and well-behaved. For comparison, you can refer to Figure 26 of the [HSC Y3 shape catalog](https://ui.adsabs.harvard.edu/abs/2022PASJ...74..421L/abstract) and Figure 12 of the [DES Y3 shape catalog analyses](https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.1282J/abstract), which display `rho` statistics obtained from on-sky data from their respective galaxy surveys (HSC Survey and DES). These figures also include the requirements that each `rho` statistic should satisfy (e.g., Eq. 37 of [Li et al. 2022](https://ui.adsabs.harvard.edu/abs/2022PASJ...74..421L/abstract)) to ensure that there are no significant PSF modeling systematic errors that would impact the scientific goals related to cosmic shear."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22b3ff6f-07cb-41f5-b884-f4ca9d623c24",
   "metadata": {},
   "source": [
    "## 6. Exercises for the learner\n",
    "\n",
    "1. Calculate the PSF ellipticity using the definition based on the second moments, as described in section 7 of [Jarvis et al. 2021](https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.1282J/abstract). Create plots similar to those in section 3 to visualize the PSF ellipticity variation across the image.\n",
    "\n",
    "2. Modify the query in section 2.1 to retrieve both point sources and extended sources. Create a size-magnitude diagram using the definition of size from section 2.1 (trace radius). Identify the region in the diagram where the size is constant, representing likely stars during the star-galaxy separation process. Compare your results with Figure 3 in [Jarvis et al. 2021](https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.1282J/abstract).\n",
    "\n",
    "3. Modify the query in section 2.1 to retrieve the required catalog fields according to the input schema for the `RhoStatistics` of `analysis_tools`. Calculate the `rho` statistics for a different band, such as the `g` band. Note that chromatic effects, including Differential Chromatic Refraction (DCR) and wavelength-dependent seeing, have a larger impact in the bluer bands. Use the `rho` statistics as a diagnostic tool to assess whether observations in these bluer bands should be included in weak lensing analyses (refer to the discussion in section 7.4 of [Jarvis et al. 2021](https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.1282J/abstract))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5c666e4-078c-4a71-8d74-7e5652e460e9",
   "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
}