DiaObject Sample Identification

Contact author: Melissa Graham
Last verified to run: 2024-08-19
LSST Science Pipelines version: Weekly 2024_32
Container Size: medium
Targeted learning level: intermediate

Description: To use the DiaObject table parameters to identify a sample of time-variable objects of interest, and to investigate the forced photometry lightcurves of a specific object of interest with the ForcedSourceOnDiaObject table.

Skills: Use the TAP service and the DP0.2 DiaObject and DiaSource tables.

LSST Data Products: TAP tables dp02_dc2_catalogs.DiaObject, DiaSource, and ForcedSourceOnDiaObject.

Packages: lsst.rsp, astropy.cosmology, astropy.stats, numpy, matplotlib

Credit: Originally developed by Melissa Graham and the Rubin Community Science Team for Data Preview 0. Updated for ForcedSourceOnDiaObject by Ryan Lau.

Get Support: Find DP0-related documentation and resources at dp0.lsst.io. Questions are welcome as new topics in the Support - Data Preview 0 Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction

This notebook demonstrates how to use the DiaObject table's lightcurve summary statistics to identify samples of objects for further study, and how to retrieve and plot lightcurves for DiaObjects of interest from the DiaSource table.

As an example, a sample of potential low-redshift Type Ia supernovae are identified and retrieved from the DiaObject table.

Notice: This notebook does not produce a pure or complete sample of SNIa, just a sample of potential well-sampled low-redshift SNIa. This is decent enough to use as a demonstration, or a starting point for a more rigorous classification process. Different science goals will have different requirements on sample identification.

There are three different types of catalogs built from the difference images:

1. DiaSource: astrometric and photometric measurements for sources detected with signal-to-noise ratio > 5 in the difference images.
2. DiaObject: derived summary parameters for DiaSources associated by sky location, including lightcurve statistics.
3. ForcedSourceOnDiaObject: point-source forced-photometry measurements on individual single-epoch visit images and difference images, based on and linked to the entries in the DiaObject table.

Learn more about the contents of the DiaObject, DiaSource, ForcedSourceOnDiaObjects tables in the DP0.2 Documentation.

This notebook uses the Table Access Protocol (TAP) service to access the catalogs, but they are also available via the Butler.

1.1. Package Imports

lsst.rsp: The LSST Science Pipelines package for RSP functionality such as the TAP service (pipelines.lsst.io).

astropy.cosmology: An open-source package of cosmology tools (the astropy cosmology documentation).

astropy.stats: An open-source package of astronomy-related statistical tools, used in this notebook to look for outliers in the forced-source lightcurves.

import numpy
import matplotlib.pyplot as plt
from astropy.cosmology import FlatLambdaCDM
from astropy.stats import sigma_clipped_stats

from lsst.rsp import get_tap_service

1.2. Define functions and parameters

Start the TAP service.

service = get_tap_service("tap")

Set the cosmology to use with the astropy.cosmology package.

cosmo = FlatLambdaCDM(H0=70, Om0=0.3)

Set a few parameters to use later, when plotting lightcurves.

plt.style.use('tableau-colorblind10')

plot_filter_labels = ['u', 'g', 'r', 'i', 'z', 'y']
plot_filter_colors = {'u': '#0c71ff', 'g': '#49be61', 'r': '#c61c00',
                      'i': '#ffc200', 'z': '#f341a2', 'y': '#5d0000'}
plot_filter_symbols = {'u': 'o', 'g': '^', 'r': 'v', 'i': 's', 'z': '*', 'y': 'p'}

2. Understand the DiaObject table's contents

DiaSource Table: Measurements for sources detected with a signal-to-noise ratio > 5 in the difference images (i.e., lightcurves). Note that these measurements are not forced photometry at the locations of all DiaObjects, but detections only.

DiaObject Table: Summary parameters for DiaSources associated by coordinate (including lightcurve summary statistics).

Option to print all of the available column names for the DiaObject table.

# results = service.search("SELECT column_name from TAP_SCHEMA.columns "
#                          "WHERE table_name = 'dp02_dc2_catalogs.DiaObject'")
# for column_name in results['column_name']:
#     print(column_name)

2.1. Review these descriptions for selected lightcurve summary statistics

The DiaObjects table includes the following lightcurve summary statistics, which are derived from the contents of the DiaSource table.

nDiaSources: The number of difference-image detections in any filter (i.e., number of DiaSources associated with a given DiaObject).

The following statistics are all based on difference-image point source (PS) flux values for each filter [f].
[f]PSFluxMin: The faintest flux.
[f]PSFluxMax: The brightest flux.
[f]PSFluxMean: The average flux.
[f]PSFluxSigma: The standard deviation of the fluxes.
[f]PSFluxMAD: The mean absolute deviation of the fluxes (i.e., the average distance from the mean).
[f]PSFluxChi2: The Chi2 statistic for the scatter of the fluxes around the mean.
[f]PSFluxNdata: The number of data points used to compute [f]PSFluxChi2.
[f]PSFluxSkew: A measure of asymmentry in the distribution of fluxes about the mean (where 0 means symmetric).
[f]PSFluxStetsonJ: A variability index developed for Cepheids (defined in Stetson 1996).
[f]PSFluxPercentile05, 25, 50, 75, 95: Derived from the cumulative distribution of flux values.

The following statistics are all based on the direct-image total (TOT) flux values for each filter [f].
[f]TOTFluxMean: The average flux.
[f]TOTFluxSigma: The standard deviation of the fluxes.

Notice: The DP0.2 DiaObject table is missing some variability characterization parameters (DMTN-118) and host association parameters (DMTN-151) which will exist for future data releases.

2.2. Plot summary-statistic histograms for a random sample of DiaObjects

In order to learn a bit more about these lightcurve summary statistics, histograms of their values for a random selection of DiaObjects are plotted.

First, retrieve a random sample of DiaObjects.

query = "SELECT TOP 100000 "\
        "ra, decl, diaObjectId, nDiaSources, "\
        "rPSFluxMin, rPSFluxMax, rPSFluxMean, rPSFluxSigma, "\
        "rPSFluxMAD, rPSFluxChi2, rPSFluxNdata, rPSFluxSkew, "\
        "rPSFluxStetsonJ, rPSFluxPercentile05, rPSFluxPercentile25, "\
        "rPSFluxPercentile50, rPSFluxPercentile75, rPSFluxPercentile95, "\
        "rTOTFluxMean, rTOTFluxSigma "\
        "FROM dp02_dc2_catalogs.DiaObject"
print(query)

