DiaObjectDuplicationIssue_LittleDemo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1bd650c0-fb61-4c2d-a127-f5927a3af6e4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-12-12T19:00:22.188286Z",
     "iopub.status.busy": "2023-12-12T19:00:22.187984Z",
     "iopub.status.idle": "2023-12-12T19:00:22.192878Z",
     "shell.execute_reply": "2023-12-12T19:00:22.192124Z",
     "shell.execute_reply.started": "2023-12-12T19:00:22.188261Z"
    }
   },
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<br>\n",
    "<b> (Little demo) Exploring the DiaObject Duplication Issues in DP0.2</b> <br>\n",
    "Contact author(s): Ryan Lau <br>\n",
    "Last verified to run: 2024-04-16 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_04 <br>\n",
    "Container Size: medium"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bac73400-1c64-43e7-a116-3d92cdce5af2",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "The purpose of this notebook is to demonstrate and characterize the DiaObject duplication issue identified in notebook DP02_07b \"Variable Stars in DP0.2,\" where the same variable star can be identified with multiple DiaOjbectIds. The key goal of this notebook is to inform how this issue may or may not affect your scientific analysis of transients and/or variables. \n",
    "\n",
    "In Section 2, we demonstrate the diaObject duplication issue using the known RR Lyrae variable presented as an example in notebook DP02_07b. This known RR Lyrae star has two diaObjectIds associated with it within a 0.5'' radius, which highlights the issue with the source association algorithm in the LSST Pipeline that should associate diaSources with diaObjects within a 0.5'' radius. \n",
    "\n",
    "Lastly, in Section 3, we conduct a broader investigation of the diaObject duplication issue. We estimate the occurrence rate of duplicate diaObjectIds using a random sample of diaObjects from the diaObject catalog. We also investigate the types of sources affected by this using the TruthSummary table. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43eef134-2a78-49d1-a361-f2ad64c8f368",
   "metadata": {},
   "source": [
    "### 1.1 Package Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90bfc9f1-27c9-404c-9dcd-8a99a2013b6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from lsst.rsp import get_tap_service\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33359ddc-28f7-48c0-bb59-6bcea675caa3",
   "metadata": {},
   "source": [
    "### 1.2 Define Functions and Parameters\n",
    "\n",
    "Setting the plot format parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dffe42a-90fd-420d-8ad5-95ee359bc3d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.style.use('tableau-colorblind10')\n",
    "params = {'axes.labelsize': 24,\n",
    "          'font.size': 20,\n",
    "          'legend.fontsize': 14,\n",
    "          'xtick.major.width': 3,\n",
    "          'xtick.minor.width': 2,\n",
    "          'xtick.major.size': 12,\n",
    "          'xtick.minor.size': 6,\n",
    "          'xtick.direction': 'in',\n",
    "          'xtick.top': True,\n",
    "          'lines.linewidth': 3,\n",
    "          'axes.linewidth': 3,\n",
    "          'axes.labelweight': 3,\n",
    "          'axes.titleweight': 3,\n",
    "          'ytick.major.width': 3,\n",
    "          'ytick.minor.width': 2,\n",
    "          'ytick.major.size': 12,\n",
    "          'ytick.minor.size': 6,\n",
    "          'ytick.direction': 'in',\n",
    "          'ytick.right': True,\n",
    "          'figure.figsize': [10, 8],\n",
    "          'figure.facecolor': 'White'}\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a29e3831-5216-426b-be23-6c53c8d93fcb",
   "metadata": {},
   "source": [
    "Start the TAP service, which we will use for all data retrieval in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96c84ce5-2e30-478e-8981-a0d0017df2f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "976590c4-4738-4446-8885-5d939317860f",
   "metadata": {},
   "source": [
    "## 2. Demonstration of DiaObject Duplication Issue on a Known RR Lyrae Variable\n",
    "\n",
    "In this section, we will present an example of the DiaObject duplication issue using the known RR Lyrae Variable that was also presented as an example variable in notebook DP02_07b. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb7638a5-0907-48d6-8755-ddbd5a8a7d95",
   "metadata": {},
   "source": [
    "### 2.1 Identifying Two DiaObjectIds associated with known variable from DiaSources\n",
    "\n",
    "As in DP02_07b, the known RR Lyrae star we will use is at position (RA, Dec) = (62.1479031, -35.799138) degrees. \n",
    "\n",
    "We define the coordinate (ra and dec) of the known position this variable, and then obtain the diaSource properties including the associated diaObjectIds, total flux, and detector coordinate positions by conducting the following search. Note that we use a narrow search radius of 0.5'', which is the radius used for associating a diaSource with a known diaObject. Lastly, we sort the diaSources by 'midPointTai' for convenience of investigating when the diaObject duplication occurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "339190de-9341-456a-b1bc-434d224e5a23",
   "metadata": {},
   "outputs": [],
   "source": [
    "ra_known_rrl = 62.1479031\n",
    "dec_known_rrl = -35.799138\n",
    "\n",
    "results = service.search(\"SELECT ra, decl, diaObjectId, diaSourceId, ccdVisitId,\"\n",
    "                             \"filterName, midPointTai, psFlux, totFlux, totFluxErr, \"\n",
    "                             \"apFlux, apFluxErr,psFluxErr, snr, x, y \"\n",
    "                             \"FROM dp02_dc2_catalogs.DiaSource \"\n",
    "                             \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\n",
    "                             \"CIRCLE('ICRS',\" + str(ra_known_rrl) + \", \"\n",
    "                             + str(dec_known_rrl) + \", 0.000139)) = 1 \", maxrec=100000)\n",
    "DiaSrcs = results.to_table()\n",
    "DiaSrcs.sort('midPointTai')\n",
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec0ebc67-7430-42a7-a145-ca4e9fd64669",
   "metadata": {},
   "source": [
    "We can identify the unique diaObjects from the resulting 'DiaSrcs' table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e023aa6-0195-40a6-83fb-ce68b0a28d96",
   "metadata": {},
   "outputs": [],
   "source": [
    "DiaObjIds = list(set(DiaSrcs['diaObjectId']))\n",
    "print(DiaObjIds)\n",
    "print('The first DiaObject assigned to this target is %s' %DiaSrcs['diaObjectId'][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28eac65c-9622-47f1-a039-d3806cb381b8",
   "metadata": {},
   "source": [
    "If the DIA source association is working propertly, there should only be one DiaObject associated with the 0.5'' search of a true astrophysical source. The occurance of a duplicate DiaObjectId indicates an issue with the source association. Note that the first DiaObjectId created for this target is 1651589610221862935. The duplicate DiaObjectId is 1651589610221864014.\n",
    "\n",
    "The following cell will identify the MJD when the duplicate DiaObject was created and show where in the DiaSrcs table this occurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efdab588-328c-4b90-8e8e-99a500e14e39",
   "metadata": {},
   "outputs": [],
   "source": [
    "mjd_dup = DiaSrcs[np.where(DiaSrcs['diaObjectId']==1651589610221864014)]['midPointTai'][0]\n",
    "DiaSrcs[np.where(DiaSrcs['midPointTai']==mjd_dup)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc97a411-0ce2-4c47-b809-9b86fb691aa3",
   "metadata": {},
   "source": [
    "Surprisingly, the DiaObject duplication seems to have occured on the same visit where a diaSource was associated with the original diaObject (1651589610221862935) with nearly identical measure DiaSource properties (excepted 'totFlux')."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7f5009c-b9b3-4c31-906a-e6c6b7d32d9c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-15T20:48:53.948198Z",
     "iopub.status.busy": "2024-02-15T20:48:53.947359Z",
     "iopub.status.idle": "2024-02-15T20:48:53.950941Z",
     "shell.execute_reply": "2024-02-15T20:48:53.950363Z",
     "shell.execute_reply.started": "2024-02-15T20:48:53.948171Z"
    }
   },
   "source": [
    "### 2.2 Light curve and coordinate properties of DiaSources associated with known variable\n",
    "\n",
    "We can plot the light curve (y-band) of both diaObjects to investigate if they perhaps exhibit distinct brightnesses and/or levels of variability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffa25207-b2d5-4560-8aac-67e259557533",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(12, 8))\n",
    "DiaObjIds = list(set(DiaSrcs['diaObjectId']))\n",
    "filt = 'y'\n",
    "\n",
    "\n",
    "for j in np.arange(len(DiaObjIds)):\n",
    "    fx = np.where((DiaSrcs['diaObjectId']==DiaObjIds[j]) & (DiaSrcs['filterName']==filt))[0]\n",
    "    ax.plot(DiaSrcs['midPointTai'][fx], DiaSrcs['totFlux'][fx],['p','v','^','o'][j], ms=10,color = ['r','g','b','y'][j], mew=2, mec=['r','g','b','y'][j],\n",
    "            alpha=0.3, label=DiaObjIds[j])\n",
    "\n",
    "#Plotting line indicating MJD where duplicate DiaObject is created\n",
    "plt.axvline(x = mjd_dup,color = 'black', ls = '--', lw = 1)\n",
    "\n",
    "ax.set_xlabel('Modified Julian Date')\n",
    "ax.set_ylabel('TotFlux')\n",
    "ax.set_title('Forced PSF flux measured on the direct image')\n",
    "ax.legend(loc='lower right')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eda32bf5-5947-4676-aaa8-b715be7ef94e",
   "metadata": {},
   "source": [
    "Other than the initial low value of the first diaSource associated with the duplicate diaObject (1651589610221864014), the light curve appears to indicate that both diaObjects are indeed associated with the same variable star.\n",
    "\n",
    "In the following cell, we plot the measured ra and dec offsets from the known coordinates of the variable star for both diaObjects (in all filters)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67f3c0b8-37dc-4b2f-b55c-e85a6a1c0a16",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(10, 10))\n",
    "DiaObjIds = list(set(DiaSrcs['diaObjectId']))\n",
    "\n",
    "\n",
    "for j in np.arange(len(DiaObjIds)):\n",
    "    fx = np.where(DiaSrcs['diaObjectId']==DiaObjIds[j])\n",
    "    ax.plot((ra_known_rrl-DiaSrcs['ra'][fx])*3600, (dec_known_rrl-DiaSrcs['decl'][fx])*3600,\n",
    "               ['p','v','^','o'][j], ms=10, mew=2, mec=['r','g','b','y'][j],\n",
    "               alpha=0.5, color='none', label='%s' % DiaObjIds[j])\n",
    "\n",
    "ax.set_xlabel('ra offset (arcsec)')\n",
    "ax.set_ylabel('dec offset (arcsec)')\n",
    "ax.legend(loc='lower left')\n",
    "\n",
    "\n",
    "#plt.subplots_adjust(wspace=.0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cca7076-f12a-4580-bc14-09105515f180",
   "metadata": {},
   "source": [
    "As shown in the offset coordinate plot, the diaSources associated with both diaObjects are within 0.5'' of eachother and therefore should not have been split into two diaObjects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfd24b22-72ca-4898-813f-f6a9eb60280b",
   "metadata": {},
   "source": [
    "## 3. Broader Investigation of DiaObject Duplication Issue\n",
    "\n",
    "In this section, we estimate the occurrence rate of the diaObject duplication issue by obtaining a random sample of diaObjects from the DiaObject catalog, and verify the types of objects that are affected by this using the TruthSummary table.\n",
    "\n",
    "### 3.1. Estimated Occurrence Rate of DiaObject Duplication Issue\n",
    "\n",
    "First, we grab a random sample of 100 DiaObjects that have at least 10 detections to help ensure they are associated with real astrophysical transients/variables and not arising from spurious diaSource detections. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1da0a1d4-102f-4256-956e-035ced51abc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "nDiaSources_min = 10\n",
    "\n",
    "results = service.search(\"SELECT TOP 100 \"\n",
    "                         \"ra, decl, diaObjectId, nDiaSources \"\n",
    "                         \"FROM dp02_dc2_catalogs.DiaObject \"\n",
    "                         \"WHERE nDiaSources > \"+str(nDiaSources_min)+\" \")\n",
    "DiaObjs = results.to_table()\n",
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1356de33-4762-4486-a5eb-b76a22b88c6f",
   "metadata": {},
   "source": [
    "In the following cell, we conduct a coordinate search through the DiaSource catalog for the sample of diaObjects obtained above to identify DiaSources within a 0.5'' radius of the diaObject coordinates. \n",
    "\n",
    "The resulting DiaSources will notably include the diaObjectId(s) they are associated with. The diaObject duplication issue will manifest as multiple diaObjectIds for a specific coordinate position.\n",
    "\n",
    "We will keep track of the number of duplicate diaObjectIds for each of the diaObject position in an array labeled \"NDup,\" which we use to create an \"nDuplicates\" column for the DiaObjs table. The following cell should take less than 30 seconds to execute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a94dcc2-403f-4df8-a2be-e67bd675b3cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "NDup = np.zeros(len(DiaObjs))\n",
    "for i in np.arange(len(DiaObjs)):\n",
    "    ra = DiaObjs['ra'][i]\n",
    "    decl = DiaObjs['decl'][i]\n",
    "    results = service.search(\"SELECT ra, decl, diaObjectId, diaSourceId, ccdVisitId,\"\n",
    "                             \"filterName, midPointTai \"\n",
    "                             \"FROM dp02_dc2_catalogs.DiaSource \"\n",
    "                             \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\n",
    "                             \"CIRCLE('ICRS',\" + str(ra) + \", \"\n",
    "                             + str(decl) + \", 0.000139)) = 1 \", maxrec=100000) #0.5'' radius coordinate search\n",
    "    DiaSrcs = results.to_table()\n",
    "    del results\n",
    "    NDup[i]=len(list(set(DiaSrcs['diaObjectId'])))-1\n",
    "DiaObjs['nDuplicates']=NDup"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb86f731-31bf-40db-8007-f6e15cca3c31",
   "metadata": {},
   "source": [
    "We can estimate the occurrence frequency of the diaObject duplication issue by calculating what fraction of the DiaObject sample has at least one duplicate diaObject."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d5af851-3af0-45ef-abd5-a7aa8f7de1c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(NDup[NDup>0])/len(NDup)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "472bae45-7b31-4b3f-a038-debf331bb7be",
   "metadata": {},
   "source": [
    "The results of the above cell indicates that >~50% of diaObjects are affected by the diaObject duplication issue, which could present significant issues for statistical analyses of transients and variables. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5533d89a-4c62-49ad-9b94-d78091d95157",
   "metadata": {},
   "source": [
    "### 3.2. Investigating DiaObject Truth Type affected by Duplication Issue\n",
    "\n",
    "In order to investigate the type of sources (i.e. variables vs transients) that are producing the duplication issue, we will grab a random sample of diaObjects associated with variable stars and SNe and identify how many diaObjects of each type have duplicates. In order to verify the diaObject type, we utilize the TruthSummary table, the DESC DC2 truth catalog as described in arXiv:2101.04855.\n",
    "\n",
    "**First, we will explore the variable stars with diaObject duplicates**. We conduct a (1''-radius) coordinate search in the TruthSummary table for the diaObjects with duplicates identified in the random sample of diaObjects obtained in Sec. 3.1. The following cell should take around 30 seconds to execute.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d41fb899-fe49-49cf-b336-6a5a979baf5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "Var = np.zeros(len(DiaObjs))\n",
    "for i in np.arange(len(DiaObjs)):\n",
    "    ra = DiaObjs['ra'][i]\n",
    "    decl = DiaObjs['decl'][i]\n",
    "    results = service.search(\"SELECT ts.ra, ts.dec, ts.is_variable, ts.truth_type \"\n",
    "                             \"FROM dp02_dc2_catalogs.TruthSummary AS ts \"\n",
    "                             \"WHERE CONTAINS(POINT('ICRS', ts.ra, ts.dec), \"\n",
    "                             \"CIRCLE('ICRS',\" + str(ra) + \", \"\n",
    "                             + str(decl) + \", 0.00028)) = 1 \", maxrec=100000)\n",
    "    SrcTruth = results.to_table()\n",
    "    if 1 in SrcTruth['is_variable'] and 2 in SrcTruth[SrcTruth['is_variable']==1]['truth_type']:\n",
    "        Var[i] = 1\n",
    "    else:\n",
    "        Var[i] = 0\n",
    "    del results\n",
    "DiaObjs['Var']=Var\n",
    "\n",
    "print('Out of the %s diaObjects that are associated with a variable star'\n",
    "      ' (i.e. is_variable = 1 and truth_type = 2), %s have duplicates.' % \n",
    "      (len(DiaObjs[DiaObjs['Var']==1]),len(DiaObjs[(DiaObjs['nDuplicates']>0) & (DiaObjs['Var']==1)])))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40cff8ef-761a-453b-adf6-279aa64b08fc",
   "metadata": {},
   "source": [
    "The results of the search in the previous cell indicate that roughly ~50% of the diaObjects associated with a variable star has duplicate diaObjects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d30e745-fe6b-4010-a91e-5703ce988c9b",
   "metadata": {},
   "source": [
    "**Now we explore the SNe with diaObject duplicates.** In order to grab a sample of Type Ia SN candidates, we utilize the search parameters described in the DP02 07a notebook \"DiaObject Sample Identification.\" (See notebook 07a for more details)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1052ed80-7efb-4b11-9ca8-542c840a2a04",
   "metadata": {},
   "outputs": [],
   "source": [
    "snia_peak_mr_min = 18.82\n",
    "snia_peak_mr_max = 22.46\n",
    "\n",
    "snia_ampl_mr_min = 1.5\n",
    "snia_ampl_mr_max = 5.5\n",
    "\n",
    "snia_peak_mg_max = 24.0\n",
    "snia_peak_mi_max = 24.0\n",
    "\n",
    "snia_nDiaSources_min = 15\n",
    "snia_nDiaSources_max = 100\n",
    "\n",
    "snia_duration_min = 50\n",
    "snia_duration_max = 300"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa22e735-461c-4569-b084-3bf83c5b1a21",
   "metadata": {},
   "source": [
    "Conducting a search for 150 diaObjects associated with Type Ia SN candidates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bbd3342-1516-4810-a8f8-9e66b6f41489",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT TOP 150 \"\n",
    "                         \"ra, decl, diaObjectId, nDiaSources, \"\n",
    "                         \"scisql_nanojanskyToAbMag(rPSFluxMin) AS rMagMax, \"\n",
    "                         \"scisql_nanojanskyToAbMag(rPSFluxMax) AS rMagMin, \"\n",
    "                         \"scisql_nanojanskyToAbMag(gPSFluxMax) AS gMagMin, \"\n",
    "                         \"scisql_nanojanskyToAbMag(iPSFluxMax) AS iMagMin \"\n",
    "                         \"FROM dp02_dc2_catalogs.DiaObject \"\n",
    "                         \"WHERE nDiaSources > \"+str(snia_nDiaSources_min)+\" \"\n",
    "                         \"AND nDiaSources < \"+str(snia_nDiaSources_max)+\" \"\n",
    "                         \"AND scisql_nanojanskyToAbMag(rPSFluxMax) > \"+str(snia_peak_mr_min)+\" \"\n",
    "                         \"AND scisql_nanojanskyToAbMag(rPSFluxMax) < \"+str(snia_peak_mr_max)+\" \"\n",
    "                         \"AND scisql_nanojanskyToAbMag(gPSFluxMax) < \"+str(snia_peak_mg_max)+\" \"\n",
    "                         \"AND scisql_nanojanskyToAbMag(iPSFluxMax) < \"+str(snia_peak_mi_max)+\" \"\n",
    "                         \"AND scisql_nanojanskyToAbMag(rPSFluxMin)\"\n",
    "                         \" - scisql_nanojanskyToAbMag(rPSFluxMax) < \"+str(snia_ampl_mr_max)+\" \"\n",
    "                         \"AND scisql_nanojanskyToAbMag(rPSFluxMin)\"\n",
    "                         \" - scisql_nanojanskyToAbMag(rPSFluxMax) > \"+str(snia_ampl_mr_min)+\" \")\n",
    "\n",
    "DiaObjsSN = results.to_table()\n",
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "514c1d90-de05-4bf3-919d-bf51825f4c57",
   "metadata": {},
   "source": [
    "In the following cell, we identify duplicates in the sample of SN candidates. The following cell takes around 30 seconds to execute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82118054-5c0f-48d4-b816-870c050964ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "NDup = np.zeros(len(DiaObjsSN))\n",
    "for i in np.arange(len(DiaObjsSN)):\n",
    "    ra = DiaObjsSN['ra'][i]\n",
    "    decl = DiaObjsSN['decl'][i]\n",
    "    results = service.search(\"SELECT ra, decl, diaObjectId, diaSourceId, ccdVisitId,\"\n",
    "                             \"filterName, midPointTai \"\n",
    "                             \"FROM dp02_dc2_catalogs.DiaSource \"\n",
    "                             \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\n",
    "                             \"CIRCLE('ICRS',\" + str(ra) + \", \"\n",
    "                             + str(decl) + \", 0.000139)) = 1 \", maxrec=100000) #0.5'' radius coordinate search\n",
    "    DiaSrcs = results.to_table()\n",
    "    del results\n",
    "    NDup[i]=len(list(set(DiaSrcs['diaObjectId'])))-1\n",
    "DiaObjsSN['nDuplicates']=NDup"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9079e1b7-8d28-4c1d-8a36-627536e8f523",
   "metadata": {},
   "source": [
    "Laslty, we utilize the TruthSummary table to confirm which candidates are true SNe and determine how many of the diaObjects associated with true SN have duplicates. The following cell takes around 1 min to execute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb14f67e-4083-4fdb-8a08-27466d69430a",
   "metadata": {},
   "outputs": [],
   "source": [
    "SN = np.zeros(len(DiaObjsSN))\n",
    "for i in np.arange(len(DiaObjsSN)):\n",
    "    ra = DiaObjsSN['ra'][i]\n",
    "    decl = DiaObjsSN['decl'][i]\n",
    "    results = service.search(\"SELECT ts.ra, ts.dec, ts.is_variable, ts.truth_type \"\n",
    "                             \"FROM dp02_dc2_catalogs.TruthSummary AS ts \"\n",
    "                             \"WHERE CONTAINS(POINT('ICRS', ts.ra, ts.dec), \"\n",
    "                             \"CIRCLE('ICRS',\" + str(ra) + \", \"\n",
    "                             + str(decl) + \", 0.00028)) = 1 \", maxrec=100000)\n",
    "    SrcTruth = results.to_table()\n",
    "    if 3 in SrcTruth['truth_type']:\n",
    "        SN[i] = 1\n",
    "    else:\n",
    "        SN[i] = 0\n",
    "    del results\n",
    "DiaObjsSN['SN']=SN\n",
    "print('Out of the %s diaObjects that are associated with a supernova (i.e. truth_type = 3), %s have duplicates.' \n",
    "      % (len(DiaObjsSN[DiaObjsSN['SN']==1]),len(DiaObjsSN[(DiaObjsSN['nDuplicates']>0) & (DiaObjsSN['SN']==1)])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aae8c34b-d400-4527-9d7f-c8eb3cdfa741",
   "metadata": {},
   "source": [
    "The results of the search in the previous cell indicate that roughly ~20% of the diaObjects associated with a SN has duplicate diaObjects. \n",
    "\n",
    "**Note.** The results of the diaObject search 3 cells above will result in a difference set of 150 diaObjects each time the search is executed. Therefore, the number of diaObjects associated with SN and those with duplicates in the previous cell will change each time the search is re-run.\n",
    "\n",
    "In conclusion, the analysis indicates that diaObjects associated with variable stars are more likely to be impacted than diaObjects associated with transients / SNe. "
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}