In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from astropy.table import Table

In [2]:
import sncosmo
# import simulate_lsst as sl

In [3]:
import sncosmo
import simulate_lsst as sl
import numpy as np
from numpy.random import normal, uniform
import matplotlib.pyplot as plt
from astropy.table import Table

In [4]:
# First try with the Simulated LSST files. These files are simulated using SNANA
# so the formatting is different. We can load the files and read them using the next few
# functions

In [5]:
# register LSST bandpasses with names like 'LSST_u', 'LSST_g', 'LSST_r', 'LSST_i' ,'LSST_z', 'LSST_Y' 
sl.getlsstbandpassobjs()

Functions to maninpulate the SNANA format


In [6]:
# Helper Functions to load the SNANA Data
def snanadatafile(snanafileroot, filetype='head', location='./'):
    '''
    obtain the name of the head of phot file of an SNANA simulation and dataset
    
    '''
    import os
    suffix = '_HEAD.FITS'
    if filetype == 'phot':
        suffix = '_PHOT.FITS'
    fname = snanafileroot + suffix
    return os.path.join(location, fname)
    
def loadSNANAData(snanafileroot, location='/.', snids=None, n=None):
    '''
    load a SNANA fits file into a list of `~astropy.Table.table` objects.
    
    
    Parameters
    ----------
    snanafileroot: string, mandatory
        root file name for the SNANA which is the prefix to '_HEAD.FITS', or '_PHOT.FITS'
    location: string, optional defaults to current working directory './' 
        directory where the head and phot files are located
    snids: integer/string, optional defaults to None
        if not None, only SN observations corresponding to SNID snid are loaded
    n: Integer, defaults to None
        if not None, only the first n SN light curves are loaded
        
        
    Returns: data
        list of `~astropy.Table.Table` each Table containing a light curve of a SN. 
        
    ..note: The column names of the SNANA data files are not reformated for SNCosmo use
    '''
    headfile = snanadatafile(snanafileroot, filetype='head', location=location)
    photfile = snanadatafile(snanafileroot, filetype='phot', location=location)
    data = sncosmo.read_snana_fits(head_file=headfile, phot_file=photfile, snids=snids, n=None)
    return data

In [7]:
def addbands(sn, lsstbands, replacement):
    '''
    add a column called 'band' to the `~astropy.Table.Table` by 
    applying the map of lsstbands to replacements to the content
    of a column called 'FLT' 
    
    Parameters
    ----------
    sn: `~astropy.Table.Table` obtained by reading an SNANA light curve
    lsstbands: list of strings, mandatory
        list of strings representing the filters in sn, which can be found
        by `np.unique(sn['FLT'])
    replacements: list of strings, mandatory
        list of strings representing the filters as registered in SNCosmo in
        the same order as lsstbands
        
    Returns
    -------
    `~astropy.Table.Table` with 'FLT' column removed and 'band' column added
    '''
    filterarray = np.zeros(len(sn), dtype='S8')
    for i, flt in enumerate(lsstbands):
        mask = sn['FLT']==flt
        filterarray[mask] = replacement[i]
        band = Table.Column(filterarray, name='band', dtype='S8')
    sn.add_column(band)
    sn.remove_column('FLT')