SELECT TOP 100000 ra, decl, diaObjectId, nDiaSources, rPSFluxMin, rPSFluxMax, rPSFluxMean, rPSFluxSigma, rPSFluxMAD, rPSFluxChi2, rPSFluxNdata, rPSFluxSkew, rPSFluxStetsonJ, rPSFluxPercentile05, rPSFluxPercentile25, rPSFluxPercentile50, rPSFluxPercentile75, rPSFluxPercentile95, rTOTFluxMean, rTOTFluxSigma FROM dp02_dc2_catalogs.DiaObject

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED

DiaObjs = job.fetch_result().to_table()

Plot the number of DiaSources per DiaObject

Plot the distribution of the number of DiaSources per DiaObject (i.e., the total number of difference-image detections in any filter; at left), and the distribution of the number of r-band DiaSources per DiaObject (at right).

Notice how the distribution is peaked at small numbers of DiaSources (or r-band detections) -- these are time-variable sources that were only detected in a few difference images, and are most likely to have a faint time-variable flux (e.g., high-z supernovae that were only detectable for a short duration).

fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)

ax[0].hist(DiaObjs['nDiaSources'], bins=50, log=True, color='gray')
ax[0].set_xlabel('nDiaSources')
ax[0].set_ylabel('log(Number of DiaObjects)')

ax[1].hist(DiaObjs['rPSFluxNdata'], bins=50, log=True, color=plot_filter_colors['r'])
ax[1].set_xlabel('rPSFluxNdata')

plt.tight_layout()
plt.show()

Figure 1: At left, a histogram (grey) showing the distribution of the number of diaSources per diaObject, across all filters. At right, a histogram (dark red) of the number of r-band PSF flux detections per diaObject.

Plot distributions of min/max PS fluxes

Plot the distribution of minimum and maximum r-band PS flux from the difference-image detections.

Notice how the PS fluxes can be negative because they are measured on the difference images. Any source that is fainter in the direct image than in the template image will be detected as a negative-flux source in the difference image.

Warning: Take care when converting difference-image fluxes into magnitudes, because they can be negative.

Notice how these distributions are peaked at low time-variable flux components -- likely the same types of faint sources which were only detectable for a short duration (as shown above).

fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)

ax[0].hist(DiaObjs['rPSFluxMin'], bins=100, log=True, color=plot_filter_colors['r'])
ax[0].set_xlabel('rPSFluxMin')
ax[0].set_ylabel('log(Number of DiaObjects)')

ax[1].hist(DiaObjs['rPSFluxMax'], bins=100, log=True, color=plot_filter_colors['r'])
ax[1].set_xlabel('rPSFluxMax')

plt.tight_layout()
plt.show()

Figure 2: The distribution of the minimum (left) and maximum (right) r-band PSF flux measurements for diaObjects.

Confirm that DiaObjects with few DiaSources have min/max PS fluxes close to zero

Plot the correlation between the number of r-band DiaSources per DiaObject, and the minimum (left) and maximum (right) difference-image PS flux values.

Notice that it is, indeed, the fainter DiaObjects that have fewer DiaSource detections, as hypothesized above.

fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)

ax[0].plot(numpy.log10(numpy.absolute(DiaObjs['rPSFluxMin'])),
           DiaObjs['rPSFluxNdata'],
           'o', ms=5, alpha=0.2, mew=0, color=plot_filter_colors['r'])
ax[0].set_xlabel('log |rPSFluxMin|')
ax[0].set_ylabel('rPSFluxNdata')

ax[1].plot(numpy.log10(numpy.absolute(DiaObjs['rPSFluxMax'])),
           DiaObjs['rPSFluxNdata'],
           'o', ms=5, alpha=0.2, mew=0, color=plot_filter_colors['r'])
ax[1].set_xlabel('log |rPSFluxMax|')
ax[1].set_ylabel('rPSFluxNdata')

plt.tight_layout()
plt.show()

Figure 3: The r-band PSF flux minimum (left) and maximum (right) for diaObjects, versus the number of r-band PSF flux detections. This plot demonstrates that diaObjects with fewer detections (smaller y-axis values) are fainter (have lower PSF flux values).

Other PS Flux statistics

It is left as an exercise for the learner to investigate typical values of the other PS Flux summary statistics, such as rPSFluxMean, rPSFluxSigma, rPSFluxMAD, rPSFluxChi2, rPSFluxSkew, and/or rPSFluxStetsonJ.

Plot distributions of TOT flux mean and sigma

Plot the distributions of the DiaObject parameters (mean and sigma) derived from the total fluxes from the direct-images.

Note that the TOT fluxes cannot be negative because they are measured on the direct images. The TOT fluxes are safe to convert to magnitudes.

Notice how the distribution is not so sharply peaked at low flux, compared to the PS Flux plots above.

fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)

ax[0].hist(DiaObjs['rTOTFluxMean'], bins=50, log=True, color=plot_filter_colors['r'])
ax[0].set_xlabel('rTOTFluxMean')
ax[0].set_ylabel('log(Number of DiaObjects)')

ax[1].hist(DiaObjs['rTOTFluxSigma'], bins=50, log=True, color=plot_filter_colors['r'])
ax[1].set_xlabel('rTOTFluxSigma')

plt.tight_layout()
plt.show()

Figure 4: The distribution of the mean (left) and standard deviation (right) of the total flux (from the 'direct' or science image) measurements for diaObjects.

2.3. Investigate the DiaSource data for one random, bright, well-sampled DiaObject

Choose a DiaObject that was detected >40 times in an r-band difference image, had an average total (direct-image) flux > 1e6, and an average PS difference-image flux > 5e5.

tx = numpy.where((DiaObjs['rPSFluxNdata'] > 40)
                 & (DiaObjs['rTOTFluxMean'] > 1000000)
                 & (DiaObjs['rPSFluxMean'] > 500000))[0]

use_index = tx[0]
use_diaObjectId = DiaObjs['diaObjectId'][tx[0]]
del tx

Retrieve the DiaSource data for this DiaObject.

query = "SELECT ra, decl, diaObjectId, diaSourceId, ccdVisitId,"\
        "filterName, midPointTai, psFlux, totFlux, totFluxErr, "\
        "apFlux, apFluxErr,psFluxErr, snr  "\
        "FROM dp02_dc2_catalogs.DiaSource "\
        "WHERE diaObjectId = "+str(use_diaObjectId)
print(query)

SELECT ra, decl, diaObjectId, diaSourceId, ccdVisitId,filterName, midPointTai, psFlux, totFlux, totFluxErr, apFlux, apFluxErr,psFluxErr, snr  FROM dp02_dc2_catalogs.DiaSource WHERE diaObjectId = 1249502605990690829

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED

DiaSrcs = job.fetch_result().to_table()

Sort the DiaSrcs table based on the mid-point time of the exposure, and create an array fx that indexes all observations obtained with the r-band filter.

DiaSrcs.sort('midPointTai')
fx = numpy.where(DiaSrcs['filterName'] == 'r')[0]
print('Number of r-band observations: ', len(fx))

Number of r-band observations:  42

Plot the lightcurves for one DiaObject

Plot the difference-image (PSFlux) and direct-image (TOTFlux) lightcurves. Mark the DiaObject minimum and maximum flux with solid lines; average flux with dashed lines; and the standard deviation in flux with dotted lines.

The purpose of this plot is to serve as a visual aid for understanding the lightcurve statistic summary parameters.

fig, ax = plt.subplots(2, figsize=(10, 8), sharey=False, sharex=False)

ax[0].plot(DiaSrcs['midPointTai'][fx], DiaSrcs['psFlux'][fx],
           plot_filter_symbols['r'], ms=15, mew=0, alpha=0.5, color=plot_filter_colors['r'])
ax[0].axhline(DiaObjs['rPSFluxMin'][use_index])
ax[0].axhline(DiaObjs['rPSFluxMax'][use_index])
ax[0].axhline(DiaObjs['rPSFluxMean'][use_index], ls='dashed')
ax[0].axhline(DiaObjs['rPSFluxMean'][use_index] - DiaObjs['rPSFluxSigma'][use_index], ls='dotted')
ax[0].axhline(DiaObjs['rPSFluxMean'][use_index] + DiaObjs['rPSFluxSigma'][use_index], ls='dotted')
ax[0].set_xlabel('Modified Julian Date')
ax[0].set_ylabel('PS Flux')
ax[0].set_title('Difference-Image PSFlux r-band Lightcurve')

ax[1].plot(DiaSrcs['midPointTai'][fx], DiaSrcs['totFlux'][fx],
           plot_filter_symbols['r'], ms=15, mew=0, alpha=0.5, color=plot_filter_colors['r'])
ax[1].axhline(DiaObjs['rTOTFluxMean'][use_index], ls='dashed')
ax[1].axhline(DiaObjs['rTOTFluxMean'][use_index] - DiaObjs['rTOTFluxSigma'][use_index], ls='dotted')
ax[1].axhline(DiaObjs['rTOTFluxMean'][use_index] + DiaObjs['rTOTFluxSigma'][use_index], ls='dotted')
ax[1].set_xlabel('Modified Julian Date')
ax[1].set_ylabel('Total Flux')
ax[1].set_title('Direct-Image TOTFlux r-band Lightcurve')

plt.tight_layout()
plt.show()

Figure 5: The r-band light curve for a random, well-sampled diaObject. The top panel uses the difference-image flux whereas the bottom panel uses the direct-image flux. Horizontal lines mark the minimum and maximum flux (solid), average flux (dashed), and the standard deviation (dotted).

Plot the distribution of flux measurements for one DiaObject

Plot the distribution (left) and normalized cumulative distribution (right) of difference-image fluxes (PSFlux), along with the relevant DiaObject characterization parameters (e.g., mean, sigma, and percentiles). Mark the average flux with dashed lines and the standard deviation in flux with dotted lines.

The purpose of this plot is to serve as a visual aid for understanding the lightcurve statistic summary parameters.

fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)

ax[0].hist(DiaSrcs['psFlux'][fx], bins=20, color=plot_filter_colors['r'],
           label='skew=%5.2f' % DiaObjs['rPSFluxSkew'][use_index])
ax[0].axvline(DiaObjs['rPSFluxMean'][use_index], ls='dashed')
ax[0].axvline(DiaObjs['rPSFluxMean'][use_index] - DiaObjs['rPSFluxSigma'][use_index], ls='dotted')
ax[0].axvline(DiaObjs['rPSFluxMean'][use_index] + DiaObjs['rPSFluxSigma'][use_index], ls='dotted')
ax[0].set_xlabel('r-band PS Flux')
ax[0].set_ylabel('Number of DiaSources')
ax[0].legend(loc='upper left')

ax[1].hist(DiaSrcs['psFlux'][fx], bins=len(fx), color=plot_filter_colors['r'],
           cumulative=True, density=True, histtype='step')
ax[1].plot(DiaObjs['rPSFluxPercentile05'][use_index], 0.05, '*', ms=10, color='black',
           label='percentiles: 0.05, 0.25, 0.50, 0.75, and 0.95')
ax[1].plot(DiaObjs['rPSFluxPercentile25'][use_index], 0.25, '*', ms=10, color='black')
ax[1].plot(DiaObjs['rPSFluxPercentile50'][use_index], 0.50, '*', ms=10, color='black')
ax[1].plot(DiaObjs['rPSFluxPercentile75'][use_index], 0.75, '*', ms=10, color='black')
ax[1].plot(DiaObjs['rPSFluxPercentile95'][use_index], 0.95, '*', ms=10, color='black')
ax[1].set_xlabel('r-band PS Flux')
ax[1].set_ylabel('Cumulative Fraction')
ax[1].legend(loc='upper left')

plt.tight_layout()
plt.show()

Figure 6: A visual demonstration of diaObject lightcurve parameters. At left, the distribution of r-band PSF difference-image flux measurements for the selected diaObject, the skew value of which is noted in the legend. Vertical lines mark the average flux (dashed) and standard deviation (dotted). Note the influence of temporal sampling, as the histogram shows many detections near the maximum PSF flux, yet in Figure 5 it is evident these are temporally clustered observations. At right, the normalized cumulative distribution of flux measurements with percentile values marked as five-point stars.

Clean up by deleting a few arrays that are no longer needed.

del DiaObjs, DiaSrcs, fx
del use_index, use_diaObjectId

3. Identify DiaObjects that are potential Type Ia Supernovae

For this example, a sample of potential low-redshift, well-sampled Type Ia supernovae (SNIa) are identified.

Compared to other types of supernovae, SNIa have homogenous lightcurves with very similar peak absolute brightnesses (about -19 mag in B-band), and similar rise and decline times (i.e., similar observed lightcurve amplitudes and durations, for a survey with a given limiting magnitude).

In LSST-like data sets with an r-band limiting magnitude of about 24.5 mag, such as the DC2 simulation, low-redshift SNIa (0.1 < z < 0.3) lightcurves would look approximately like those plotted in the image below.

