In this notebook, we'll do a crude estimate of the mass of the binary that produced GW170817. We'll make the assumption in this notebook that the component objects have equal mass and aren't spinning.
import sys
!{sys.executable} -m pip install pycbc ligo-common --no-cache-dir
Requirement already satisfied: pycbc in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages/PyCBC-9rc6331-py3.7-linux-x86_64.egg (9rc6331) Requirement already satisfied: lalsuite in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (6.70) Requirement already satisfied: ligo-common in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (1.0.3) Requirement already satisfied: numpy>=1.16.0 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (1.19.0) Requirement already satisfied: Mako>=1.0.1 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (1.1.3) Requirement already satisfied: cython>=0.29 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (0.29.20) Requirement already satisfied: decorator>=3.4.2 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (4.4.2) Requirement already satisfied: matplotlib>=1.5.1 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (3.2.2) Requirement already satisfied: pillow in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (7.1.2) Requirement already satisfied: h5py>=2.5 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (2.10.0) Requirement already satisfied: jinja2 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (2.11.2) Requirement already satisfied: mpld3>=0.3 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (0.5.1) Requirement already satisfied: lscsoft-glue>=1.59.3 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (2.0.0) Requirement already satisfied: emcee==2.2.1 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (2.2.1) Requirement already satisfied: requests>=1.2.1 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (2.24.0) Requirement already satisfied: beautifulsoup4>=4.6.0 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (4.9.1) Requirement already satisfied: six>=1.10.0 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (1.15.0) Requirement already satisfied: ligo-segments in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (1.2.0) Requirement already satisfied: tqdm in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (4.46.1) Requirement already satisfied: gwdatafind in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (1.0.4) Requirement already satisfied: astropy>=2.0.3 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (4.0.1.post1) Requirement already satisfied: scipy>=0.16.0 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pycbc) (1.5.0) Requirement already satisfied: python-dateutil in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from lalsuite) (2.8.1) Requirement already satisfied: MarkupSafe>=0.9.2 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from Mako>=1.0.1->pycbc) (1.1.1) Requirement already satisfied: cycler>=0.10 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from matplotlib>=1.5.1->pycbc) (0.10.0) Requirement already satisfied: kiwisolver>=1.0.1 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from matplotlib>=1.5.1->pycbc) (1.2.0) Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from matplotlib>=1.5.1->pycbc) (2.4.7) Requirement already satisfied: pyOpenSSL in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from lscsoft-glue>=1.59.3->pycbc) (19.1.0) Requirement already satisfied: idna<3,>=2.5 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from requests>=1.2.1->pycbc) (2.9) Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from requests>=1.2.1->pycbc) (1.25.9) Requirement already satisfied: certifi>=2017.4.17 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from requests>=1.2.1->pycbc) (2020.6.20) Requirement already satisfied: chardet<4,>=3.0.2 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from requests>=1.2.1->pycbc) (3.0.4) Requirement already satisfied: soupsieve>1.2 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from beautifulsoup4>=4.6.0->pycbc) (2.0.1) Requirement already satisfied: cryptography>=2.8 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (2.9.2) Requirement already satisfied: cffi!=1.11.3,>=1.8 in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from cryptography>=2.8->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (1.14.0) Requirement already satisfied: pycparser in /home/ahnitz/projects/PyCBC-Tutorials/env/lib/python3.7/site-packages (from cffi!=1.11.3,>=1.8->cryptography>=2.8->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (2.20)
!wget -nc https://dcc.ligo.org/public/0146/P1700349/001/H-H1_LOSC_CLN_4_V1-1187007040-2048.gwf
!wget -nc https://dcc.ligo.org/public/0146/P1700349/001/L-L1_LOSC_CLN_4_V1-1187007040-2048.gwf
--2020-06-27 14:35:13-- https://dcc.ligo.org/public/0146/P1700349/001/H-H1_LOSC_CLN_4_V1-1187007040-2048.gwf Resolving dcc.ligo.org (dcc.ligo.org)... 131.215.125.144 Connecting to dcc.ligo.org (dcc.ligo.org)|131.215.125.144|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 57824232 (55M) Saving to: ‘H-H1_LOSC_CLN_4_V1-1187007040-2048.gwf’ H-H1_LOSC_CLN_4_V1- 100%[===================>] 55.14M 1.01MB/s in 52s 2020-06-27 14:36:06 (1.06 MB/s) - ‘H-H1_LOSC_CLN_4_V1-1187007040-2048.gwf’ saved [57824232/57824232] --2020-06-27 14:36:06-- https://dcc.ligo.org/public/0146/P1700349/001/L-L1_LOSC_CLN_4_V1-1187007040-2048.gwf Resolving dcc.ligo.org (dcc.ligo.org)... 131.215.125.144 Connecting to dcc.ligo.org (dcc.ligo.org)|131.215.125.144|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 62070438 (59M) Saving to: ‘L-L1_LOSC_CLN_4_V1-1187007040-2048.gwf’ L-L1_LOSC_CLN_4_V1- 100%[===================>] 59.19M 3.76MB/s in 36s 2020-06-27 14:36:44 (1.65 MB/s) - ‘L-L1_LOSC_CLN_4_V1-1187007040-2048.gwf’ saved [62070438/62070438]
%matplotlib inline
import pylab
from pycbc.filter import highpass
from pycbc.catalog import Merger
from pycbc.frame import read_frame
merger = Merger("GW170817")
strain, stilde = {}, {}
for ifo in ['L1', 'H1']:
# We'll download the data and select 256 seconds that includes the event time
ts = read_frame("{}-{}_LOSC_CLN_4_V1-1187007040-2048.gwf".format(ifo[0], ifo),
'{}:LOSC-STRAIN'.format(ifo),
start_time=merger.time - 224,
end_time=merger.time + 32,
check_integrity=False)
# Read the detector data and remove low frequency content
strain[ifo] = highpass(ts, 15)
# Remove time corrupted by the high pass filter
strain[ifo] = strain[ifo].crop(4, 4)
# Also create a frequency domain version of the data
stilde[ifo] = strain[ifo].to_frequencyseries()
from pycbc.psd import interpolate, inverse_spectrum_truncation
psds = {}
for ifo in ['L1', 'H1']:
# Calculate a psd from the data. We'll use 2s segments in a median - welch style estimate
# We then interpolate the PSD to the desired frequency step.
psds[ifo] = interpolate(strain[ifo].psd(2), stilde[ifo].delta_f)
# We explicitly control how much data will be corrupted by overwhitening the data later on
# In this case we choose 2 seconds.
psds[ifo] = inverse_spectrum_truncation(psds[ifo], int(2 * strain[ifo].sample_rate),
low_frequency_cutoff=15.0,
trunc_method='hann')
pylab.loglog(psds[ifo].sample_frequencies, psds[ifo], label=ifo)
pylab.xlim(20, 1024)
pylab.ylim(1e-47, 1e-42)
pylab.legend()
<matplotlib.legend.Legend at 0x7f796d182f50>
from pycbc.waveform import get_fd_waveform
from pycbc.filter import matched_filter
from pycbc.conversions import mass1_from_mchirp_q
import numpy
# We will try different component masses and see which gives us the largest
chirp_mass = numpy.arange(1.1966, 1.1996, .0001)
masses = mass1_from_mchirp_q(chirp_mass, 1)
# Variables to store when we've found the max
hmax, smax, tmax, mmax, nsnr = None, {}, {}, 0, 0
snrs = []
for m in masses:
#Generate a waveform with a given component mass; assumed equal mass, nonspinning
hp, hc = get_fd_waveform(approximant="TaylorF2",
mass1=m, mass2=m,
f_lower=20, delta_f=stilde[ifo].delta_f)
hp.resize(len(stilde[ifo]))
# Matched filter the data and find the peak
max_snr, max_time = {}, {}
for ifo in ['L1', 'H1']:
snr = matched_filter(hp, stilde[ifo], psd=psds[ifo], low_frequency_cutoff=20.0)
# The complex SNR at the peak
snr = snr.time_slice(merger.time - 1, merger.time + 1)
_, idx = snr.abs_max_loc()
max_snr[ifo] = snr[idx]
# The time of the peak
max_time[ifo] = float(idx) / snr.sample_rate + snr.start_time
network_snr = (abs(numpy.array(list(max_snr.values()))) ** 2.0).sum() ** 0.5
snrs.append(max_snr)
# Keep track of only the loudest peak
if network_snr > nsnr:
tmax, hmax, mmax, smax = max_time, hp, m, max_snr
nsnr = network_snr
# See the SNR as a function of the component mass. Notice where this peaks as it gives us
# an estimate of what the parameters of the source system are. Note that masses
# here are in the *detector* frame, so if the source is located far away, it will in
# fact correspond to a lighter system due to cosmological redshift.
print("We found the best Mass1=Mass2 was %2.2f solar masses (detector frame)" % mmax)
We found the best Mass1=Mass2 was 1.38 solar masses (detector frame)
# We can see how the estimate from each IFO agrees
import numpy
network_snrs = [abs(numpy.array(list(s.values()))**2.0).sum()**0.5 for s in snrs]
for ifo in ['L1', 'H1']:
ind_snrs = [abs(s[ifo]) for s in snrs]
pylab.plot(masses, ind_snrs, label=ifo)
pylab.plot(masses, network_snrs, label='Network')
pylab.ylabel('Signal-to-noise')
pylab.xlabel('Mass of each neutron star (Detector frame Solar Masses)')
pylab.legend()
<matplotlib.legend.Legend at 0x7f796ca59c10>
for ifo in ['L1', 'H1']:
# Whiten the data
hoft = (stilde[ifo] / psds[ifo] ** 0.5).to_timeseries()
# Select the time around the event
zoom = hoft.time_slice(merger.time - 30, merger.time + 2)
# Calculate the qtransform (a kind of time-frequency representation similar to a spectrogram)
times, freqs, power = zoom.qtransform(.01, logfsteps=100, frange=(20, 512), qrange=(110, 110))
pylab.figure(figsize=(15,3))
pylab.pcolormesh(times - merger.time, freqs, power**0.5, vmin=0, vmax=6)
pylab.ylim(20, 512)
pylab.title('Interferometer: %s' % ifo)
pylab.xlabel('Time (s)')
pylab.ylabel('Frequency (Hz)')
pylab.yscale('log')
from pycbc.filter import sigma
for ifo in ['L1', 'H1']:
# Shift the template to the maximum time at this sample rate
dt = tmax[ifo] - stilde[ifo].start_time
inj = hmax.cyclic_time_shift(dt)
# Scale the template to the SNR and phase we measured above
inj /= sigma(hmax, psd=psds[ifo], low_frequency_cutoff=20.0) # This scales the template to unit SNR
inj *= smax[ifo] # This scales the template to the SNR / phase we found
# Subtract from the data
stilde2 = stilde[ifo] - inj
# Whiten the data
hoft = (stilde2 / psds[ifo] ** 0.5).to_timeseries()
# Select the time around the event
zoom = hoft.time_slice(merger.time - 30, merger.time + 2)
# Calculate the qtransform (a kind of time-frequency representation similar to a spectrogram)
times, freqs, power = zoom.qtransform(.01, logfsteps=100, frange=(20, 512), qrange=(110, 110))
pylab.figure(figsize=(15,3))
pylab.pcolormesh(times - merger.time, freqs, power**0.5, vmin=0, vmax=6)
pylab.ylim(20, 512)
pylab.title('Interferometer: %s' % ifo)
pylab.xlabel('Time (s)')
pylab.ylabel('Frequency (Hz)')
pylab.yscale('log')
# Note that a better match to the data can be found if you allow more
# freedom in the possible source parameters than we
# have here.