Notebook for a mocked up GRB afterglow lightcurve.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import healpy as hp
import lsst.sims.maf.slicers as slicers
import lsst.sims.maf.metrics as metrics
import lsst.sims.maf.metricBundles as metricBundles
import lsst.sims.maf.db as db
In [2]:
import matplotlib as mpl
mpl.rcParams['figure.autolayout'] = False
On-axis afterglows decay as $F \sim t^{-\alpha}$ until a jet break, with $\alpha \sim 1$.
This requires subclassing TransientMetric, which assumes linear rise & decays in magnitude. We therefore have writted GRBTransientMetric that synthesizes the appropriate lightcurve shape.
We will use the rough numbers in 2011PASP..123.1034J: a Gaussian with $\mu$=15.35 mag, $\sigma=1.59$ in R-band at 1 minute after explosion, $\alpha = 1.0$.
The right way to do this is not to assume fiducial values, but to simulate a cosmological population of events.
In [3]:
from mafContrib import GRBTransientMetric, TripletBandMetric
We'll calculate detection fractions for one, two, and three detections. For testing we may want to use one year at a time:
Now let's plot the lightcurves. Because we include a random draw from the brightness distribution these plots will change each time we re-execute the cell.
In [13]:
surveyDuration=10
if surveyDuration < 10:
year = 3
sqlconstraint = 'night between %f and %f' % ((365.25*year,365.25*(year+1)))
else:
sqlconstraint = ''
In [14]:
colors = ['b','g','r','purple','y','magenta','k']
filterNames = ['u','g','r','i','z','y']
transDuration = 100 # days--controls how frequently we inject new transients.
transMetric = GRBTransientMetric(alpha= 1.,surveyDuration=surveyDuration, transDuration=transDuration)
times = np.logspace(-3,2,500)
for filterName, color in zip(filterNames,colors):
filters = np.array([filterName]*times.size)
lc = transMetric.lightCurve(times % transDuration,filters)
plt.plot(times,lc, color, label=filterName)
plt.xlabel('time (days)')
plt.ylabel('mags')
plt.gca().invert_yaxis()
plt.xscale('log')
plt.legend()
Out[14]:
In [15]:
# Pick a slicer
slicer = slicers.HealpixSlicer(nside=16)
summaryMetrics = [metrics.MinMetric(), metrics.MeanMetric(), metrics.MaxMetric(),
metrics.MedianMetric(), metrics.RmsMetric(),
metrics.PercentileMetric(percentile=25), metrics.PercentileMetric(percentile=75)]
# Configure some metrics
metricList = []
# What fraction of transients are detected at least once?
metricList.append(transMetric)
# create metrics for 2 and 3 detections
m2 = GRBTransientMetric(alpha= 1., surveyDuration=surveyDuration, nPerFilter=2, metricName='GRBTransientMetric_k2')
#m2 = TripletBandMetric()
metricList.append(m2)
m3 = GRBTransientMetric(alpha= 1., surveyDuration=surveyDuration, nPerFilter=3, metricName='GRBTransientMetric_k3')
metricList.append(m3)
In [19]:
# Set the database and query
#runName = 'minion_1016' # baseline cadence
runName = 'enigma_1281' # NEO with triples
bDict={}
for i,metric in enumerate(metricList):
bDict[i] = metricBundles.MetricBundle(metric, slicer, sqlconstraint,
runName=runName, summaryMetrics=summaryMetrics)
In [20]:
opsdb = db.OpsimDatabase(runName + '_sqlite.db')
outDir = 'Transients'
resultsDb = db.ResultsDb(outDir=outDir)
bgroup = metricBundles.MetricBundleGroup(bDict, opsdb, outDir=outDir, resultsDb=resultsDb)
bgroup.runAll()
bgroup.plotAll(closefigs=False)
In [21]:
# enigma_1281
for key in sorted(bDict):
if bDict[key].metric.name.endswith('anyNfilters'):
bDict[key].computeSummaryStats(resultsDb=resultsDb)
print bDict[key].metric.name, bDict[key].summaryValues
In [18]:
# minion_1016
for key in sorted(bDict):
if bDict[key].metric.name.endswith('anyNfilters'):
bDict[key].computeSummaryStats(resultsDb=resultsDb)
print bDict[key].metric.name, bDict[key].summaryValues
So half of these afterglows are detected in one epoch, 11% in two, and 3% in three epochs (in any filter).
Let's try to understand these detection efficiencies by looking more generically at how often we have three observations in a single filter within a relevant interval to detect a GRB.
We'll start by determining how long LSST can observe afterglows parameterized as above--this is be band-dependent, as both LSST's limiting magnitude and the afterglow color vary. For simplicity we consider just $r$ band for now.
In [19]:
colors = ['b','g','r','purple','y','magenta','k']
filterNames = ['u','g','r','i','z','y']
transDuration = 300 # days--controls how frequently we inject new transients.
transMetric = GRBTransientMetric(alpha= 1.,surveyDuration=1, transDuration=transDuration)
def find_time_at_limiting_mag(limiting_mag,times,lc):
return np.interp(limiting_mag,lc,times)
limit_times = []
lc_mins = []
#times = np.linspace(0,transDuration+1,500)
times = np.logspace(-2,2,500)
for i in range(300):
# since we haven't included band-dependent afterglow brightness yet,
# look at the brightness distribution as a function of time in a single band
filters = np.array(['r']*times.size)
lc = transMetric.lightCurve(times,filters)
lc_mins.append(lc[-1])
limit_times.append(find_time_at_limiting_mag(24.7,times,lc))
# plt.plot(times,lc)
#plt.xlabel('time (days)')
#plt.ylabel('mags')
#plt.gca().invert_yaxis()
#plt.xscale('log')
plt.hist(limit_times)
plt.xlabel('Time after trigger last detectable by LSST (days)')
plt.ylabel('Number of simulated bursts')
print np.percentile(limit_times,50)
So there is a substantial tail of afterglows detectable to tens or even hundreds of days after the burst. Let's look at the cumulative histogram:
In [20]:
plt.hist(limit_times,normed=True,cumulative=True,bins=np.logspace(-1,2,25))
plt.xscale('log')
plt.xlabel('Time after trigger last detectable by LSST (days)')
plt.ylabel('Cumulative Fraction')
Out[20]:
It is important to remember that we have not included a jet break in our lightcurve model here--in reality we expect the afterglow decay to decay faster at late times. Still, it does appear that with appropriately-chosen cadences it should be possible to detect afterglows in multiple visits.