Figure 7: Example light curves of SNe Ia at redshift 0.1 (green) and 0.3 (orange), from Nugent et al. (2002; templates).

In this example, the expected peak brightnesses, amplitude, and duration for the lightcurves of low-redshift SNIa are used to construct a TAP query on the DiaObjects table to identify a sample of potential SNIa.

Warning: The following is a very rough, back-of-the-envelope way to identify a bunch of DiaObjects with parameters that suggest they are SNIa-like. The resulting sample is not scientifically useful (not pure, not complete) – it is just a few simple cuts that will get you a bunch of DiaObjects for which we can demonstrate how to plot lightcurves using the DiaSource table. Real SN-finding algorithms are photometric classifiers that fit the lightcurves with templates, rigorously, and provide probabilistic classifications. The text in Section 1 of the notebook describes this.

3.1. Establish lightcurve parameter constraints to identify potential SNIa

Define the desired redshift boundaries, the known approximate peak absolute magnitude for SNeIa (-19 mag), and the desired peak range to use to create the TAP query. This range roughly includes the intrinsic diversity in SNIa brightness and color, and host-galaxy reddening.

redshift_min = 0.1
redshift_max = 0.3

snia_peak_mag = -19.0
snia_peak_mag_range = 0.5

Use the astropy.cosmology package to convert redshift to distance modulus. Define the minimum and maximum peak apparent r-band magnitudes -- allowing for the intrinsic diversity range specified above -- to use in the TAP query.

snia_peak_mr_min = cosmo.distmod(redshift_min).value + snia_peak_mag - snia_peak_mag_range
snia_peak_mr_max = cosmo.distmod(redshift_max).value + snia_peak_mag + snia_peak_mag_range
print('The minimum and maximum apparent r-band magnitudes '
      'to use in the TAP query are %5.2f and %5.2f mag.' %
      (snia_peak_mr_min, snia_peak_mr_max))

The minimum and maximum apparent r-band magnitudes to use in the TAP query are 18.82 and 22.46 mag.

Define maximum magnitudes in the g- and i-bands to use in the TAP query. The point of this is to simply enforce detection in at least the three filters g, r, and i. With knowledge of SNIa colors, this could be made more constraining, but for the purposes of this example these values are fine.

snia_peak_mg_max = 24.0
snia_peak_mi_max = 24.0

Define the r-band minimum and maximum lightcurve amplitudes to use in the TAP query (i.e., the difference between the brightest and faintest detections in the difference image, in magnitudes). Well-sampled, low-redshift SNIa should be observed to change by at least 1.5 mag (z=0.3), and up to 5.5 mag (z=0.1).

snia_ampl_mr_min = 1.5
snia_ampl_mr_max = 5.5

Define the minimum and maximum number of DiaSources (i.e., difference-image detections) to use in the TAP query. The goal of this example is to identify potential well-sampled Type Ia supernovae, and here a minimum of 15 detections is used. Since the DC2 dataset was simulated using a baseline observing strategy (and does not include deep drilling fields), there are no more than 100 visits per year per field. Any DiaObject with >100 DiaSources had a duration >1 year, and is not a SNIa.

snia_nDiaSources_min = 15
snia_nDiaSources_max = 100

Define the minimum and maximum lightcurve duration, in days. The duration is the time between the first and last difference-image detection in any filter. As seen in the template lightcurve plot above, SNIa at redshifts 0.1 < z < 0.3 will have durations of 50 to 300 days.

snia_duration_min = 50
snia_duration_max = 300

Notice: Of the above parameters defined to identify potential SNIa, only the lightcurve duration is not represented in the DiaObjects table, and cannot be included as a constraint in the TAP query used below. Instead, the lightcurve durations are calculated and used in Section 3.3.

3.2. Retrieve a sample of potentially SNIa-like DiaObjects

Only retrieve 1000 DiaObjects for this example. When the TAP query completes, transfer the results to an astropy table.

This TAP query takes about a minute.

Note the use of the function scisql_nanojanskyToAbMag to retrieve magnitudes directly rather than having to convert them from fluxes after the query.

query = "SELECT TOP 1000 "\
        "ra, decl, diaObjectId, nDiaSources, "\
        "scisql_nanojanskyToAbMag(rPSFluxMin) AS rMagMax, "\
        "scisql_nanojanskyToAbMag(rPSFluxMax) AS rMagMin, "\
        "scisql_nanojanskyToAbMag(gPSFluxMax) AS gMagMin, "\
        "scisql_nanojanskyToAbMag(iPSFluxMax) AS iMagMin "\
        "FROM dp02_dc2_catalogs.DiaObject "\
        "WHERE nDiaSources > "+str(snia_nDiaSources_min)+" "\
        "AND nDiaSources < "+str(snia_nDiaSources_max)+" "\
        "AND scisql_nanojanskyToAbMag(rPSFluxMax) > "+str(snia_peak_mr_min)+" "\
        "AND scisql_nanojanskyToAbMag(rPSFluxMax) < "+str(snia_peak_mr_max)+" "\
        "AND scisql_nanojanskyToAbMag(gPSFluxMax) < "+str(snia_peak_mg_max)+" "\
        "AND scisql_nanojanskyToAbMag(iPSFluxMax) < "+str(snia_peak_mi_max)+" "\
        "AND scisql_nanojanskyToAbMag(rPSFluxMin)"\
        " - scisql_nanojanskyToAbMag(rPSFluxMax) < "+str(snia_ampl_mr_max)+" "\
        "AND scisql_nanojanskyToAbMag(rPSFluxMin)"\
        " - scisql_nanojanskyToAbMag(rPSFluxMax) > "+str(snia_ampl_mr_min)

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED

DiaObjs = job.fetch_result().to_table()

Calculate the r-band magnitude amplitude as the difference between the brightest and faintests difference-image magnitude, and add it to the table.

DiaObjs['rMagAmp'] = DiaObjs['rMagMin'] - DiaObjs['rMagMax']

Option: display the table.

# DiaObjs

Plot histograms of brightest magnitude and lightcurve amplitude

Plot histograms to characterize the brightest detected r-band magnitudes (left) and the r-band amplitudes (right) for the DiaObjects in this initial sample of potentially SNIa-like time-domain objects.

Notice that the brightest r-band magnitude distribution (left) is bimodal, and that the distribution of r-band amplitudes declines until about 3.5 mag and then extends as a tail of brighter time-variable objects to amplitudes of 5.5 mag. The underlying cause of both of these effects is likely the fact that the only time-variable sources in DC2 are stars (mostly brighter) and Type Ia supernovae (mostly fainter) -- but confirming this is left as an exercise for the learner.

fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)

ax[0].hist(DiaObjs['rMagMin'], bins=20, color=plot_filter_colors['r'])
ax[0].set_xlabel('Brightest Detected r-band Magnitude')
ax[0].set_ylabel('Number of Potential SNIa')

ax[1].hist(DiaObjs['rMagAmp'], bins=20, color=plot_filter_colors['r'])
ax[1].set_xlabel('Amplitude in r-band Magnitude')

plt.tight_layout()
plt.show()

Figure 8: At left, the distribution of the brightest difference-image magnitude for all the potential SNIa. At right, the distribution of the difference-image magnitude amplitude for all the potential SNIa.

3.3. Calculate lightcurve duration and identify potential SNIa

The lightcurve duration -- time between the first and last detected DiaSource in any filter -- is not included in the DiaObject table. It is calculated below, using all of the DiaSources for each DiaObject.

Time is reported in the DiaSource table as midPointTai, which is in the SI unit of "TAI" (International Atomic Time), and is presented in days (in particular, as "Modified Julian Days").

First, create a string list of the diaObjectId that is formatted like "(1249546586455802795, 1248684569339625641, 1248684569339626125)".

all_ids = DiaObjs['diaObjectId']
temp_list = '('
for i, id in enumerate(all_ids):
    temp_list = temp_list + str(id)
    if i < len(all_ids)-1:
        temp_list = temp_list + ', '
    elif i == len(all_ids)-1:
        temp_list = temp_list + ')'

Option to print the 1000-element long list.

# print(temp_list)

Retrieve the midPointTai for all observations (diaSources) for all the diaObjects.

query = "SELECT diaObjectId, midPointTai "\
        "FROM dp02_dc2_catalogs.DiaSource "\
        "WHERE diaObjectId IN " + temp_list

# print(query)

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED

results = job.fetch_result().to_table()

Calculate the durations for each diaObject and add them to the DiaObjs table.

DiaObjs['duration'] = numpy.zeros(len(DiaObjs), dtype='float')

for j, DiaObjId in enumerate(DiaObjs['diaObjectId']):
    tx = numpy.where(results['diaObjectId'] == DiaObjId)[0]
    DiaObjs['duration'][j] = numpy.max(results['midPointTai'][tx]) - numpy.min(results['midPointTai'][tx])
    del tx

Clean up.

del results

Select only DiaObjects that have lightcurve durations within the specified range for SNIa.

tx = numpy.where((DiaObjs['duration'] > snia_duration_min)
                 & (DiaObjs['duration'] < snia_duration_max))[0]
print(len(tx))

42

Plot a histogram of lightcurve durations, and a scatter plot of duration vs. amplitude, for the potential SNIa

Plot a histogram of lightcurve durations (left) and a plot of duration versus r-band amplitude (right) to further characterize this sample of potential SNIa.

Notice the spread in duration, and how it is correlated with lightcurve amplitude, as expected for SNIa.

Tracking down the origin of the two outliers is left as an exercise for the learner: if they are truly SNIa observed for >100 days, then they should have a brighter amplitude.

fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharey=False, sharex=False)

ax[0].hist(DiaObjs['duration'][tx], bins=20, color='gray')
ax[0].set_xlabel('Lightcurve Duration [days]')
ax[0].set_ylabel('Number of Potential SNIa')

ax[1].plot(DiaObjs['duration'][tx], DiaObjs['rMagAmp'][tx], 'o', color='gray')
ax[1].set_xlabel('Lightcurve Duration [days]')
ax[1].set_ylabel('Amplitude in r-band [mag]')

plt.tight_layout()
plt.show()

Figure 9: At left, the distribution of the lightcurve duration, in days, for potential SNIa. At right, the lightcurve duration versus r-band amplitude.

3.4. Plot multi-band lightcurves for the first 20 potential SNIa

Notice that all of these lightcurves do resemble Type Ia supernovae! That is a promising sign that the methodology of using the DiaObject summary parameters, along with the lightcurve duration (which must be calculate separately), can identify an initial sample of potential SNIa.

Keep in mind that this is not necessarily a pure or complete sample of SNIa, just a sample of potential well-sampled low-redshift SNIa - decent enough to use as a demonstration, or a starting point for a more rigorous classification process. Different science goals will have different requirements on sample identification.

Below, the lightcurves plotted are based on the signal-to-noise ratio > 5 detections from the DiaSource table only. In Section 4, lightcurves based on forced photometry are used instead.

fig, ax = plt.subplots(5, 4, figsize=(10, 10), sharey=False, sharex=False)

x = 0

for i in range(5):
    for j in range(4):
        results = service.search("SELECT ra, decl, diaObjectId, diaSourceId, "
                                 "filterName, midPointTai, "
                                 "scisql_nanojanskyToAbMag(psFlux) AS psAbMag "
                                 "FROM dp02_dc2_catalogs.DiaSource "
                                 "WHERE diaObjectId = "+str(DiaObjs['diaObjectId'][tx[x]]))
        results = results.to_table()

        for f, filt in enumerate(plot_filter_labels):
            fx = numpy.where(results['filterName'] == filt)[0]
            ax[i, j].plot(results['midPointTai'][fx], results['psAbMag'][fx],
                          plot_filter_symbols[filt], ms=10, mew=0, alpha=0.5,
                          color=plot_filter_colors[filt])
            del fx

        ax[i, j].invert_yaxis()
        ax[i, j].set_title(DiaObjs['diaObjectId'][tx[x]])

        if i == 4:
            ax[i, j].xaxis.set_label_text('MJD (days)')
        if j == 0:
            ax[i, j].yaxis.set_label_text('mag')

        x += 1
        del results

plt.tight_layout()
plt.show()

Figure 10: Multi-band lightcurves for 20 potential SNIa. Legend: $u$, blue circle; $g$, green triangle; $r$, dark red inverted triangle; $i$, golden square; $z$, pink star; $y$, brown pentagram.

4. Forced photometry lightcurves

Above, the DiaSource table is used to plot lightcurves, but it only includes SNR>5 detections in the difference images. If your science goal requires lower-SNR measurements, use the forced photometry in the ForcedSourceOnDiaObjects table instead, as demonstrated below.

Let's say there is a specific object whose lightcurve looks particularly interesting. In order to explore the flux of the object in all sets of differences images (not just those with detected sources), the ForcedSourceOnDiaObject catalog can be utilized.

ForcedSourceOnDiaObject Table: point-source forced-photometry measurements on individual single-epoch visit images and difference images, based on and linked to the entries in the DiaObject table.

For more information, see the schema of the ForcedSourceOnDiaObject catalog.

Option: print all of the available column names and descriptions for the ForcedSourceOnDiaObject table.

# results = service.search("SELECT column_name, description from TAP_SCHEMA.columns "
#                          "WHERE table_name = 'dp02_dc2_catalogs.ForcedSourceOnDiaObject'")
# for column_name, description in zip(results['column_name'], results['description']):
#     print(column_name, '-', description)

4.1. Retrieve photometry for a DiaObject

As an example, let's look at the forced photometry on the difference images of the object 1250953961339360185, which has an interesting lightcurve. Note that it is unlikely this object is shown in the 20 SNIa candidates above since they will change every time the search is run. That object should however be in the list of 1000 SNIa candidates that we started with.

DiaObjID = 1250953961339360185
print(DiaObjID)

1250953961339360185

Here, we will define the query to extract the information from the ForcedSourceOnDiaObject table on our DiaObject of interest in order to construct a lightcurve that includes forced photometry (i.e., sub-threshold measurements with SNR<5).

Below, a JOIN to the CcdVisit table is used to obtain the exposure time mid-point in modified julian date (MJD) format (expMidptMJD), which is necessary to construct the lightcurve.

The ForcedSourceOnDiaObject contains forced photometry on both the difference image (psfDiffFlux, psfDiffFluxErr) and the processed visit image (PVI), also called the "direct" image (psfFlux, psfFluxErr). Both types of fluxes are retrieved below for our DiaObject of interest.

query = "SELECT fsodo.coord_ra, fsodo.coord_dec, "\
        "fsodo.diaObjectId, fsodo.ccdVisitId, fsodo.band, "\
        "fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, "\
        "fsodo.psfFlux, fsodo.psfFluxErr, "\
        "cv.expMidptMJD "\
        "FROM dp02_dc2_catalogs.ForcedSourceOnDiaObject as fsodo "\
        "JOIN dp02_dc2_catalogs.CcdVisit as cv ON cv.ccdVisitId = fsodo.ccdVisitId "\
        "WHERE fsodo.diaObjectId = "+str(DiaObjID)
print(query)

SELECT fsodo.coord_ra, fsodo.coord_dec, fsodo.diaObjectId, fsodo.ccdVisitId, fsodo.band, fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, fsodo.psfFlux, fsodo.psfFluxErr, cv.expMidptMJD FROM dp02_dc2_catalogs.ForcedSourceOnDiaObject as fsodo JOIN dp02_dc2_catalogs.CcdVisit as cv ON cv.ccdVisitId = fsodo.ccdVisitId WHERE fsodo.diaObjectId = 1250953961339360185

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED

Save the results in a table "FrcdSrc" (short for forced source, to distinguish it from the results from DiaSource saved as DiaSrc, below).

FrcdSrc = job.fetch_result().to_table()

# FrcdSrc

Also retrieve the DiaSource measurements for the DiaObject of interest. Recall that this will only return flux measurements from the difference images with SNR>5 detection.

query = "SELECT ra, decl, diaObjectId, diaSourceId, psFlux, psFluxErr, "\
        "filterName, midPointTai, SNR, ccdVisitId, "\
        "scisql_nanojanskyToAbMag(psFlux) AS psAbMag "\
        "FROM dp02_dc2_catalogs.DiaSource "\
        "WHERE diaObjectId = "+str(DiaObjID)
print(query)

SELECT ra, decl, diaObjectId, diaSourceId, psFlux, psFluxErr, filterName, midPointTai, SNR, ccdVisitId, scisql_nanojanskyToAbMag(psFlux) AS psAbMag FROM dp02_dc2_catalogs.DiaSource WHERE diaObjectId = 1250953961339360185

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)

Job phase is COMPLETED

DiaSrc = job.fetch_result().to_table()

Unlike the DiaSource catalog, the ForcedSourceOnDiaObject catalog does not provide the SNR for the flux measurements. Since it provides a useful diagnostic, we define 'diffSNR' for our object by the ratio of the absolute value of 'psfDiffFlux' and 'psfDiffFluxErr.'

FrcdSrc['DiffSNR'] = abs(FrcdSrc['psfDiffFlux']/FrcdSrc['psfDiffFluxErr'])

Option to show that the DiffSNR column has been added to our "FrcdSrc" table.

# FrcdSrc

4.2. Plot the lightcurves

For our DiaObject of interest, compare the lightcurves that include forced photometry from the ForcedSourceOnDiaObject table (open markers, left plots) with lightcurves that include only SNR>5 detections from the DiaSource table (filled markers, right plots).

Notice: The forced photometry lightcurves (left) include flux measurements from before and after the source was identified in the DIAsource catalogs. The shaded gray regions indicate the dates spanned by the detections only (DiaSources only).

Define the MJD ranges for the "full" lightcurve from the forced photometry, and for the "detected" part of the lightcurve from the DiaSource table.

mjd_full_start = numpy.min(FrcdSrc['expMidptMJD'])
mjd_full_end = numpy.max(FrcdSrc['expMidptMJD'])

mjd_det_start = numpy.min(DiaSrc['midPointTai'])
mjd_det_end = numpy.max(DiaSrc['midPointTai'])

4.2.1. Lightcurve for the full survey

fig, ax = plt.subplots(1, 2, figsize=(12, 4), sharey=True, sharex=False)

for f, filt in enumerate(plot_filter_labels):
    fx = numpy.where(FrcdSrc['band'] == filt)[0]
    ax[0].plot(FrcdSrc['expMidptMJD'][fx], FrcdSrc['psfDiffFlux'][fx],
               plot_filter_symbols[filt], ms=10, mew=2, mec=plot_filter_colors[filt],
               alpha=0.5, color='none', label=filt)
    del fx

ax[0].set_xlabel('Modified Julian Date')
ax[0].set_ylabel('Difference-Image Flux')
ax[0].set_title('Forced Photometry (ForcedSourceOnDiaObject)')
ax[0].set_xlim(mjd_full_start, mjd_full_end)
ax[0].legend(loc='upper left')
ax[0].axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')

for f, filt in enumerate(plot_filter_labels):
    fx = numpy.where(DiaSrc['filterName'] == filt)[0]
    ax[1].plot(DiaSrc['midPointTai'][fx], DiaSrc['psFlux'][fx],
               plot_filter_symbols[filt], ms=10, mew=0, color=plot_filter_colors[filt],
               alpha=0.5, label=filt)
    del fx

ax[1].set_xlabel('Modified Julian Date')
ax[1].set_title('SNR>5 Detections (DiaSource)')
ax[1].set_xlim(mjd_full_start, mjd_full_end)
ax[1].legend(loc='upper left')
ax[1].axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')
plt.subplots_adjust(wspace=.0)

plt.show()

Figure 11: The multi-band lightcurve for the chosen diaObject based on forced photometry in all difference images (left) and measurement photometry for signal-to-noise ratio $>$ 5 difference-image detections only (right).

4.2.2. Lightcurve for epochs with SNR>5 detections

fig, ax = plt.subplots(1, 2, figsize=(12, 4), sharey=True, sharex=False)

for f, filt in enumerate(plot_filter_labels):
    fx = numpy.where(FrcdSrc['band'] == filt)[0]
    ax[0].plot(FrcdSrc['expMidptMJD'][fx], FrcdSrc['psfDiffFlux'][fx],
               plot_filter_symbols[filt], ms=10, mew=2, mec=plot_filter_colors[filt],
               alpha=0.5, color='none', label=filt)
    del fx
ax[0].set_xlabel('Modified Julian Date')
ax[0].set_ylabel('Difference-Image Flux')
ax[0].set_title('Forced Photometry (ForcedSourceOnDiaObject)')
ax[0].set_xlim(mjd_det_start-50, mjd_det_end+50)
ax[0].legend(loc='upper left')
ax[0].axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')

for f, filt in enumerate(plot_filter_labels):
    fx = numpy.where(DiaSrc['filterName'] == filt)[0]
    ax[1].plot(DiaSrc['midPointTai'][fx], DiaSrc['psFlux'][fx],
               plot_filter_symbols[filt], ms=10, mew=0, color=plot_filter_colors[filt],
               alpha=0.5, label=filt)
    del fx

ax[1].set_xlabel('Modified Julian Date')
ax[1].set_title('SNR>5 Detections (DiaSource)')
ax[1].set_xlim(mjd_det_start-50, mjd_det_end+50)
ax[1].legend(loc='upper left')
ax[1].axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')
plt.subplots_adjust(wspace=.0)

plt.show()

Figure 12: As in Figure 11, but zoomed-in on the dates for which the diaObject was detected.

Notice: In the above plots, the forced photometry (psfDiffFlux from ForcedSourceOnDiaObject, open symbols, left panels) is nearly identical to the detection photometry (psFlux from DiaSource, filled symbols, right panels), but the ForcedSourceOnDiaObject contains additional measurements from epochs in which the source in the difference image had a SNR<5. No such cutoff is applied to the ForcedSourceOnDiaObject catalog, in which photometry from all visits is included.

4.3. Identify "missing measurements" from SNR<5 epochs

Within the "gray zone" above (defined by epochs in which the source was detected with SNR>5 in the difference images), identify how many measurements there are in the DiaSource and ForcedSourceOnDiaObject tables, and compare.

cov = numpy.where((FrcdSrc['expMidptMJD'] >= mjd_det_start)
                  & (FrcdSrc['expMidptMJD'] <= mjd_det_end))

print('Within the date range of SNR>5 detections, there were:')
print(len(DiaSrc), 'measurements in DiaSource catalog of %.i' % (DiaObjID))
print(len(FrcdSrc[cov]), 'measurements in ForcedSourceOnDiaObject catalog of %.i' % (DiaObjID))

Within the date range of SNR>5 detections, there were:
27 measurements in DiaSource catalog of 1250953961339360185
40 measurements in ForcedSourceOnDiaObject catalog of 1250953961339360185

Identify only the SNR<5 measurements from the ForcedSourceOnDiaObject table by looking for unique "ccdVisitId" values that are not in the DiaSource results. Note that we are still applying the condition where the measurements are taken within the date range of the DiaSource detections.

First, identify the unique visit IDs.

UniqueVID = set(FrcdSrc[cov]['ccdVisitId']) - set(DiaSrc['ccdVisitId'])

Next, construct a table of ForcedSourceOnDiaObject "missing measurements" that are not in the DiaSource catalog

UniqueObj = []
for x in UniqueVID:
    UniqueObj += numpy.where(FrcdSrc['ccdVisitId'] == x)
UniqueObj = numpy.concatenate(UniqueObj)

Inspect the table of "missing measurements".

FrcdSrc[UniqueObj]

Table length=13

coord_ra	coord_dec	diaObjectId	ccdVisitId	band	psfDiffFlux	psfDiffFluxErr	psfFlux	psfFluxErr	expMidptMJD	DiffSNR
deg	deg				nJy	nJy	nJy	nJy	d
float64	float64	int64	int64	object	float64	float64	float64	float64	float64	float64
60.2901401	-44.142051	1250953961339360185	960109025	g	397.7962834	109.4827468	360.7160608	107.5246745	60993.0627302	3.63341526429441
60.2901401	-44.142051	1250953961339360185	948985095	u	-226.9062832	223.3410915	-253.8417292	216.7644011	60970.1579432	1.0159629904020595
60.2901401	-44.142051	1250953961339360185	936368072	y	4098.2779628	1452.0982926	3700.1518782	1435.3723413	60952.1606072	2.8223144284964228
60.2901401	-44.142051	1250953961339360185	936421079	y	2472.0998484	1413.2807205	2400.1562592	1404.5707216	60952.1842852	1.7491923667687224
60.2901401	-44.142051	1250953961339360185	936422155	y	4647.6798641	1429.1860639	4530.4114916	1420.2903694	60952.1847302	3.2519767590073547
60.2901401	-44.142051	1250953961339360185	936410027	y	5923.5763416	1407.7463626	5710.2399879	1402.4808996	60952.1793702	4.207843471646132
60.2901401	-44.142051	1250953961339360185	937668075	y	1975.4528672	1362.9556253	1879.4946357	1353.4784136	60954.3511812	1.4493889826861994
60.2901401	-44.142051	1250953961339360185	936428011	y	2615.2255212	1419.1242198	2424.3486295	1412.3991488	60952.1874072	1.8428446817492614
60.2901401	-44.142051	1250953961339360185	948940154	u	452.2146039	248.4894382	456.5023539	241.2969625	60970.1377562	1.8198544259093583
60.2901401	-44.142051	1250953961339360185	936431155	y	3180.2418489	1401.5102644	3107.1975747	1388.7839996	60952.1887502	2.2691534480209405
60.2901401	-44.142051	1250953961339360185	963902134	u	151.6832421	207.2712358	117.0506357	201.4125375	60998.0695872	0.7318103812839813
60.2901401	-44.142051	1250953961339360185	956243158	y	3042.9740454	1353.8394814	2913.5862492	1345.4704937	60980.0924622	2.247662361163579
60.2901401	-44.142051	1250953961339360185	961630074	u	-29.0326013	220.3748454	-46.4890129	212.8287779	60995.0636662	0.13174190206374614

Finally, verify why these "missing measurements" are not in the DiaSource catalog (their SNR is <5) by printing the maximum SNR.

numpy.max(FrcdSrc[UniqueObj]['DiffSNR'])

4.207843471646132

4.4. Investigate potential SNR<5 variability

Going back to the forced photometry lightcurve of our object of interest, it looked like there may have been some past variability in the z band (purple star) at around MJD 60250. Let's investigate this point. This is for demonstration purposes only since there are no simulated precursor outbursts for the supernovae in DP0.2.

Let's check if it is at least greater than 3$\sigma$ of the baseline (median) measurements using "sigma_clipped_stats".

Define the filter in which flaring (precursor activity) is investigated.

flare_filt = 'z'

Specify the median and standard deviation in the "flare" filter band for the flux measurements.

_, med, std = sigma_clipped_stats(FrcdSrc['psfDiffFlux'][numpy.where(FrcdSrc['band'] == flare_filt)[0]])

print('Median psfDiffFlux for all z-band forced sources: ', med, FrcdSrc['psfDiffFlux'].unit)
print('Standard deviation in psfDiffFlux for all z-band forced sources: ',std, FrcdSrc['psfDiffFlux'].unit)

Median psfDiffFlux for all z-band forced sources:  -44.376622 nJy
Standard deviation in psfDiffFlux for all z-band forced sources:  568.5471832642495 nJy

Define where any outlier variability (above 3$\sigma$) might be in the "flare" band, BEFORE or AFTER the "gray zone" defined by the epochs in which our DiaObject of interest was detected with SNR>5 in any filter.

First, index the flare-band epochs "before" and "after" within FrcdSrc.

ex = numpy.where(((FrcdSrc['expMidptMJD'] < mjd_det_start)
                  | (FrcdSrc['expMidptMJD'] > mjd_det_end))
                 & (FrcdSrc['band'] == flare_filt))[0]
print(len(ex))

77

Next, index which of these is >3$\sigma$.

ix = numpy.where(FrcdSrc['psfDiffFlux'][ex] > med+3*std)[0]
print('Number of >3𝜎 z-band measurements (outside SNR>5 detections):',len(ix))

Number of >3𝜎 z-band measurements (outside SNR>5 detections): 1

Inspect the single >3$\sigma$ measurement in the forced photometry that occurred outside of the "gray zone".

As a reminder, 3$\sigma$ is defined from the distribution of forced photometry measurements and is not the photometric signal-to-noise ratio (i.e., DiffSNR).

FrcdSrc[ex[ix]]

Table length=1

coord_ra	coord_dec	diaObjectId	ccdVisitId	band	psfDiffFlux	psfDiffFluxErr	psfFlux	psfFluxErr	expMidptMJD	DiffSNR
deg	deg				nJy	nJy	nJy	nJy	d
float64	float64	int64	int64	object	float64	float64	float64	float64	float64	float64
60.2901401	-44.142051	1250953961339360185	476269099	z	2517.1316512	1251.6465398	2667.6647459	1241.6250563	60275.0197362	2.0110562935780663

Statistics check: recall that $\pm 3\sigma$ encompasses 99.7% of values in a normally-distributed random sample. So out of the 77 z-band forced photometry measurements we could have reasonably expected 0.23 at the $\sigma$ > 3 level, and the fact that we found one is consistent with small-number statistics, i.e., consistent with noise.

Again, since this is from DP0.2, we know this "flaring activity" is just noise but we use it to demonstrate how the ForcedSourceOnDIAObject catalog can be applied to search for pre- or post-transient variability.

As a final demonstration, plot the point(s) on the lightcurve that exceed 3$\sigma$ and are outside the DiaSource date coverage. Below, we plot the 3$\sigma$ flux range in red and the apparent variability outlier(s) as a yellow-filled plot marker.

fig, ax = plt.subplots(figsize=(7, 4))

for f, filt in enumerate(plot_filter_labels):
    fx = numpy.where(FrcdSrc['band'] == filt)[0]
    ax.plot(FrcdSrc['expMidptMJD'][fx], FrcdSrc['psfDiffFlux'][fx],
            plot_filter_symbols[filt], ms=10, mew=2, mec=plot_filter_colors[filt],
            alpha=0.5, color='none', label=filt)
    del fx

for i in ix:
    ax.plot(FrcdSrc['expMidptMJD'][ex[i]], FrcdSrc['psfDiffFlux'][ex[i]],
            plot_filter_symbols['z'], ms=14, mew=2, mec='black', alpha=1.0,
            color='yellow', label='z-band >3$\sigma$ FrcdSrc psfDiffFlux')

ax.hlines(y=med, xmin=numpy.min(FrcdSrc['expMidptMJD']), xmax=numpy.max(FrcdSrc['expMidptMJD']),
          linewidth=2, color='black', ls='--', label='median z-band FrcdSrc psfDiffFlux')
ax.axhspan(med, med+3*std, alpha=0.3, color='tab:red', label='stddev of z-band FrcdSrc psfDiffFlux')

ax.set_xlabel('Modified Julian Date')
ax.set_ylabel('Difference-Image Flux')
ax.set_title('Forced Photometry (ForcedSourceOnDiaObject)')
ax.legend(loc='upper left', ncol=2)
ax.set_xlim(mjd_full_start, mjd_full_end)
ax.axvspan(mjd_det_start, mjd_det_end, alpha=0.3, color='gray')
plt.show()

Figure 13: The same lightcurve as in the left panel of Figure 11, but with statistics for the $z$-band forced flux: median (dashed line) and standard deviation (pink).

5. Exercises for the learner

From here, there are a variety of exercises that a learner could undertake.

1. Change the query constraints and identify a sample of high-redshift SNIa.
2. Tighten the query contraints to identify only SNIa with pre-peak observations.
3. Add error bars to the lightcurves. Magnitude errors can be retrieved during a TAP survey with, e.g., scisql_nanojanskyToAbMagSigma(psFlux, psFluxErr) as psAbMagErr.
4. Apply a lightcurve-template fitter, or try a photometric classifier on the sample.
5. Investigate the forced source photometry lightcurves of other SNIa candidates.