DIA Analysis: Supernovae from the Run 1.2p Test

Michael Wood-Vasey Last Verified to Run: 2019-07-17

After completing this Notebook, the user will be able to

  1. Get a list of simulated SN from the Run 1.2p truth catalog
  2. Select SNe within the test patch in the Run 1.2p DIA Test
  3. Plot lightcurves from the DIA analysis
  4. Plot lightcurves from the truth catalog
  5. Compare the input and recovered lightcurves.

See the Truth GCR Variables for a basic overview of the Truth information and how to access it for variables. https://github.com/LSSTDESC/DC2-analysis/blob/master/tutorials/truth_gcr_variables.ipynb

In [1]:
# Inject gcr-catalogs that supports DIA source into path.
import os
import math
import sys

import numpy as np
import pandas as pd

from astropy.coordinates import SkyCoord
import astropy.units as u
In [2]:
import lsst.afw.display as afwDisplay
import lsst.afw.geom as afwGeom
from lsst.daf.persistence import Butler
from lsst.geom import SpherePoint
import lsst.geom
In [3]:
import GCRCatalogs
In [4]:
%matplotlib inline

import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
In [5]:
repo = '/global/cscratch1/sd/rearmstr/new_templates/diffim_template'
In [6]:
butler = Butler(repo)
In [7]:
diaSrc = GCRCatalogs.load_catalog('dc2_dia_source_run1.2p_test')
diaObject = GCRCatalogs.load_catalog('dc2_dia_object_run1.2p_test')
truth_cat = GCRCatalogs.load_catalog('dc2_truth_run1.2_variable_summary')
truth_lc = GCRCatalogs.load_catalog('dc2_truth_run1.2_variable_lightcurve')

(We presently will get a warning from the catalog reader in the initalization above because there is no u-band in the subtractions.)

Select SNe from the truth catalog

In [8]:
columns = ['ra', 'dec', 'redshift', 'uniqueId', 'galaxy_id', 'sn']
truth_all_sne = pd.DataFrame(truth_cat.get_quantities(columns, filters=[f'sn == 1']))

You'll get a dataset.value deprecation warning. Don't worry about this. The Data Access Team will fix this someday.

In [9]:
print(f'We found {len(truth_all_sne)} simulated SNe.')
We found 76689 simulated SNe.

Select SNe from the truth catalog in the DIA test region

In [10]:
tract = 4849
patch = (6, 6)
In [11]:
skymap = butler.get('deepCoadd_skyMap')
tract_info = skymap[tract]
In [12]:
foo = tract_info.getPatchInfo(patch)
In [13]:
bar = foo.getOuterSkyPolygon(tract_info.getWcs())
In [14]:
tract_box = afwGeom.Box2D(tract_info.getBBox())
tract_pos_list = tract_box.getCorners()
wcs = tract_info.getWcs()
corners = wcs.pixelToSky(tract_pos_list)
corners = np.array([[c.getRa().asDegrees(), c.getDec().asDegrees()] for c in corners])
/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/ipykernel/__main__.py:1: FutureWarning: Call to deprecated function (or staticmethod) Box2D. (Replaced by lsst.geom.Box2D (will be removed before the release of v20.0))
  if __name__ == '__main__':
In [15]:
print(corners)
[[ 54.72776196 -29.78294629]
 [ 52.93572545 -29.78294586]
 [ 52.949112   -28.22767319]
 [ 54.71437636 -28.2276736 ]]
In [16]:
ra = corners[:, 0]
dec = corners[:, 1]
min_ra, max_ra = np.min(ra), np.max(ra)
min_dec, max_dec = np.min(dec), np.max(dec)
In [17]:
area_cut = [f'ra > {min_ra}', f'ra < {max_ra}', f'dec > {min_dec}', f'dec < {max_dec}']
sn_cut = ['sn == 1']
all_cuts = area_cut + sn_cut
In [18]:
columns = ['ra', 'dec', 'redshift', 'uniqueId', 'galaxy_id', 'sn']
truth_sne = pd.DataFrame(truth_cat.get_quantities(columns, filters=all_cuts))
In [19]:
print(f'We found {len(truth_sne)} simulated SNe.')
We found 8508 simulated SNe.
In [20]:
avg_dec = (min_dec + max_dec)/2
size = 8
dec_size = size
ra_size = dec_size * np.cos(np.deg2rad(avg_dec))
aspect_ratio = dec_size / ra_size

