Plotting gravitational wave strain data

When LIGO/Virgo announces a gravitational wave, both the samples and gravitational wave data are released to the public. Through pesummary’s pesummary.gw.file.strain module, we can not only plot the gravitational wave data, but we can also compare it to the maximum likelihood waveform from the parameter estimation analysis. Below we show an example for GW150914,

from pesummary.gw.file.strain import StrainData
from pesummary.io import read
import subprocess

# First we download the GW150914 posterior samples. We use the posterior
# samples generated by the `bilby` python package
subprocess.check_call(
    "curl https://dcc.ligo.org/public/0168/P2000193/002/pesummary_samples.zip -o pesummary_samples.zip", shell=True
)
subprocess.check_call("unzip pesummary_samples.zip", shell=True)

# Next we fetch the LIGO Livingston data around the time of GW150914
L1_data = StrainData.fetch_open_data('L1', 1126259446, 1126259478)

# We then read in the posterior samples and generate the maximum likelihood
# waveform in the time domain
samples = read("pesummary_samples/GW150914_bilby_pesummary.dat").samples_dict
maxL = samples.maxL_td_waveform("IMRPhenomPv2", 1. / 2048, 20., project="L1")

# Next we plot the data
fig = L1_data.plot(
    type="td", merger_time=1126259462.4, window=(-0.1, 0.1),
    template={"L1": maxL}, bandpass_frequencies=[40., 300.]
)
fig.savefig("GW150914.png")
fig.close()

and GW190814,

from pesummary.gw.file.strain import StrainData
from pesummary.gw.fetch import fetch_open_samples

# First we download the GW190814 posterior samples and generate the maximum
# likelihood waveform in the time domain
data = fetch_open_samples("GW190814", unpack=True, path="GW190814.h5")
samples = data.samples_dict["C01:SEOBNRv4PHM"]
maxL = samples.maxL_td_waveform("SEOBNRv4PHM", 1. / 4096, 20., project="L1")

# # Next we fetch the LIGO Livingston data around the time of GW190814
L1_data = StrainData.fetch_open_data('L1', 1249852257.01 - 20, 1249852257.01 + 5)

# Next we plot the data
fig = L1_data.plot(
    type="td", merger_time=1249852257.01, window=(-0.1, 0.04),
    template={"L1": maxL}, bandpass_frequencies=[50., 300.]
)
fig.savefig("GW190814.png")
fig.close()