In [8]:
def reformat_SNANASN(sn, lsstbands=None, replacements=None):
    '''
    reformat an SNANA light curve for use with SNCosmo
    
    Parameters
    ----------
    sn: `~astropy.Table.Table`, mandatory
        representing SNANA light curve
    lsstbands: list of strings, optional defaults to None
        list of unique strings in any of the 'FLT' column of SNANA files
    replacements: list of strings, optional defaults to None
        list of unique strings of the same size as lsstbands, and indexed in the 
        same order representing the keys in the sncosmo.bandpass registry for the
        same filters
    
    
    Returns
    -------
    `astropy.Table.Table` of the SNANA light curve reformatted for SNCosmo 
    '''
    #rename cols to names SNCosmo understands
    sn.rename_column("FLUXCAL",'flux')
    sn.rename_column("FLUXCALERR", 'fluxerr')
    #Add in SNANA magic ZP and sys
    sn["ZP"] = 27.5
    sn["ZPSYS"] = 'ab'
    # sn.rename_column('FLT', 'band')
    
    #Set up a truth dictionary from the metadata to set sim models
    truth ={}
    truth["c"] = sn.meta["SIM_SALT2c"]
    truth["x0"] = sn.meta["SIM_SALT2x0"]*10**(-0.4 * 0.27)
    truth["x1"] = sn.meta["SIM_SALT2x1"]
    truth["t0"] = sn.meta["SIM_PEAKMJD"]
    truth["mwebv"] = sn.meta["SIM_MWEBV"]
    truth["z"] = sn.meta["REDSHIFT_FINAL"]
    if replacements is not None:
        addbands(sn, lsstbands, replacements)
    return sn, truth

Use the functions for a fit


In [9]:
lsstdeeplocation='/Users/rbiswas/data/SNDATA/SIM/PUBLIC_LSSTDEEP'

In [10]:
data = loadSNANAData(snanafileroot='LSST_Ia', location=lsstdeeplocation)

In [11]:
lsstbands = ['u', 'g', 'r', 'i', 'z', 'Y']
replacement = ['LSST_' + flt for flt in lsstbands]
print replacement, lsstbands


['LSST_u', 'LSST_g', 'LSST_r', 'LSST_i', 'LSST_z', 'LSST_Y'] ['u', 'g', 'r', 'i', 'z', 'Y']

In [12]:
source = sncosmo.get_source('salt2-extended')
dust = sncosmo.CCM89Dust()
model = sncosmo.Model(source=source,
                      effects=[dust, dust],
                      effect_names=['host', 'mw'],
                      effect_frames=['rest', 'obs'])

In [13]:
sn, truth = reformat_SNANASN(data[1], lsstbands, replacement)
snr = sn['flux']/sn['fluxerr']
plt.hist(snr[snr > 3.])


Out[13]:
(array([ 4.,  1.,  2.,  2.,  1.,  1.,  1.,  1.,  1.,  1.]),
 array([  3.07398057,   3.9444258 ,   4.81487103,   5.68531625,
          6.55576148,   7.42620671,   8.29665194,   9.16709716,
         10.03754239,  10.90798762,  11.77843285]),
 <a list of 10 Patch objects>)

In [14]:
model.set(**truth)

In [ ]:
fit_res, fit_model = sncosmo.fit_lc(sn, model, vparam_names=['t0', 'x0', 'x1', 'c'], 
                                    bounds={'c':(-0.3, 0.3), 'x1':(-3.0, 3.0)}, minsnr=3.0)

In [ ]:
nest_res, nest_model = sncosmo.nest_lc(sn, model, vparam_names=['t0', 'x0', 'x1', 'c'], 
                                    bounds={'c':(-0.3, 0.3), 'x1':(-3.0, 3.0)}, guess_amplitude_bound=True, minsnr=3.0, 
                                    verbose=True)


 iter=   257 logz=-288.258610

In [ ]:
fig = sncosmo.plot_lc(sn, model=(model, fit_model, nest_model), model_label=('SIM', 'Fit', 'nest'))

Other Functions


In [ ]:
def lc(SNID, filename, location="./"):
    from astropy.io import fits
    summaryfile = location + '/' +filename + "_HEAD.FITS"
    photofile = location + '/' + filename + "_PHOT.FITS"

    
    summary = Table(fits.open(summaryfile)[1].data)
    
    snpointers = summary[np.array(map(int,summary["SNID"]))==SNID]
    ptrobs_min = snpointers["PTROBS_MIN"][0]
    ptrobs_max = snpointers["PTROBS_MAX"][0]
        
    lc = Table(fits.open(photofile)[1].data[ptrobs_min: ptrobs_max])
    lc.rename_column("FLUXCAL", "flux")
    lc.rename_column("FLUXCALERR","fluxerr")
    lc.rename_column("ZEROPT","ZP")
    #lc = photofile
    return lc