# fig = plt.figure(figsize=(size, size))
ax = plt.gca()
ax.set_aspect(aspect_ratio)

patch_region = Polygon(corners, color='red', fill=False)

ax.scatter(truth_sne['ra'], truth_sne['dec'])
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.set_xlim(plt.xlim()[::-1])
ax.add_patch(patch_region);
ax.set_title(f'tract: {tract}, patch: {patch}');

A Lightcurve

Search for diaObjects that match input simulated SNe.

In [21]:
i = 0
sn = truth_sne.iloc[0]
In [22]:
print(sn)
galaxy_id    3.010259e+11
dec         -2.838755e+01
sn           1.000000e+00
redshift     2.376888e-01
ra           5.300728e+01
uniqueId     3.082506e+14
Name: 0, dtype: float64

Oops, we need to fix up the dtypes somewhere. Those shouldn't all be floats.

Get the lightcurve

In [23]:
columns = ['obshistid', 'mjd', 'mag', 'filter']

sn_lc = pd.DataFrame(truth_lc.get_quantities(columns,
                                             native_filters=[f'uniqueId == {sn["uniqueId"]}']))
In [24]:
sn_lc.rename(columns={'filter': 'filter_code'}, inplace=True)
In [25]:
# Translate filter codes to filter names
filter_names = ['u', 'g', 'r', 'i', 'z', 'y']
sn_lc['filter'] = [filter_names[f] for f in sn_lc['filter_code']]
sn_lc = sn_lc.sort_values('mjd')
In [26]:
def plot_lightcurve(df, plot='mag', flux_col_names=None,
                    title=None, marker='o', linestyle='none',
                    colors=None, label_prefix='',
                    **kwargs):
    """Plot a lightcurve from a DataFrame.
    """
    # At lexigraphical order, if not wavelength order.
    # Assume fixed list of filters.
    filter_order = ['u', 'g', 'r', 'i', 'z', 'y']

    if colors is None:
        colors = {'u': 'violet', 'g': 'indigo', 'r': 'blue', 'i': 'green', 'z': 'orange', 'y': 'red'}
    
    if flux_col_names is not None:
        flux_col, flux_err_col = flux_col_names
    else:
        if plot == 'flux':
            flux_col = 'psFlux'
            flux_err_col = 'psFluxErr'
        else:
            flux_col = 'mag'
            flux_err_col = 'mag_err'
        
    for filt in filter_order:
        this_filter = df.query(f'filter == "{filt}"')
        if this_filter.empty:
            continue
        # This if sequence is a little silly.
        plot_kwargs = {'linestyle': linestyle, 'marker': marker, 'color': colors[filt],
                       'label': f'{label_prefix} {filt}'}
        plot_kwargs.update(kwargs)

        if flux_err_col in this_filter.columns:
            plt.errorbar(this_filter['mjd'], this_filter[flux_col], this_filter[flux_err_col],
                         **plot_kwargs)
                        
        else:
            if marker is None:
                plt.plot(this_filter['mjd'], this_filter[flux_col], **plot_kwargs)

            else:
                plot_kwargs.pop('linestyle')
                plt.scatter(this_filter['mjd'], this_filter[flux_col], **plot_kwargs)



    plt.xlabel('MJD')

    if plot == 'flux':
        plt.ylabel('psFlux [nJy]')
    else:
        # Ensure that y-axis decreases as one goes up
        # Because plot_lightcurve could be called several times on the same axis,
        # simply inverting is not correct.  We have to reverse a sorted list.
        plt.ylim(sorted(plt.ylim(), reverse=True))
        plt.ylabel('mag [AB]')

    if title is not None:
        plt.title(title)
    plt.legend()
In [27]:
plt.figure(figsize=(12, 8))
plot_lightcurve(sn_lc, plot='mag', linestyle='-', marker=None)

Match to DIAObject Catalog

In [28]:
ra, dec = sn['ra'], sn['dec']
print(ra, dec)
53.00727812651823 -28.387552632959427
In [29]:
# Match on RA, Dec
sn_position = SkyCoord(sn['ra'], sn['dec'], unit='deg')

diaObject_cat = diaObject.get_quantities(['ra', 'dec', 'diaObjectId'])
diaObject_positions = SkyCoord(diaObject_cat['ra'], diaObject_cat['dec'], unit='deg')

idx, sep2d, _ = sn_position.match_to_catalog_sky(diaObject_positions)
print(f'Index: {idx} is {sep2d.to(u.arcsec)[0]:0.6f} away')
Index: 22 is 0.028836 arcsec away
In [30]:
diaObjectId = diaObject_cat['diaObjectId'][idx]

How did we do with the lightcurve?

In [31]:
sn_lc
Out[31]:
filter_code mjd mag obshistid filter
0 0 60554.288025 26.432639 657401 u
1 0 60555.258471 25.581090 658296 u
2 0 60556.250832 25.019064 659233 u
3 0 60557.246086 24.602021 660201 u
4 0 60558.249432 24.266763 660882 u
5 0 60559.245016 23.986753 661527 u
6 2 60567.219180 21.397669 664557 r
7 1 60567.228879 21.461969 664577 g
8 3 60567.234411 21.431454 664587 i
9 4 60567.244111 21.719140 664607 z
10 5 60567.256310 21.978863 664633 y
11 2 60578.271042 21.354127 671659 r
12 1 60578.280324 21.694316 671678 g
13 3 60578.285857 21.346081 671688 i
14 4 60578.295556 21.825918 671708 z
15 5 60578.307755 22.263662 671734 y
16 2 60581.174302 21.450943 674105 r
17 3 60581.189152 21.472849 674134 i
18 4 60581.198851 21.983958 674154 z
19 5 60581.211050 22.407560 674180 y
20 0 60582.204745 24.257999 675069 u
21 0 60583.180252 24.374555 675910 u
22 0 60584.176755 24.494031 676739 u
23 0 60585.171706 24.608092 677595 u
24 0 60586.172558 24.727928 678467 u
25 0 60587.166237 24.855251 679317 u
26 0 60588.163054 24.992185 680167 u
27 2 60593.164657 22.093300 684469 r
28 1 60593.173940 23.101931 684488 g
29 3 60593.179472 22.013917 684498 i
... ... ... ... ... ...
43 4 60608.151733 22.571821 696296 z
44 5 60608.163933 22.559708 696322 y
45 2 60611.119189 23.094521 698622 r
46 1 60611.128471 24.440580 698641 g
47 3 60611.134004 22.555700 698651 i
48 4 60611.143703 22.708548 698671 z
49 5 60611.155902 22.649910 698697 y
50 0 60613.097971 26.392147 700187 u
51 0 60614.092032 26.435882 700968 u
52 0 60615.087792 26.482756 701754 u
53 0 60616.088773 26.532553 702553 u
54 0 60617.085683 26.583687 703320 u
55 0 60620.131032 26.722439 704240 u
56 2 60621.122783 23.508854 704746 r
57 1 60621.132065 24.736442 704765 g
58 3 60621.137598 22.996806 704775 i
59 4 60621.147297 23.206187 704795 z
60 5 60621.159496 23.115017 704821 y
61 2 60624.078591 23.592920 706923 r
62 3 60624.093405 23.095192 706952 i
63 4 60624.103104 23.344373 706972 z
64 5 60624.115303 23.258472 706998 y
65 2 60627.276161 23.667327 709647 r
66 3 60627.291393 23.191451 709677 i
67 4 60627.301092 23.484316 709697 z
68 5 60627.313291 23.379125 709723 y
69 2 60632.036052 23.747787 712881 r
70 3 60632.051295 23.312274 712911 i
71 4 60632.060994 23.667778 712931 z
72 5 60632.073193 23.491572 712957 y

73 rows × 5 columns

In [32]:
# We can't use a direct filters = match in the GCR wrapper for the diaSrc table.
# So we have to use a lambda function here to match the ID
dia_lc = pd.DataFrame(diaSrc.get_quantities(['visit', 'mjd', 'psFlux', 'psFluxErr', 'mag', 'mag_err', 'filter'],
                                            filters=[(lambda x: x == diaObjectId, 'diaObjectId')]))
dia_lc = dia_lc.sort_values('mjd')                    
In [33]:
plt.figure(figsize=(12, 8))
plot_lightcurve(sn_lc, plot='mag', linestyle='-', marker=None, label_prefix='Sim')
plot_lightcurve(dia_lc, plot='mag', label_prefix='DIA')

