Sample release for GW190425¶
This notebook serves as a basic introduction to loading and viewing data released in associaton with the publication titled GW190425: Observation of a compact binary coalescence with total mass $\sim 3.4 M_{\odot}$ avaliable through DCC and arXiv.
The data used in these tutorials will be downloaded from the public DCC page LIGO-P2000026.
The released data file can be read in using the PESummary or h5py libraries. For this notebook we'll start with simple stuff using h5py. Then we'll use PESummary v0.3.0 to read the data files as well as for plotting. For general instructions on how to manipulate the data file and/or read this data file with h5py, see the PESummary docs.
# import useful python packages
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import h5py
Some simple stuff with "vanilla" h5py
# read in the data
fn = "GW190425_posterior_samples.h5"
data = h5py.File(fn,'r')
# print out top-level data structures
print('Top-level data structures:',data.keys())
# print out parametrized waveform family names ("approximants" in LIGO jargon).
# HS, LS = high-spin prior and low-spin prior, respectively
print('approximants:',data['approximant'].keys())
# extract posterior samples for one of the approximants
posterior_samples = data['posterior_samples']['PhenomDNRT-HS']
print('data structures in posterior_samples:',posterior_samples.keys())
pnames = [item.decode("utf-8") for item in posterior_samples['parameter_names']]
print('parameter names:',pnames)
# extract the samples data into an numpy array:
samples = np.array(posterior_samples['samples']).T
# get samples for one of the parameters
ind = pnames.index('luminosity_distance')
dL = samples[ind]
print('dL shape, mean, std =',dL.shape,dL.mean(),dL.std())
# smooth it
from scipy.stats.kde import gaussian_kde
hs = gaussian_kde(dL)
# histogram, and overlay the smoothed PDF
plt.figure()
h, b, o = plt.hist(dL,bins=100)
hsmoothed = hs(b)*len(dL)*(b[1]-b[0])
plt.plot(b,hsmoothed)
plt.xlabel('luminosity distance')
plt.ylabel('posterior PDF')
plt.show()
# release memory for the data
del data
Now use PESummary v0.3.0 to read the data files as well as for plotting.
# import ligo-specific python packages.
# pesummary is a ligo-specific python package for reading and plotting the results of Bayesian parameter estimation.
# Install with "pip install pesummary" , and make sure you have version >= 0.3.0.
import pesummary
from pesummary.gw.file.read import read
print(pesummary.__version__)
There are 6 different approximants that were used to analyze GW190425 and they are all stored in the data file.
fn = "GW190425_posterior_samples.h5"
data = read(fn)
labels = data.labels
print(labels)
To illustrate the data structure we'll pick one approximant by random and plot its respective data.
samples_dict = data.samples_dict
posterior_samples = samples_dict["PhenomPNRT-HS"]
prior_samples = data.priors["samples"]["PhenomPNRT-HS"]
parameters = posterior_samples.keys()
print(parameters)
As an example, we'll show the different posterior distributions derived for a single waveform and the posterior distribution derived using the different approximants for the luminosity_distance parameter.
from pesummary.core.plots.plot import _1d_histogram_plot
from pesummary.gw.plots.latex_labels import GWlatex_labels
parameter = "luminosity_distance"
latex_label = GWlatex_labels[parameter]
fig = _1d_histogram_plot(
parameter, posterior_samples[parameter], latex_label, prior=prior_samples[parameter]
)
fig.set_size_inches(12, 8)
plt.show()
from pesummary.core.plots.plot import _1d_comparison_histogram_plot
samples = []
for label in labels:
samples.append(samples_dict[label][parameter])
colors = ['b', 'r', 'k', 'y', 'orange', 'g']
fig = _1d_comparison_histogram_plot(parameter, samples, colors, latex_label, labels, kde=True)
fig.set_size_inches(12, 8)
plt.show()
Make a corner plot:
from pesummary.gw.plots.plot import _make_corner_plot
fig = _make_corner_plot(posterior_samples, GWlatex_labels)
plt.show()
#plt.savefig(fn+'_corner.png')
The psds that were used for each analysis can also be extracted from this file and plotted
from pesummary.gw.plots.plot import _psd_plot
psd = data.psd["PhenomPNRT-HS"]
ifos = list(psd.keys())
frequencies, strains = [], []
for ifo in ifos:
frequencies.append(np.array(psd[ifo]).T[0])
strains.append(np.array(psd[ifo]).T[1])
fig = _psd_plot(frequencies, strains, labels=ifos, fmin=19.4)
fig.set_size_inches(12, 8)
plt.show()
The calibration envelopes that were used in this analysis can also be extracted from this file and plotted
from pesummary.gw.plots.plot import _calibration_envelope_plot
prior = data.priors["calibration"]["PhenomPNRT-HS"]
calibration = data.calibration["PhenomPNRT-HS"]
frequencies = np.arange(20., 1024., 1. / 4)
calibration_data, prior_data = [], []
for ifo in ifos:
calibration_data.append(np.array(calibration[ifo]))
prior_data.append(np.array(prior[ifo]))
fig = _calibration_envelope_plot(frequencies, calibration_data, ifos, prior=prior_data)
fig.set_size_inches(16.5, 10.5)
plt.show()
The configuration file that were used for each analysis can also be extracted from this file
config = data.config["PhenomPNRT-HS"]
for i in config.keys():
print("[{}]".format(i))
for key, item in config[i].items():
print("{}={}".format(key, item[0].decode("utf-8")))
print("\n")