This notebook demonstrates the Python library rtpipe, which was used to make the first interferometric detection of a Fast Radio Burst on this data set. The results are presented in Chatterjee et al 2017, "The direct localization of a fast radio burst and its host" (http://dx.doi.org/10.1038/nature20797).
Prerequisites to run this analysis:
pip install rtpipeThe original analysis was run in quasi-real-time on a cluster at the Very Large Array using realfast. This demonstration is tuned to run on a laptop, find one burst from FRB 121102, and perform a rough analysis. For a more exact reproduction of the results in the publication, see the other notebooks in https://github.com/caseyjlaw/FRB121102.
This analysis can be used to reproduce the original detection or refine that analysis. You can also download all 9 data sets to create visualizations of all bursts, such as those at https://youtu.be/i3x0sBCQ_c8 or the sound files at http://dx.doi.org/10.7910/DVN/QSWJE6.
In [22]:
%matplotlib inline
In [23]:
from __future__ import print_function
import pylab as pl
import numpy as np
import rtpipe
import rtlib_cython as rtlib
In [24]:
# useful function for data visualization in notebook
def correctdata(st, data, u=None, v=None, w=None, corr='ph,dm', lm = None):
data2 = data.copy()
if 'ph' in corr:
if lm and len(u) and len(v):
l1, m1 = lm
rtlib.phaseshift_threaded(data2, st, l1, m1, u, v)
else:
print('No phase center provided')
if 'dm' in corr:
rtlib.dedisperse_par(data2, st['freq'], st['inttime'], st['dmarr'][0], [0, st['nbl']])
return data2
1) Go to dataverse, download the tarred data file for the first burst (detected 23 August 2016) and appropriate calibration (*.GN) file. E.g.,
- 16A-459_TEST_1hr.57623.72670021991.cut.tar
- 16A-459_TEST_1hr.57623.72670021991.GN
2) Untar data (e.g., 'tar xvf 16A-459_TEST_1hr.57623.72670021991.cut.tar')
3) Move to directory with sdmfile and gainfile
In [25]:
sdmfile = '16A-459_TEST_1hr.57623.72670021991.cut'
gainfile = '16A-459_TEST_1hr.57623.72670021991.GN'
# a segment defines a piece of data in time. some tuning required for other FRB 121102 data files.
segment = 0
In [26]:
sdm = rtpipe.parsesdm.getsdm(sdmfile)
bdfscan = 6
The rtpipe library has a few different interfaces. First, we show the end-to-end, real-time analysis that uses Python's multiprocessing library with shared memory. After that, we demonstrate another set of functions to run all the steps of the analysis and visualize the results.
All analysis begins by defining a dictionary with metadata from the observation, preferences for the search pipeline, and parameters derived from the joint consideration of both. The command below is tuned to find the known FRB with modest computational expense.
In [27]:
st = rtpipe.RT.set_pipeline(sdmfile, bdfscan, dmarr=[567.], dtarr=[1],
flaglist=[('badap', 3., 0.2)], uvoversample=1.5,
gainfile=gainfile, flagantsol=True, timesub='mean',
npix_max=1024, logfile=False, savecands=False, savenoise=False)
In [28]:
rtpipe.RT.pipeline(st, segment)
The previous command should find the first burst from FRB 121102 with SNR=29. This data was recorded on 23 August 2016 and the burst was found 5 hours later by realfast.
Rather than using this pipelined approach, we can use another interface to get calibrated, flagged data for playing.
In [29]:
data = rtpipe.RT.pipeline_reproduce(st, candloc=[segment,0,0,0,0], product='data')
print('Data has shape of (integrations, baselines, channels, polarizations): {}'.format(data.shape))
In [30]:
u, v, w = rtpipe.parsesdm.get_uvw_segment(st, segment)
burst_int = 400 # the burst starts in this integration at the top of the band
datadm = correctdata(st, data, corr='dm')
im = rtpipe.RT.sample_image(st, datadm, u, v, w, i=burst_int)
l,m = rtpipe.RT.calc_lm(st, im)
ra_center, dec_center = st['radec']
ra, dec = rtpipe.reproduce.source_location(ra_center, dec_center, l, m)
print('\nDirty image:')
print('Peak significance: {0:.1f}'.format(im.max()/im.std()))
print('Peak location (J2000): {0}, {1}'.format(ra, dec))
In [31]:
# show me!
fig = pl.figure(figsize=(12,5))
# full field
ax = fig.add_subplot(121)
pl.imshow(im.transpose(), interpolation='nearest', cmap='magma_r')
pl.xlabel('RA offset (pix)', fontsize=20)
pl.ylabel('Dec offset (pix)', fontsize=20)
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=3, color='w')
ax.yaxis.set_tick_params(width=3, color='w')
# zoom
ax = fig.add_subplot(122)
zoom = (450, 650, 550, 350)
pl.imshow(im.transpose()[zoom[3]:zoom[2], zoom[0]:zoom[1]], interpolation='nearest',
cmap='magma_r', extent=zoom)
pl.xlabel('RA offset (pix)', fontsize=20)
pl.ylabel('Dec offset (pix)', fontsize=20)
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=3, color='w')
ax.yaxis.set_tick_params(width=3, color='w')
In [32]:
# info to phase data
(l1, m1) = (-3.85e-04, 5.45e-04)
# phase up
dataph = correctdata(st, data, u, v, w, corr='ph', lm = (-3.85e-04,5.45e-04))
sgram = dataph[burst_int-50:burst_int+100].mean(axis=3).mean(axis=1).real.transpose()
# phase up and dedisperse
dataphdm = correctdata(st, data, u, v, w, corr='ph,dm', lm = (-3.85e-04,5.45e-04))
sgramd = dataphdm[burst_int-50:burst_int+100].mean(axis=3).mean(axis=1).real.transpose()
In [33]:
fig = pl.figure(figsize=(12,5))
# dynamic spectrum without dedispersion
ax = fig.add_subplot(121)
pl.imshow(sgram, interpolation='nearest', origin='lower', cmap='viridis',
extent=(0, len(sgram)*0.005, 2.5, 3.5), vmax=0.8*sgram.max())
pl.xlabel('Time (seconds)', fontsize=20)
pl.ylabel('Frequency (GHz)', fontsize=20)
ax.set_title('No dedispersion', fontsize=20)
# dynamic spectrum corrected for dispersion
ax = fig.add_subplot(122)
pl.imshow(sgramd, interpolation='nearest', origin='lower', cmap='viridis',
extent=(0, len(sgram)*0.005, 2.5, 3.5), vmax=0.8*sgram.max())
pl.xlabel('Time (seconds)', fontsize=20)
pl.ylabel('Frequency (GHz)', fontsize=20)
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=3, color='w')
ax.yaxis.set_tick_params(width=3, color='w')
ax.set_title('Dedispersed', fontsize=20)
Out[33]:
In [34]:
spectra = dataphdm[burst_int].mean(axis=0).real
fig = pl.figure(figsize=(15,8))
ax = fig.add_subplot(111)
pl.plot(spectra[:,0], 'c.', label='RR')
pl.plot(spectra[:,1], 'm.', label='LL')
pl.xlabel('Channel', fontsize=20)
pl.ylabel('Amplitude (system units)', fontsize=20)
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=3, color='w')
ax.yaxis.set_tick_params(width=3, color='w')
pl.legend(fontsize=14)
Out[34]:
In total, the VLA found 9 bursts from FRB 121102 from August through September 2016. Eight other data sets (visibility cutout plus calibration file) can be downloaded from the dataverse at https://doi.org/10.7910/DVN/TLDKXG.
You can also play with the next generation search pipeline at http://github.com/realfastvla/rfpipe.