sim_date_range = [min(sn_lc['mjd']), max(sn_lc['mjd'])]
sim_date_delta = sim_date_range[1] - sim_date_range[0]
buffer_fraction = 0.05
plot_date_range = sim_date_range
plot_date_range[0] -= sim_date_delta * buffer_fraction
plot_date_range[1] += sim_date_delta * buffer_fraction

plt.xlim(plot_date_range)
plt.ylim(25, 19);

The shape seems good. But the calibration is magnitudes off. It does look like it's a constant magnitude offset.

In [34]:
sn_lc.query('(60567 < mjd) & (mjd < 60568)')
Out[34]:
filter_code mjd mag obshistid filter
6 2 60567.219180 21.397669 664557 r
7 1 60567.228879 21.461969 664577 g
8 3 60567.234411 21.431454 664587 i
9 4 60567.244111 21.719140 664607 z
10 5 60567.256310 21.978863 664633 y
In [35]:
dia_lc.query('(60567 < mjd) & (mjd < 60568)')
Out[35]:
mjd visit mag_err mag psFlux psFluxErr filter
1 60567.219782 664557 0.009033 19.610946 51954.298318 430.379216 r
2 60567.235013 664587 0.010489 19.698708 47920.015598 461.672183 i
3 60567.244713 664607 0.012521 19.666216 49375.749296 568.633628 z

The MJD for the same visits (recall visit == obshistid) are slightly different between the truth catalog and the DIA lightcurve:

In [36]:
(60567.219180 - 60567.219782) * 24 * 3600
Out[36]:
-52.0128000061959

Why the 52 second difference?

Let's match on obshistid == visit to get a matched set of dataframes that we can compare.

