This notebook shows how to load the output for Eric's survey simulator ztf_sim and generate Macronova lightcurves for it using an SED from Rosswog et al. (2016). (Check out the other notebooks for examples how to simulate other transients.)
Note: You need to download Eric's newest sample output here. The link was also included in Eric's email, so you will likely only need to change the path below.
Furthermore you'll require the dust map from Schlegel, Finkbeiner & Davis (1998) for full functionality. It can be found here.
In [1]:
import os
home_dir = os.environ.get('HOME')
# Please enter the filename of the ztf_sim output file you would like to use. The example first determines
# your home directory and then uses a relative path (useful if working on several machines with different usernames)
survey_file = os.path.join(home_dir, 'data/ZTF/test_schedule_v6.db')
# Please enter the path to where you have placed the Schlegel, Finkbeiner & Davis (1998) dust map files
# You can also set the environment variable SFD_DIR to this path (in that case the variable below should be None)
sfd98_dir = os.path.join(home_dir, 'data/sfd98')
In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import simsurvey
import sncosmo
from astropy.cosmology import Planck15
import simsurvey_tools as sst
In [3]:
# Load the ZTF CCD corners and filters
ccds = sst.load_ztf_ccds()
sst.load_ztf_filters()
In [4]:
# Load simulated survey from file (download from ftp://ftp.astro.caltech.edu/users/ebellm/one_year_sim_incomplete.db)
# Currently DES filters are used as proxies for ZTF filters
plan = simsurvey.SurveyPlan(load_opsim=survey_file, band_dict={'g': 'ztfg', 'r': 'ztfr', 'i': 'desi'}, ccds=ccds)
mjd_range = (plan.cadence['time'].min() - 30, plan.cadence['time'].max() + 30)
In [5]:
# To review the pointing schedule, you can use this table
plan.pointings
Out[5]:
In [6]:
# Load phase, wavelengths and flux from file
phase, wave, flux = sncosmo.read_griddata_ascii('data/macronova_sed_wind20.dat')
# Create a time series source
source = sncosmo.TimeSeriesSource(phase, wave, flux)
# Create the model that combines SED and propagation effects
dust = sncosmo.CCM89Dust()
model = sncosmo.Model(source=source,
effects=[dust],
effect_names=['host'],
effect_frames=['rest'])
In [7]:
def random_parameters(redshifts, model,
r_v=2., ebv_rate=0.11,
**kwargs):
cosmo = Planck15
# Amplitude
amp = []
for z in redshifts:
d_l = cosmo.luminosity_distance(z).value * 1e5
amp.append(d_l**-2)
return {
'amplitude': np.array(amp),
'hostr_v': r_v * np.ones(len(redshifts)),
'hostebv': np.random.exponential(ebv_rate, len(redshifts))
}
The transient generator combines model and ditribution, and randomly draws all parameters needed to simulate the lightcurves.
(Note that here we also set the volumetric rate as funtion of z. For macronovae, a good guess would be $5\cdot 10^{-7}~\textrm{Mpc}^{-3}~\textrm{yr}^{-1}$ but this would only results in a couple of observed macronovae. For this example we'll use a 100 times larger rate.)
In [8]:
transientprop = dict(lcmodel=model,
lcsimul_func=random_parameters)
tr = simsurvey.get_transient_generator((0.0, 0.05),
ratefunc=lambda z: 5e-5,
dec_range=(-30,90),
mjd_range=(mjd_range[0],
mjd_range[1]),
transientprop=transientprop,
sfd98_dir=sfd98_dir)
In [9]:
survey = simsurvey.SimulSurvey(generator=tr, plan=plan)
lcs = survey.get_lightcurves(
#progress_bar=True, notebook=True # If you get an error because of the progress_bar, delete this line.
)
In [10]:
len(lcs.lcs)
Out[10]:
In [11]:
lcs[0]
Out[11]:
In [ ]:
In [12]:
lcs.save('lcs_tutorial_mne.pkl')
In [13]:
lcs = simsurvey.LightcurveCollection(load='lcs_tutorial_mne.pkl')
You can inspect the lightcurves manually. This example should return the lightcurve with the most points with S/N > 5.
In [14]:
_ = sncosmo.plot_lc(lcs[0])
The two figures below show how early the MNe are detected and at what redshifts. The simulation input parameters of transients that were not detected are also kept, so can check completeness.
In [15]:
plt.hist(lcs.stats['p_det'], lw=2, histtype='step', range=(0,10), bins=20)
plt.xlabel('Detection phase (observer-frame)', fontsize='x-large')
_ = plt.ylabel(r'$N_{MNe}$', fontsize='x-large')
In [16]:
plt.hist(lcs.meta_full['z'], lw=1, histtype='step', range=(0,0.05), bins=20, label='all')
plt.hist(lcs.meta['z'], lw=2, histtype='step', range=(0,0.05), bins=20, label='detected')
plt.xlabel('Redshift', fontsize='x-large')
plt.ylabel(r'$N_{MNe}$', fontsize='x-large')
plt.xlim((0, 0.05))
plt.legend()
Out[16]:
In [ ]: