In [1]:
import os
In [2]:
import numpy as np
In [3]:
from collections import Iterable
In [4]:
%matplotlib inline
import matplotlib.pyplot as plt
# import seaborn as sns
# sns.set()
In [5]:
import sncosmo
In [6]:
from astropy.table import Column
In [7]:
import transientSources.transients as ts
example_data = os.path.join(os.path.split(ts.__file__)[0], 'example_data')
In [8]:
location = os.path.join(example_data, '2007uy/cfa_2007uy')
The representation of the astrophysical object in terms of the spectra provided is in TransientObject
In [9]:
help(ts.TransientObject)
In [10]:
mysource = ts.TransientObject.fromDataDir(name='Karpenka-2007uy', spectralDataDir=location)
mysource is an instance of TransientObject that keeps track of how we built it. The attribute of mysource
that is most interesting for our purpose is SNCosmoSource. This is the source which we will use to build
SNCosmo Light curves. This is built from the restframe 'mangled' spectra provided by Natasha, and their phases
which do not know where the 'peak' is. 2007uy is of Type Ib.
In [11]:
print mysource.SNCosmoSource
In [12]:
peak = mysource.phasePeak
print peak
Without being told what the peak phase is, TransientObject calculates the peak phase to be at the restFrame 'besselB' band peak. Note that in the above instance, the peak has been calculated to 9.47, and so the spectra which originally had a date of 0. now has a phase of -9.47.
However,one can also set this peak if they have prior knowledge
In [13]:
mysource.setphasePeak(3.)
In [14]:
print mysource.SNCosmoSource
Note that now, the 0 phase is mapped to -3. in recognition of the fact that peak is at phase 3.0 of the supplied days.
Now we will set it back (assuming we don't know what the peak really is
In [15]:
mysource.setphasePeak(peak)
In [16]:
print mysource.SNCosmoSource
We can examine the minimum and maximum rest frame phases and wavelength
In [17]:
mysource.SNCosmoSource.minphase(), mysource.SNCosmoSource.maxphase()
Out[17]:
In [18]:
mysource.SNCosmoSource.minwave(), mysource.SNCosmoSource.maxwave()
Out[18]:
This is one of the disadvantages related to single spectra (as opposed to ensembles) due to the lack of data (the rest of the notebook will show advantages of this approach in this regime of low data).
In [19]:
model = sncosmo.Model(source=mysource.SNCosmoSource)
In [20]:
print model
In [63]:
#Set parameters to match the photometry
# t0 is a guess provided for the peak of the observed photometry.
In [21]:
t0 =54478
model.set(z=0.0007, t0=54478.)
earlymodel = sncosmo.Model(source=mysource.SNCosmoSource)
earlymodel.set(t0=54477., z=0.0007)
In [22]:
waves = model._source._wave * (1.0 + 0.0007)
In [23]:
fluxes = model.flux(time=np.arange(-5.+ t0, 10. + t0), wave=waves)
In [24]:
fig = plt.figure()
ax1 = fig.add_subplot(221)
ax1.plot(waves, fluxes[0])
ax2 = fig.add_subplot(222)
for flux in fluxes:
ax2.plot(waves, flux)
The SEDs are single object, and are therefore noisy. It is possible to smooth them, if necessary. The right figure shows a bunch of the spectra. This is easier to see in animation.
In [25]:
filtlist = ['bessellB', 'bessellV', 'bessellr', 'besselli']
filtfig, ax = plt.subplots()
for filt in filtlist:
bp = sncosmo.get_bandpass(filt)
ax.plot(bp.wave, bp.trans)
ax.plot(waves, fluxes[0]*1e15)
The Bessell filters I am using and the spectra (multiplied by a constant to show on the figure). B, V are overlapped by the spectrum, i and r have support outside the wavlength range of spectrum, and therefore band fluxes cannot be computed.
We will compare the evolution of the spectral series with Ia (SALT2) and other templates for Ib /Ibc. Nugent Template is a template for Type Ibc built from an ensemble of objects. The other three templates are the Nugent template matched to the photometry of these SN: 2005hl, 2005hm and 2006jo. These were first referenced in Sako, 2011, and the code to do this was probably made by Saurabh Jha.
In [69]:
# saltsource = sncosmo.Model(source='salt2')._source
# nugent_1bc = sncosmo.Model(source='nugent-sn1bc')._source
# sn_2005_hl = sncosmo.Model(source='s11-2005hl')._source
# sn_2005_hm = sncosmo.Model(source='s11-2005hm')._source
# sn_2006_jo = sncosmo.Model(source='s11-2006jo')._source
In [70]:
# test = sncosmo.animate_source([mysource.SNCosmoSource, saltsource, nugent_1bc, sn_2005_hl, sn_2005_hm, sn_2006_jo],
# label=['Karpenka_2007_uy', 'salt2', 'NugentIbc', '2005_hl', '2005_hm', '2006_jo'],
# fname='Karpenka_comparison.mp4')
Given the current directory structure, we have a getPhotometry() method, that loads the photometry into
an astropy Table. This is illustrated below. Note that the fluxes are in units of ergs/cm^2/sec.
In [26]:
photometryDirectory = os.path.join(example_data, '2007uy/')
photometry = mysource.getPhotometry(photometryDirectory, filterList=['B', 'V', 'i', 'r'])
Here is a comparison of synthetic photometry photometry from SNCosmo's interpolated time series spectra, and the photometry. Note that SNCosmo does not calculate fluxes in ergs/cm^2/sec but in counts/cm^2/sec. So, the comparison is done with a branch of SNCosmo in https://github.com/rbiswas4/sncosmo/tree/CC_templates (27fbf6db25d1474ed1fb4cfbe7fd8ff680e9e4dc). So for the rest of this notebook to work you need to install this particular branch of SNCosmo.
There is a factor of 1000, that I have not managed to track down yet that is different between the photometry and the synthetic photometry flux values.
In [27]:
def renorm_bp(band):
b = sncosmo.get_bandpass(band)
norm = np.sum(b.trans * b.dwave)
bp = sncosmo.Bandpass(wave=b.wave, trans=b.trans/norm)
return bp
In [34]:
def _plotPhotometry(band, model=model):
photometry = mysource.getPhotometry(photometryDirectory, filterList=[band])
col = Column(name='ren', data=map(renorm_bp, np.asarray(photometry['band'])) )
photometry.add_column(col)
fig, ax = plt.subplots()
# bp = sncosmo.get_bandpass()
if isinstance(model, Iterable):
for m in model:
ax.plot(photometry['time'], m.bandflux(time=np.asarray(photometry['time']), band=np.asarray(photometry['ren'])))
else:
ax.plot(photometry['time'], model.bandflux(time=np.asarray(photometry['time']), band=np.asarray(photometry['ren'])), 'rs')
ax.errorbar(photometry['time'], photometry['flux'], photometry['fluxerr'], fmt='k.')
ax.grid(True)
return fig
In [35]:
figB = _plotPhotometry('B', model=[model, earlymodel])
In [ ]:
In [30]:
figV= _plotPhotometry('V', model=[model, earlymodel])
The difference in the two bands between the red (photometry files) and the blue (synthetic photometry) from the SNCosmo source probably shows the problems with interpolation of the spectra. In the coming days, we will check a few different things which we will list in a document and issues in this repository.
In [180]:
t0
Out[180]:
In [43]:
fig, ax = plt.subplots()
photometry = mysource.getPhotometry(photometryDirectory, filterList=['B'])
col = Column(name='ren', data=map(renorm_bp, np.asarray(photometry['band'])) )
photometry.add_column(col)
ax.plot(photometry['time'], model.bandflux(time=np.asarray(photometry['time']), band=np.asarray(photometry['ren'])))
# Code to plot the band flux at a single epoch of 54481.
ax.plot(54481., model.bandflux(time=54481., band=renorm_bp('BessellB')),'o')
ax.grid(True)
Let us look at the photometry comparison in B band a little more closely. Our first guess would be that the bivariate spline interpolation we used to go from the mangled spectra available at a discrete set of times (which we will refer to as the native grid of the model) to the spectra in between the times is the point of failure.
This should mean that the mismatch between the photometry points and the simulated points should be small if the photometry is at times close to the native grid points. the failure should also be bad beyond the minimum and maximum times corresponding to the native grid.
In [181]:
obsTimes = photometry[photometry['band'] == 'bessellB']['time']
In [182]:
nativeTimes = mysource.SNCosmoSource._phase * (1.0 + model.get('z')) + model.get('t0')
earlyNativeTimes = mysource.SNCosmoSource._phase * (1.0 + earlymodel.get('z')) + earlymodel.get('t0')
In [183]:
ax = figB.axes[0]
map(lambda x: ax.axvline(x, linestyle='dashed', lw =2.), nativeTimes)
Out[183]:
In [184]:
map(lambda x: ax.axvline(x, linestyle='dashed', color='g', lw =2.), earlyNativeTimes)
Out[184]:
In [185]:
figB.set_size_inches(12, 8)
In [186]:
ax = figB.axes[0]
ax.set_xlim(54460,54520 )
Out[186]:
In [187]:
figB
Out[187]:
This figure shows the photometry (black points with errors), and two curves with the same model except the peak position is one day apart. The vertical lines map the native phases of the model source to this plot.
So, this is quite strange we can see that there are points with appreciable difference with the photometry even though they are sampled at the same times. If the blue (t0 = 544 There are two possibilities that could explain this:
In [213]:
figV.set_size_inches(12,8)
In [214]:
ax = figV.axes[0]
map(lambda x: ax.axvline(x, linestyle='dashed', lw =2., color='b'), nativeTimes)
map(lambda x: ax.axvline(x, linestyle='dashed', lw =2., color='g'), earlyNativeTimes)
Out[214]:
The comparison seems better for V band.
In [215]:
figV
Out[215]:
In [188]:
from utils import plotutils as pu
In [203]:
def _comparetonative(mysource):
_phases = mysource.SNCosmoSource._phase
_wave = mysource.SNCosmoSource._wave
fig, ax0, ax1 = pu.settwopanel(setdifflimits=(0.05, 2.0))
for i, _phase in enumerate(_phases):
ax0.plot(_wave, mysource.SNCosmoSource.flux(phase=_phase, wave=_wave))
ax0.plot(_wave, mysource.mangled_spectra[i][:, 1], ls = 'dashed')
ax1.plot(_wave, mysource.SNCosmoSource.flux(phase=_phase, wave=_wave) / mysource.mangled_spectra[i][:, 1])
In [205]:
_comparetonative(mysource)