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()
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
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
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]:
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)
In [ ]:
fig = sncosmo.plot_lc(sn, model=(model, fit_model, nest_model), model_label=('SIM', 'Fit', 'nest'))
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