In [37]:
joint_lc = pd.merge(sn_lc, dia_lc, left_on='obshistid', right_on='visit', suffixes=('_sim', '_dia'))
In [38]:
joint_lc['mjd'] = joint_lc['mjd_dia']
joint_lc['filter'] = joint_lc['filter_sim']
joint_lc['delta_mag'] = joint_lc['mag_dia'] - joint_lc['mag_sim']
In [39]:
joint_lc
Out[39]:
filter_code mjd_sim mag_sim obshistid filter_sim mjd_dia visit mag_err mag_dia psFlux psFluxErr filter_dia mjd filter delta_mag
0 2 60567.219180 21.397669 664557 r 60567.219782 664557 0.009033 19.610946 51954.298318 430.379216 r 60567.219782 r -1.786722
1 3 60567.234411 21.431454 664587 i 60567.235013 664587 0.010489 19.698708 47920.015598 461.672183 i 60567.235013 i -1.732747
2 4 60567.244111 21.719140 664607 z 60567.244713 664607 0.012521 19.666216 49375.749296 568.633628 z 60567.244713 z -2.052925
3 2 60578.271042 21.354127 671659 r 60578.271644 671659 0.008895 19.717774 47085.850480 385.539974 r 60578.271644 r -1.636353
4 1 60578.280324 21.694316 671678 g 60578.280926 671678 0.020358 21.298140 10983.583655 205.931007 g 60578.280926 g -0.396176
5 3 60578.285857 21.346081 671688 i 60578.286459 671688 0.010054 19.572688 53817.660622 498.132297 i 60578.286459 i -1.773393
6 4 60578.295556 21.825918 671708 z 60578.296158 671708 0.012146 19.498688 57613.595367 644.364773 z 60578.296158 z -2.327230
7 5 60578.307755 22.263662 671734 y 60578.308357 671734 0.020698 19.714040 47248.056561 900.609948 y 60578.308357 y -2.549622
8 2 60581.174302 21.450943 674105 r 60581.174904 674105 0.011141 19.933641 38596.129793 395.892130 r 60581.174904 r -1.517302
9 1 60593.173940 23.101931 684488 g 60593.174542 684488 0.074385 23.067606 2152.571205 147.471570 g 60593.174542 g -0.034324
10 5 60593.201370 22.608499 684544 y 60593.201972 684544 0.026534 20.276792 28137.340049 687.511028 y 60593.201972 y -2.331706
11 2 60596.178013 22.270623 686469 r 60596.178615 686469 0.041427 21.337737 10590.229808 404.050880 r 60596.178615 r -0.932887
12 1 60596.187712 23.396907 686489 g 60596.188314 686489 0.159190 23.230260 1853.087931 271.697872 g 60596.188314 g -0.166647
13 3 60596.193244 22.072659 686499 i 60596.193846 686499 0.025884 20.575759 21364.683865 509.285207 i 60596.193846 i -1.496901
14 4 60596.202943 22.472107 686519 z 60596.203545 686519 0.021702 20.212251 29860.670592 596.795208 z 60596.203545 z -2.259856
15 5 60596.215142 22.573407 686545 y 60596.215744 686545 0.033466 20.332999 26717.792021 823.473242 y 60596.215744 y -2.240409
16 2 60605.134462 22.778590 693871 r 60605.135064 693871 0.048667 22.051219 5489.243765 246.036487 r 60605.135064 r -0.727371
17 5 60605.171175 22.502005 693946 y 60605.171777 693946 0.033979 20.480992 23313.259952 729.486718 y 60605.171777 y -2.021012
18 2 60608.127220 22.936279 696247 r 60608.127822 696247 0.059354 22.403800 3967.161602 216.864375 r 60608.127822 r -0.532479
19 1 60608.136502 24.292038 696266 g 60608.137104 696266 0.173746 24.204006 755.784122 120.944561 g 60608.137104 g -0.088033
20 3 60608.142034 22.400639 696276 i 60608.142636 696276 0.035801 21.450470 9545.790759 314.726180 i 60608.142636 i -0.950169
21 4 60608.151733 22.571821 696296 z 60608.152335 696296 0.022675 20.578074 21319.168260 445.118346 z 60608.152335 z -1.993746
22 5 60608.163933 22.559708 696322 y 60608.164535 696322 0.034725 20.659315 19782.167896 632.615017 y 60608.164535 y -1.900393
23 2 60611.119189 23.094521 698622 r 60611.119791 698622 0.071238 22.413331 3932.489517 258.012421 r 60611.119791 r -0.681190
24 1 60611.128471 24.440580 698641 g 60611.129073 698641 0.168484 24.326891 674.906577 104.731225 g 60611.129073 g -0.113689
25 3 60611.134004 22.555700 698651 i 60611.134606 698651 0.044743 21.594204 8362.159167 344.553900 i 60611.134606 i -0.961496
26 5 60611.155902 22.649910 698697 y 60611.156504 698697 0.038859 20.833420 16851.248510 603.064822 y 60611.156504 y -1.816491
27 2 60621.122783 23.508854 704746 r 60621.123385 704746 0.072668 22.687504 3054.910661 204.460426 r 60621.123385 r -0.821351
28 5 60621.159496 23.115017 704821 y 60621.160098 704821 0.060528 21.273221 11238.585132 626.510025 y 60621.160098 y -1.841796
29 3 60624.093405 23.095192 706952 i 60624.094007 706952 0.068156 22.045376 5518.862250 346.433388 i 60624.094007 i -1.049815
30 4 60624.103104 23.344373 706972 z 60624.103706 706972 0.042003 21.195894 12068.193213 466.842536 z 60624.103706 z -2.148479
31 5 60624.115303 23.258472 706998 y 60624.115905 706998 0.071504 21.380532 10180.921902 670.470176 y 60624.115905 y -1.877940
32 3 60627.291393 23.191451 709677 i 60627.291995 709677 0.140924 22.364977 4111.582997 533.663023 i 60627.291995 i -0.826474
33 5 60627.313291 23.379125 709723 y 60627.313893 709723 0.120211 21.638405 8028.566688 888.906650 y 60627.313893 y -1.740720
34 2 60632.036052 23.747787 712881 r 60632.036654 712881 0.136071 22.897444 2517.807350 315.546726 r 60632.036654 r -0.850344
35 4 60632.060994 23.667778 712931 z 60632.061596 712931 0.071836 21.486596 9233.395339 610.905922 z 60632.061596 z -2.181181
In [40]:
plot_lightcurve(joint_lc, flux_col_names=('delta_mag', 'mag_err'))

Hmmm.... so not really a constant offset in magnitude. Nevertheless, I can only suspect there's some calibration or flux interpretation error somewhere.

If one wants to look at the postage stamps of this SN, get started with the examples in dia_source_object_stamp.ipynb

In [ ]: