Astrophysical variable sources including transients like supernovae are sources that are (roughly) fixed in position but change appreciably at human time scales, such that their measurements at different points of time in a survey are expected to be different. The key science of such sources has to do with interpreting such data.
In a photometric or imaging survey, the raw data associated with such measurements is in the form of a time series of images, along with metadata. This is often extracted into a catalog summarizing the inferred brightness of the source in the filter band and the uncertainties of the inferred. The uncertainties in this process are dominated by the sky noise which is Poisson distributed to a good approximation. Hence, the format of this recorded data is usually of the form
timestamp, flux, flux uncertainty, zero point, magnitude system
Note that the flux is really the band flux obtained in a filter band
The zero point, and magnitude system are methods of specifying the units of the reported flux, and flux uncertainty. flux may be reported so that the following equation with bandflux and std_bandflux in the same physical units holds. The std_bandflux depends on the transmissions and the magnitude system used.
bandmag := -2.5 log10(bandflux / std_bandflux) = -2.5 log10(flux) + zero point
so that a zero point of zero, implies that flux is bandflux in units of std_bandflux
the LightCurve class is meant to be a standard internal representation of this kind of catalog data in our software. This means that
In [56]:
import sncosmo
import analyzeSN as ans
import numpy as np
In [2]:
from analyzeSN import LightCurve
In [50]:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
The class is instantiated using a constructor that has a mandatory input of a dataframe. and optional inputs. The dataframe must contain the light curve information, and the minimum set of columns required are:
mjd, band, flux, flux_err, zp, zpsys or recognized aliases. Here, we load some data and demonstrate this
In [3]:
ex_data = sncosmo.load_example_data().to_pandas()
ex_data.head()
Out[3]:
This was a SNCosmo example dataset, loaded into pandas.DataFrame. Note that the column representing temporal information is time. Now, we will use this to instantiate our class.
In [4]:
lc = LightCurve(ex_data)
We get the representation in our convention by accessing the attribute lightCurve, which is a dataFrame with columns of the standard name. It could have more columns, but the minimal columns have to be there. If we supplied a dataFrame without having the minimal number of columns, it will raise a valueError, as described later. From the lightCurve attribute, we note that time has changed to our standard name mjd
In [5]:
lc.lightCurve.head()
Out[5]:
What would have happened if an essential column was missing? to check let us do the following:
In [6]:
lc_tmp = ex_data.copy()
del lc_tmp['band']
lc_tmp.head()
Out[6]:
In [7]:
lc_bad = LightCurve(lc_tmp)
This raises a ValueError pointing out that the column band was missing. However, it correctly recognizes that while
mjd was missing, it had the aliased column time instead and did not report it missing.
SNCosmo has a data format using astropy.table.Table. Since, we use SNCosmo to plot, or fit light curves, it is essential that we can regenerate the SNCosmo format data easily. We do this by the following method
In [8]:
SNCosmoLC = lc.snCosmoLC()
In [9]:
SNCosmoLC
Out[9]:
We can now plot the data using the usual SNCosmo plotting method, showing the data in g, r, and z bands.
In [10]:
sncosmo.plot_lc(SNCosmoLC)
Out[10]:
In [11]:
from analyzeSN import SNANASims
In [12]:
from astropy.table import Table
In [13]:
megacamBandNames = 'ugriz'
megacamRegisteredNames = tuple('mega_' + band for band in megacamBandNames)
In [14]:
snana_eg = SNANASims.fromSNANAfileroot(snanafileroot='snana_fits',
location='../analyzeSN/example_data/',
coerce_inds2int=False,
SNANABandNames=megacamBandNames,
registeredBandNames=megacamRegisteredNames)
In [15]:
Table.read(snana_eg.photFile)[45:50]
Out[15]:
Let us look at what the array in the fits file looks like:
In [16]:
snana_eg.bandNames
Out[16]:
In [17]:
snana_eg.newbandNames
Out[17]:
In [18]:
lcInstance = snana_eg.get_SNANA_photometry(snid='03D1aw')
In [19]:
lcInstance.lightCurve[['mjd', 'band', 'flux', 'fluxerr', 'zp', 'zpsys']]
Out[19]:
In [20]:
lcInstance.bandNameDict
Out[20]:
In [21]:
snana_eg.headData['REDSHIFT_FINAL']
Out[21]:
In [22]:
import sncosmo
To infer the light curve model parameters, we need a model. We will use the SALT model from SNCosmo
In [35]:
dust = sncosmo.CCM89Dust()
model = sncosmo.Model(source='salt2', effects=[dust],
effect_names=['host', 'mw'], effect_frames=['rest', 'obs'])
We will not try to infer the redshift first
In [36]:
model.set(z=0.5817)
We have two main methods of trying to infer the model parameters: the first one is a maximum likelihood estimate
In [37]:
resfit = sncosmo.fit_lc(lcInstance.snCosmoLC(), model, vparam_names=['t0', 'x0', 'x1', 'c'], modelcov=False)
In [38]:
res = sncosmo.mcmc_lc(lcInstance.snCosmoLC(), model, vparam_names=['t0', 'x0', 'x1', 'c'], modelcov=False)
In [39]:
from analyzeSN import ResChar
In [40]:
reschar = ResChar.fromSNCosmoRes(res)
In [33]:
reschar.salt_samples()['t0'].mean()
Out[33]:
In [42]:
reschar.salt_samples()['c'].mean()
Out[42]:
In [44]:
resfit[0].parameters[2]
Out[44]:
In [43]:
resfit[0].parameters[2]
Out[43]:
In [49]:
reschar.salt_samples()['c'].std()
Out[49]:
In [70]:
fig, ax = plt.subplots(1, 4)
reschar.salt_samples()['c'].hist(bins=np.arange(-0.6, 0.6, 0.03), lw=2., histtype='step', normed=1, ax=ax[0])
reschar.salt_samples()['x1'].hist(bins=np.arange(-4., 4., 0.2), lw=2., histtype='step', normed=1, ax=ax[1])
reschar.salt_samples()['t0'].hist(bins=10, lw=2., histtype='step', normed=1, ax=ax[2])
reschar.salt_samples()['mu'].hist(bins=10, lw=2., histtype='step', normed=1, ax=ax[3])
Out[70]:
In [80]:
sns.distplot(reschar.salt_samples()['t0'], rug=False, color='r', hist_kws=dict(histtype='step', lw=2., alpha=1., color='k'))
Out[80]:
In [47]:
resfit[0].errors['c']
Out[47]:
In [ ]:
10**(-0.4 * 0.27)
In [ ]:
-2.5 * np.log10(f)
In [ ]:
10.0**(0.27 / 2.5 )
In [ ]: