Sample release for GW190521

This notebook serves as a basic introduction to loading and viewing data released in association with the publications:
GW190521: A Binary Black Hole Merger with a Total Mass of 150 Msun avaliable through PRL, arXiv, and LIGO-P2000020.
Properties and astrophysical implications of the 150 Msun binary black hole merger GW190521 avaliable through ApjL, arXiv, and LIGO-P2000021.

The data used in this notebook must be downloaded from the public LIGO DCC page LIGO-P2000158.

Data release for GW190521

In public DCC page LIGO-P2000158, you will find several data files:

  • GW190521_posterior_samples.h5 - contains parameter estimation (PE) samples from many different runs:
    • NRSur7dq4 is our "preferred" waveform model (referred to as NRSur PHM in the papers). These samples contain the info for Table 1 and Figures 2, 3, 4, 5 of LIGO-P2000020.
    • IMRPhenomPv3HM and SEOBNRv4PHM are alternative waveform models for comparison. Those samples, along with those from NRSur7dq4, contain the info for Table 1 and Figures 1, 2, 3, 4, 5, 6, 8, 9 of LIGO-P2000021.
  • GW190521_studies_posterior_samples.h5 - contains parameter estimation (PE) samples from runs used for specialized studies; these should not be used in preference to the previous runs for most applications.
    • The three sets of PE samples for NRSur7dq4_Lmax2, Lmax3, and Lmax4 are used for Figure 7 of LIGO-P2000021.
    • The set of PE samples directly using NR is used for Figure 8 of LIGO-P2000021, and nowhere else in either paper. It has luminosity_distance set to 1 Mpc and redshift close to 0 for all samples, and hence all masses labeled _source are not actually corrected for redshift. The radiated energy also does not have the redshift correction. The remnant BH properties and peak luminosity are not computed using spin evolution and are thus slightly less accurate than for other runs.
    • The set of PE samples for NRSur7dq4_nospin and NRSur7dq4_spinaligned are used for computing Bayes Factors for spin and for higher multipoles, and for the hierarchical analysis described in Section 5.2.1 and Figures 11 and 12 of LIGO-P2000021.
  • GW190521_Ringdown_samples.tgz is a tarball containing nine more h5 files with posterior samples for the ringdown analysis described in Section 3.2 and Figure 9 of LIGO-P2000021.
  • GW190521_Implications_figure_data.tgz is a tarball containing additional data needed to make Figures 5, 10, 11, 12, and 13 in LIGO-P2000021, including skymaps fits files.
  • GW190521_discovery_Fig1.tgz is a tarball containing additional data needed to make Figure 1 in LIGO-P2000020.
  • GW190521_discovery_figs_pdf.tgz - tarball containing all the figures from the GW190521 discovery paper LIGO-P2000020, in pdf.
  • GW190521_Implications_Figures_pdf.tgz - tarball containing all the figures in the GW190521 implications paper LIGO-P2000021, in pdf.
  • GW190521_md5sums.txt - containing md5sums for each of the h5 files.

In this notebook we will primarily use the data in the file GW190521_posterior_samples.h5 .

The 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.7.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.

In [1]:
%matplotlib inline

# import useful python packages
import numpy as np
import matplotlib.pyplot as plt
import h5py

Some simple stuff with "vanilla" h5py

In [2]:
# read in the data
fn = "GW190521_posterior_samples.h5"
data = h5py.File(fn,'r')
In [3]:
# browse the entire file
print('Contents of',fn,':')
for k1 in data.keys():
    for k2 in data[k1].keys():
        try:
            if (data[k1][k2].size == 1):
                print(k1,k2,data[k1][k2][0])
            else:
                print(k1,k2,data[k1][k2].shape)
        except:
            try:
                print(k1,k2,data[k1][k2].keys())
            except:
                try:
                    print(k1,k2,data[k1][k2])
                except:
                    print(k1,k2,'No data')
    
Contents of GW190521_posterior_samples.h5 :
IMRPhenomPv3HM approximant b'IMRPhenomPv3HM'
IMRPhenomPv3HM calibration_envelope <KeysViewHDF5 ['H1', 'L1', 'V1']>
IMRPhenomPv3HM config_file <KeysViewHDF5 ['analysis', 'bayeswave', 'condor', 'data', 'datafind', 'engine', 'input', 'lalinference', 'ligo-skymap-from-samples', 'ligo-skymap-plot', 'mpi', 'paths', 'ppanalysis', 'resultspage', 'singularity', 'skyarea', 'statevector']>
IMRPhenomPv3HM injection_data <KeysViewHDF5 ['injection_values']>
IMRPhenomPv3HM meta_data <KeysViewHDF5 ['meta_data', 'other', 'sampler']>
IMRPhenomPv3HM posterior_samples (49263,)
IMRPhenomPv3HM priors <KeysViewHDF5 ['calibration', 'samples']>
IMRPhenomPv3HM psds <KeysViewHDF5 ['H1', 'L1', 'V1']>
IMRPhenomPv3HM skymap <KeysViewHDF5 ['data', 'meta_data']>
IMRPhenomPv3HM version b'3.0'
NRSur7dq4 approximant b'NRSur7dq4'
NRSur7dq4 calibration_envelope <KeysViewHDF5 ['H1', 'L1', 'V1']>
NRSur7dq4 config_file <KeysViewHDF5 ['analysis', 'bayeswave', 'condor', 'data', 'datafind', 'engine', 'input', 'lalinference', 'ligo-skymap-from-samples', 'ligo-skymap-plot', 'mpi', 'paths', 'ppanalysis', 'resultspage', 'singularity', 'skyarea', 'statevector']>
NRSur7dq4 injection_data <KeysViewHDF5 ['injection_values']>
NRSur7dq4 meta_data <KeysViewHDF5 ['meta_data', 'other', 'sampler']>
NRSur7dq4 posterior_samples (65723,)
NRSur7dq4 priors <KeysViewHDF5 ['calibration', 'samples']>
NRSur7dq4 psds <KeysViewHDF5 ['H1', 'L1', 'V1']>
NRSur7dq4 version b'3.0'
SEOBNRv4PHM approximant b'SEOBNRv4PHM'
SEOBNRv4PHM calibration_envelope <KeysViewHDF5 []>
SEOBNRv4PHM config_file <KeysViewHDF5 []>
SEOBNRv4PHM injection_data <KeysViewHDF5 ['injection_values']>
SEOBNRv4PHM meta_data <KeysViewHDF5 ['meta_data', 'sampler']>
SEOBNRv4PHM posterior_samples (19822,)
SEOBNRv4PHM priors <KeysViewHDF5 ['samples']>
SEOBNRv4PHM psds <KeysViewHDF5 ['H1', 'L1', 'V1']>
SEOBNRv4PHM version b'No version information found'
history command_line b'Generated by running the following script: summarycombine --webdir /home/alan.weinstein/public_html/O3/PE/GW190521/combined --samples /home/alan.weinstein/public_html/O3/PE/GW190521/NRSur7dq4_combined/modified_posterior_samples.h5 /home/alan.weinstein/projects/GW190521/BWPSD-rerun/PhenomPv3HM_combined/modified_posterior_samples.h5 /home/charlie.hoy/public_html/public_release/S190521g/v2/SEOBNRv4PHM_combined/modified_posterior_samples.h5 --gw --no_conversion'
history creator b'alan.weinstein'
history gps_creation_time 1280850188.790526
history program b'summarycombine'
version environment b'/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37-20200806'
version packages (489,)
version pesummary b'0.7.0'
In [4]:
# print out top-level data structures: waveform models
print(' ')
print('Top-level data structures:',data.keys())

# extract posterior samples for preferred waveform model
waveform = 'NRSur7dq4'
print(' ')
print('Data products for preferred waveform model,',waveform,':',data[waveform].keys())

# parameter names:
posterior_samples = data[waveform]['posterior_samples']
pnames = posterior_samples.dtype.fields.keys()
print('parameter names:',pnames)

# extract the samples data into an numpy array:
samples = np.array(posterior_samples).T
 
Top-level data structures: <KeysViewHDF5 ['IMRPhenomPv3HM', 'NRSur7dq4', 'SEOBNRv4PHM', 'history', 'version']>
 
Data products for preferred waveform model, NRSur7dq4 : <KeysViewHDF5 ['approximant', 'calibration_envelope', 'config_file', 'injection_data', 'meta_data', 'posterior_samples', 'priors', 'psds', 'version']>
parameter names: dict_keys(['H1_matched_filter_abs_snr', 'H1_matched_filter_snr_angle', 'H1_optimal_snr', 'H1_spcal_amp_0', 'H1_spcal_amp_1', 'H1_spcal_amp_2', 'H1_spcal_amp_3', 'H1_spcal_amp_4', 'H1_spcal_amp_5', 'H1_spcal_amp_6', 'H1_spcal_amp_7', 'H1_spcal_amp_8', 'H1_spcal_amp_9', 'H1_spcal_phase_0', 'H1_spcal_phase_1', 'H1_spcal_phase_2', 'H1_spcal_phase_3', 'H1_spcal_phase_4', 'H1_spcal_phase_5', 'H1_spcal_phase_6', 'H1_spcal_phase_7', 'H1_spcal_phase_8', 'H1_spcal_phase_9', 'L1_matched_filter_abs_snr', 'L1_matched_filter_snr_angle', 'L1_optimal_snr', 'L1_spcal_amp_0', 'L1_spcal_amp_1', 'L1_spcal_amp_2', 'L1_spcal_amp_3', 'L1_spcal_amp_4', 'L1_spcal_amp_5', 'L1_spcal_amp_6', 'L1_spcal_amp_7', 'L1_spcal_amp_8', 'L1_spcal_amp_9', 'L1_spcal_phase_0', 'L1_spcal_phase_1', 'L1_spcal_phase_2', 'L1_spcal_phase_3', 'L1_spcal_phase_4', 'L1_spcal_phase_5', 'L1_spcal_phase_6', 'L1_spcal_phase_7', 'L1_spcal_phase_8', 'L1_spcal_phase_9', 'V1_matched_filter_abs_snr', 'V1_matched_filter_snr_angle', 'V1_optimal_snr', 'V1_spcal_amp_0', 'V1_spcal_amp_1', 'V1_spcal_amp_2', 'V1_spcal_amp_3', 'V1_spcal_amp_4', 'V1_spcal_amp_5', 'V1_spcal_amp_6', 'V1_spcal_amp_7', 'V1_spcal_amp_8', 'V1_spcal_amp_9', 'V1_spcal_phase_0', 'V1_spcal_phase_1', 'V1_spcal_phase_2', 'V1_spcal_phase_3', 'V1_spcal_phase_4', 'V1_spcal_phase_5', 'V1_spcal_phase_6', 'V1_spcal_phase_7', 'V1_spcal_phase_8', 'V1_spcal_phase_9', 'azimuth', 'cosalpha', 'cos_theta_jn', 'deltalogL', 'deltaloglH1', 'deltaloglL1', 'deltaloglV1', 'log_likelihood', 'log_prior', 'logw', 'network_matched_filter_snr', 'network_optimal_snr', 'phase', 'phi_12', 'phi_jl', 'mass_ratio', 't0', 'geocent_time', 'ra', 'dec', 'luminosity_distance', 'psi', 'chirp_mass', 'a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'inverted_mass_ratio', 'mass_1', 'mass_2', 'total_mass', 'symmetric_mass_ratio', 'iota', 'spin_1x', 'spin_1y', 'spin_1z', 'spin_2x', 'spin_2y', 'spin_2z', 'phi_1', 'phi_2', 'chi_eff', 'chi_p', 'final_mass', 'final_spin', 'final_kick', 'cos_tilt_1', 'cos_tilt_2', 'redshift', 'comoving_distance', 'mass_1_source', 'mass_2_source', 'total_mass_source', 'chirp_mass_source', 'final_mass_source', 'radiated_energy', 'H1_time', 'L1_time', 'V1_time', 'H1_matched_filter_snr', 'L1_matched_filter_snr', 'V1_matched_filter_snr', 'cos_iota', 'peak_luminosity'])
In [5]:
# get samples for one of the parameters
dL = samples['luminosity_distance']
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()
dL shape, mean, std = (65723,) 5265.994269405796 1464.1349618502252
In [6]:
# plot the PSDs used in this analysis
H1psd = data[waveform]['psds']['H1']
L1psd = data[waveform]['psds']['L1']
V1psd = data[waveform]['psds']['V1']
colors = {'H1':'r','L1':'c','V1':'m'}

plt.figure()
for ifo in data[waveform]['psds'].keys():
    fr = data[waveform]['psds'][ifo][()].T[0]
    asd = np.sqrt(data[waveform]['psds'][ifo][()].T[1])
    plt.loglog(fr,asd,colors[ifo],label=ifo)
plt.xlabel('frequency (Hz)')
plt.ylabel('strain ASD (strain/rtHz)')
plt.xlim([20,512])
plt.ylim([1.e-24,1.e-19])
plt.grid('both')
plt.legend()
plt.show()

The main results in the paper are based on these waveform models:
['NRSur7dq4','IMRPhenomPv3HM', 'SEOBNRv4PHM']
where the preferred waveform model is 'NRSur7dq4' .

Some of the paper content is drawn from other waveform models:
[ 'NR', , 'NRSur7dq4_Lmax2', 'NRSur7dq4_Lmax2_coprecessing_modes', 'NRSur7dq4_Lmax3', 'NRSur7dq4_Lmax4', 'NRSur7dq4_aligned']
which we won't explore here.

Now extract the parameter values for each of the three waveform models, for Table 1

In [7]:
# compile and print out the data for Table 1 of the implications paper
from IPython.display import HTML, display

# this function returns confidence intervals for a parameter 
def confidence_range(posterior,parameter,level):
    
    posterior_for_parameter = posterior[parameter]
    
    median = np.percentile(posterior_for_parameter,50)
    lower_bound = - median + np.percentile(posterior_for_parameter,(100-level)/2.)
    upper_bound = - median + np.percentile(posterior_for_parameter,100-(100-level)/2.)
    
    # round to 4 digit precision
    nsig = -int(np.floor(np.log10(np.abs(upper_bound)))) + (4 - 1)
    median = np.around(median,nsig)
    lower_bound = np.around(lower_bound,nsig)
    upper_bound = np.around(upper_bound,nsig)
    
    return(median,upper_bound,lower_bound)

################
# read in the posterior samples, get medians and 90% CIs:
fn = "GW190521_posterior_samples.h5"
data = h5py.File(fn,'r')

waveforms = ['NRSur7dq4', 'IMRPhenomPv3HM', 'SEOBNRv4PHM']

inspiral_parameters = ['mass_1_source', 'mass_2_source', 'total_mass_source','chirp_mass_source', 'mass_ratio',
                       'mass_1', 'mass_2', 'total_mass','chirp_mass', 
                       'a_1', 'a_2', 'tilt_1', 'tilt_2' , 'chi_eff', 'chi_p','final_mass_source', 'final_spin',
                       'radiated_energy', 'peak_luminosity','luminosity_distance', 'redshift', 'ra', 'dec'] 

results = []
for parameter in inspiral_parameters:   
    result = [parameter]
    for waveform in waveforms:
        samples = data[waveform]['posterior_samples']
        try:
            r = confidence_range(samples,parameter,90)
            result.append(r)
        except: 
            g=0
    
    results.append(result)
    # print some latex code. It will need to be hacked...
    #print(f'{result[0]} & ${result[1]}'+'^{'+f'{result[2]}'+'}_{'+f'{result[3]}'+'}$ \\\\')

# print a nicely formatted table for reading
print('Parameter',waveforms[0],waveforms[1],waveforms[2])
    
display(HTML(
       '<table><tr>{}</tr></table>'.format(
           '</tr><tr>'.join(
               '<td>{}</td>'.format('</td><td>'.join(str(_) for _ in row)) for row in results)
           )
    ))
    
    
Parameter NRSur7dq4 IMRPhenomPv3HM SEOBNRv4PHM
mass_1_source(84.83, 21.07, -13.9)(90.33, 22.81, -15.66)(98.9, 42.08, -18.79)
mass_2_source(65.86, 16.96, -17.79)(64.94, 15.78, -18.26)(71.13, 21.01, -27.9)
total_mass_source(149.75, 29.15, -17.01)(154.44, 25.07, -15.9)(170.3, 36.49, -23.16)
chirp_mass_source(64.04, 13.02, -7.98)(65.46, 10.99, -7.46)(71.3, 15.01, -9.92)
mass_ratio(0.7942, 0.1852, -0.2911)(0.7253, 0.2431, -0.2874)(0.7423, 0.2286, -0.424)
mass_1(152.09, 31.69, -18.83)(153.97, 39.9, -21.7)(159.12, 86.91, -24.01)
mass_2(119.8, 20.95, -32.14)(112.73, 25.18, -35.15)(116.99, 28.49, -45.06)
total_mass(272.59, 26.24, -27.38)(267.02, 37.28, -32.31)(279.54, 53.61, -35.99)
chirp_mass(116.84, 12.13, -14.27)(113.93, 15.13, -17.28)(118.13, 16.76, -17.5)
a_1(0.6921, 0.2744, -0.6211)(0.6506, 0.3174, -0.5748)(0.8002, 0.1844, -0.5826)
a_2(0.7278, 0.2436, -0.6351)(0.5276, 0.4226, -0.4764)(0.538, 0.4113, -0.4795)
tilt_1(1.423, 1.111, -0.934)(1.394, 1.113, -0.856)(1.4168, 0.8573, -0.7928)
tilt_2(1.4899, 0.9986, -0.9672)(1.539, 1.093, -1.008)(1.618, 1.057, -1.041)
chi_eff(0.0812, 0.273, -0.3633)(0.059, 0.3074, -0.3932)(0.0635, 0.3396, -0.3474)
chi_p(0.6775, 0.2498, -0.3666)(0.5977, 0.3263, -0.4371)(0.7392, 0.2135, -0.3988)
final_mass_source(142.17, 28.02, -16.37)(147.23, 23.28, -14.69)(162.48, 35.25, -22.02)
final_spin(0.7201, 0.09005, -0.11561)(0.7197, 0.1121, -0.1472)(0.7432, 0.1208, -0.1404)
radiated_energy(7.617, 2.205, -1.861)(7.216, 2.739, -2.156)(7.801, 2.817, -2.347)
peak_luminosity(3.6755, 0.6671, -0.9107)(3.5475, 0.7448, -1.0733)(3.5232, 0.7528, -1.4058)
luminosity_distance(5297.0, 2359.0, -2552.0)(4578.0, 1569.0, -1619.0)(3965.0, 1990.0, -1826.0)
redshift(0.8177, 0.2849, -0.3421)(0.7259, 0.197, -0.2195)(0.6446, 0.2541, -0.2599)
ra(3.457, 2.782, -3.414)(3.533, 2.716, -3.498)(3.645, 2.594, -3.583)
dec(-0.782, 1.508, -0.411)(-0.874, 1.57, -0.328)(-0.572, 1.315, -0.608)

Compute probabilities that the source frame masses are in the nominal PISN mass gap

In [8]:
waveforms = ['NRSur7dq4', 'IMRPhenomPv3HM', 'SEOBNRv4PHM']
nsig = 1
for waveform in waveforms:
    samples = data[waveform]['posterior_samples']
    m1 = samples['mass_1_source']
    m2 = samples['mass_2_source']
    Nsamp = m1.size
    print(waveform,'Nsamples',m1.size)
    print('P(m1<=65) =',round(m1[m1<=65.].size/Nsamp*100.,nsig),'%')
    print('P(m1<=50) =',round(m1[m1<=50.].size/Nsamp*100.,nsig),'%')
    print('P(m1<=120) =',round(m1[m1<=120.].size/Nsamp*100.,nsig),'%')
    print('P(m2<=65) =',round(m2[m2<=65.].size/Nsamp*100.,nsig),'%')
    print('P(m2<=50) =',round(m2[m2<=50.].size/Nsamp*100.,nsig),'%')
    print('P(m2<=120) =',round(m2[m2<=120.].size/Nsamp*100.,nsig),'%')

    mm = m1[((m1<65.) | (m1>120.)) & ((m2<65.) | (m2>120.))]
    print('P(m1,m2 not in [65,120]) = ',round(mm.size/m1.size*100.,nsig),'%')
    mm = m1[((m1>65.) & (m1<120.)) | ((m2>65.) & (m2<120.))]
    print('P(m1 or m2 in [65,120]) = ',round(mm.size/m1.size*100.,nsig),'%')
    print(' ')
NRSur7dq4 Nsamples 65723
P(m1<=65) = 0.3 %
P(m1<=50) = 0.0 %
P(m1<=120) = 99.1 %
P(m2<=65) = 46.2 %
P(m2<=50) = 6.6 %
P(m2<=120) = 100.0 %
P(m1,m2 not in [65,120]) =  1.0 %
P(m1 or m2 in [65,120]) =  99.0 %
 
IMRPhenomPv3HM Nsamples 49263
P(m1<=65) = 0.0 %
P(m1<=50) = 0.0 %
P(m1<=120) = 97.7 %
P(m2<=65) = 50.2 %
P(m2<=50) = 8.8 %
P(m2<=120) = 100.0 %
P(m1,m2 not in [65,120]) =  2.0 %
P(m1 or m2 in [65,120]) =  98.0 %
 
SEOBNRv4PHM Nsamples 19822
P(m1<=65) = 0.0 %
P(m1<=50) = 0.0 %
P(m1<=120) = 88.0 %
P(m2<=65) = 34.8 %
P(m2<=50) = 10.1 %
P(m2<=120) = 100.0 %
P(m1,m2 not in [65,120]) =  9.8 %
P(m1 or m2 in [65,120]) =  90.2 %
 
In [9]:
# we also need some numbers from HM runs, so read those in:
studies_samples_fn = "GW190521_studies_posterior_samples.h5"
studies_samples_data = h5py.File(studies_samples_fn,'r')

hm_waveforms = ['NRSur7dq4_Lmax2','NRSur7dq4_Lmax4']
hm_parameters = ['mass_1_source', 'mass_2_source', 'mass_ratio', 'chi_eff', 'chi_p','luminosity_distance', 'redshift','iota']
hm_results = []
for parameter in hm_parameters:   
    result = [parameter]
    for waveform in hm_waveforms:
        samples = studies_samples_data[waveform]['posterior_samples']
        try:
            r = confidence_range(samples,parameter,90)
            result.append(r)
        except: 
            g=0
    
    hm_results.append(result)

# print a nicely formatted table for reading
print('Parameter',hm_waveforms[0],hm_waveforms[1])
    
display(HTML(
       '<table><tr>{}</tr></table>'.format(
           '</tr><tr>'.join(
               '<td>{}</td>'.format('</td><td>'.join(str(_) for _ in row)) for row in hm_results)
           )
    ))
    
Parameter NRSur7dq4_Lmax2 NRSur7dq4_Lmax4
mass_1_source(104.22, 33.22, -25.7)(83.49, 19.53, -13.03)
mass_2_source(63.52, 23.22, -20.23)(65.5, 16.51, -15.42)
mass_ratio(0.6178, 0.3268, -0.2706)(0.8019, 0.1746, -0.2664)
chi_eff(0.0608, 0.2743, -0.2898)(0.0676, 0.2801, -0.3031)
chi_p(0.7505, 0.2073, -0.3944)(0.6609, 0.2726, -0.4004)
luminosity_distance(3651.0, 2440.0, -1785.0)(5392.0, 2199.0, -2358.0)
redshift(0.6022, 0.3131, -0.2598)(0.829, 0.2654, -0.3125)
iota(1.428, 1.26, -1.024)(1.013, 1.902, -0.804)
In [10]:
# compute Bayes Factors for spin precession and HMs.

def get_BF_data(BF_data):
    BFn = BF_data['meta_data']['sampler']['ln_bayes_factor'][0]
    lnz = BF_data['meta_data']['sampler']['ln_evidence'][0]
    H   = BF_data['meta_data']['other']['information_nats'][0]
    nsamp = BF_data['meta_data']['sampler']['nsamples'][0]
    nlive = BF_data['meta_data']['other']['number_live_points'][0]
    return BFn, lnz, H, nsamp, nlive

# # get BF(s/n), ln(evidence), information gain (nats) Nsamples and live points for the reference model:
BF1, lnz1, H1, nsamp1, nlive1 = get_BF_data(data['NRSur7dq4'])
print('NRSur PHM','BFs/n, lnz, H, nsamp, nlive =', BF1, lnz1, round(H1,2), nsamp1, nlive1,'\n')

# models to compare:
BFinputs = {'aligned spin': studies_samples_data['NRSur7dq4_spinaligned'], 
            'No spin':      studies_samples_data['NRSur7dq4_nospin'], 
            'Only l2m2':    studies_samples_data['NRSur7dq4_l2m2']}

print("GW190521 Bayes factor of NRSur PHM vs:")
for k in BFinputs.keys():
    # get BF(s/n), ln(evidence), information gain (nats) Nsamples and live points for alternative models
    BF2, lnz2, H2, nsamp2, nlive2 = get_BF_data(BFinputs[k])
    print(k,'BFs/n, lnz, H, nsamp, nlive =', BF2, lnz2, round(H2,2), nsamp2, nlive2)
    # compute bayes factor between models
    bf = np.exp(lnz1 - lnz2) 
    # compute theoretical error in the natural log of BF (Skilling p.844)
    delta_ln_bf = np.sqrt(H1/nlive1+ H2/nlive2)
    # compute BF upper and lower limits
    bf_p = bf*np.exp(delta_ln_bf)
    bf_m = bf*np.exp(-delta_ln_bf)
    # print out the log10 BFs
    print('\t log_10 BF = %.3f +%.3f -%.3f' % (np.log10(bf), np.log10(bf_p/bf), np.log10(bf/bf_m)))
NRSur PHM BFs/n, lnz, H, nsamp, nlive = 84.49 -11938.02 21.71 65723 2048 

GW190521 Bayes factor of NRSur PHM vs:
aligned spin BFs/n, lnz, H, nsamp, nlive = 82.06 -11940.45 20.78 63594 2048
	 log_10 BF = 1.055 +0.063 -0.063
No spin BFs/n, lnz, H, nsamp, nlive = 82.37 -11940.14 20.45 55894 2048
	 log_10 BF = 0.921 +0.062 -0.062
Only l2m2 BFs/n, lnz, H, nsamp, nlive = 85.36 -11937.15 22.47 72060 2048
	 log_10 BF = -0.378 +0.064 -0.064
In [11]:
# Done with simple stuff using `h5py` ; release memory for the data
del data
del studies_samples_data

Now use PESummary (v0.7.0 or greater) to read the data files as well as for plotting.

In [12]:
# 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.7.0.
import pesummary
from pesummary.gw.file.read import read
print(pesummary.__version__)
0.7.0

There are 3 different waveform models that were used to analyze GW190521 and they are all stored in the data file.

In [13]:
fn = "GW190521_posterior_samples.h5"
data = read(fn)
labels = data.labels
print(labels)
['IMRPhenomPv3HM', 'NRSur7dq4', 'SEOBNRv4PHM']

To illustrate the data structure we'll pick the preferred waveform model and plot its parameter samples.

In [14]:
waveform = "NRSur7dq4"
samples_dict = data.samples_dict
posterior_samples = samples_dict[waveform]
prior_samples = data.priors["samples"][waveform]
parameters = posterior_samples.keys()
print(parameters)
dict_keys(['H1_matched_filter_abs_snr', 'H1_matched_filter_snr_angle', 'H1_optimal_snr', 'H1_spcal_amp_0', 'H1_spcal_amp_1', 'H1_spcal_amp_2', 'H1_spcal_amp_3', 'H1_spcal_amp_4', 'H1_spcal_amp_5', 'H1_spcal_amp_6', 'H1_spcal_amp_7', 'H1_spcal_amp_8', 'H1_spcal_amp_9', 'H1_spcal_phase_0', 'H1_spcal_phase_1', 'H1_spcal_phase_2', 'H1_spcal_phase_3', 'H1_spcal_phase_4', 'H1_spcal_phase_5', 'H1_spcal_phase_6', 'H1_spcal_phase_7', 'H1_spcal_phase_8', 'H1_spcal_phase_9', 'L1_matched_filter_abs_snr', 'L1_matched_filter_snr_angle', 'L1_optimal_snr', 'L1_spcal_amp_0', 'L1_spcal_amp_1', 'L1_spcal_amp_2', 'L1_spcal_amp_3', 'L1_spcal_amp_4', 'L1_spcal_amp_5', 'L1_spcal_amp_6', 'L1_spcal_amp_7', 'L1_spcal_amp_8', 'L1_spcal_amp_9', 'L1_spcal_phase_0', 'L1_spcal_phase_1', 'L1_spcal_phase_2', 'L1_spcal_phase_3', 'L1_spcal_phase_4', 'L1_spcal_phase_5', 'L1_spcal_phase_6', 'L1_spcal_phase_7', 'L1_spcal_phase_8', 'L1_spcal_phase_9', 'V1_matched_filter_abs_snr', 'V1_matched_filter_snr_angle', 'V1_optimal_snr', 'V1_spcal_amp_0', 'V1_spcal_amp_1', 'V1_spcal_amp_2', 'V1_spcal_amp_3', 'V1_spcal_amp_4', 'V1_spcal_amp_5', 'V1_spcal_amp_6', 'V1_spcal_amp_7', 'V1_spcal_amp_8', 'V1_spcal_amp_9', 'V1_spcal_phase_0', 'V1_spcal_phase_1', 'V1_spcal_phase_2', 'V1_spcal_phase_3', 'V1_spcal_phase_4', 'V1_spcal_phase_5', 'V1_spcal_phase_6', 'V1_spcal_phase_7', 'V1_spcal_phase_8', 'V1_spcal_phase_9', 'azimuth', 'cosalpha', 'cos_theta_jn', 'deltalogL', 'deltaloglH1', 'deltaloglL1', 'deltaloglV1', 'log_likelihood', 'log_prior', 'logw', 'network_matched_filter_snr', 'network_optimal_snr', 'phase', 'phi_12', 'phi_jl', 'mass_ratio', 't0', 'geocent_time', 'ra', 'dec', 'luminosity_distance', 'psi', 'chirp_mass', 'a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'inverted_mass_ratio', 'mass_1', 'mass_2', 'total_mass', 'symmetric_mass_ratio', 'iota', 'spin_1x', 'spin_1y', 'spin_1z', 'spin_2x', 'spin_2y', 'spin_2z', 'phi_1', 'phi_2', 'chi_eff', 'chi_p', 'final_mass', 'final_spin', 'final_kick', 'cos_tilt_1', 'cos_tilt_2', 'redshift', 'comoving_distance', 'mass_1_source', 'mass_2_source', 'total_mass_source', 'chirp_mass_source', 'final_mass_source', 'radiated_energy', 'H1_time', 'L1_time', 'V1_time', 'H1_matched_filter_snr', 'L1_matched_filter_snr', 'V1_matched_filter_snr', 'cos_iota', 'peak_luminosity'])

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.

In [15]:
parameter = "luminosity_distance"
fig = posterior_samples.plot(parameter, type="hist", prior=prior_samples[parameter])
fig.set_size_inches(12, 8)
plt.show()
In [16]:
parameter = "luminosity_distance"
waveforms = ['NRSur7dq4', 'IMRPhenomPv3HM', 'SEOBNRv4PHM']
fig = samples_dict.plot(parameter, type="hist", kde=True, labels=waveforms)
fig.set_size_inches(12, 8)
plt.show()

Make a corner plot (this takes some time):

In [17]:
fig = posterior_samples.plot(type="corner", parameters=["mass_1_source", "mass_2_source", 
                                                        "chirp_mass_source", "total_mass_source"])
plt.show()

The psds that were used for each analysis can also be extracted from this file and plotted

In [18]:
psd = data.psd[waveform]
fig = psd.plot(fmin=11.)
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

In [19]:
from pesummary.gw.plots.plot import _calibration_envelope_plot

prior = data.priors["calibration"]["NRSur7dq4"]
calibration = data.calibration["NRSur7dq4"]
ifos = list(calibration.keys())
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(12, 8)
plt.show()

The configuration file that were used for each analysis can also be extracted from this file

In [20]:
config = data.config["NRSur7dq4"]
for i in config.keys():
    print("[{}]".format(i))
    for key, item in config[i].items():
        print("{}={}".format(key, item))   # .decode("utf-8")
    print("\n")
[analysis]
coherence-test=False
dataseed=1234
engine=lalinferencenest
ifos=['H1', 'L1', 'V1']
nparallel=4
osg=False
randomseed=4321
roq=False
service-url=https://gracedb.ligo.org/api/
singularity=False
upload-to-gracedb=False


[bayeswave]


[condor]
accounting_group=ligo.prod.o3.cbc.pe.lalinference
accounting_group_user=max.isi
coherencetest=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/lalinference_coherence_test
combinePTMCMCh5script=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/cbcBayesCombinePTMCMCh5s
computeroqweights=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/lalinference_compute_roq_weights
datafind=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/gw_data_find
gracedb=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/gracedb
lalinferencebambi=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/lalinference_bambi
lalinferencedatadump=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/lalinference_datadump
lalinferencemcmc=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/lalinference_mcmc
lalinferencenest=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/lalinference_nest
lalsuite-install=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37/
ligo-skymap-from-samples=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/ligo-skymap-from-samples
ligo-skymap-plot=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/ligo-skymap-plot
ligolw_print=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/ligolw_print
mergeMCMCscript=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/cbcBayesMCMC2pos
mergeNSscript=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/lalinference_nest2pos
mpirun=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/mpirun
mpiwrapper=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/lalinference_mpi_wrapper
pos_to_sim_inspiral=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/cbcBayesPosToSimInspiral
ppanalysis=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/cbcBayesPPAnalysis
processareas=/usr/bin/process_areas
queue=Priority_PE
resultspage=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/cbcBayesPostProc
segfind=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37//bin/ligolw_segment_query


[data]
channels={'H1': 'H1:DCS-CALIB_STRAIN_CLEAN_SUB60HZ_C01', 'L1': 'L1:DCS-CALIB_STRAIN_CLEAN_SUB60HZ_C01','V1': 'V1:Hrec_hoft_16384Hz'}


[datafind]
types={'H1': 'H1_HOFT_CLEAN_SUB60HZ_C01', 'L1': 'L1_HOFT_CLEAN_SUB60HZ_C01','V1': 'V1Online'}
url-type=file


[engine]
H1-psd=/home/katerina.chatziioannou/BayesWaveO3/GW190521/PSDFix/glitch_median_PSD_forLI_H1.dat
H1-spcal-envelope=/home/ling.sun/public_html/Calibration/Uncertainty/O3C01/LHO/Aug-24-2019_O3_LHO_GPSTime_1242442884_C01_RelativeResponseUncertainty_FinalResults.txt
L1-psd=/home/katerina.chatziioannou/BayesWaveO3/GW190521/PSDFix/glitch_median_PSD_forLI_L1.dat
L1-spcal-envelope=/home/ling.sun/public_html/Calibration/Uncertainty/O3C01/LLO/Aug-24-2019_O3_LLO_GPSTime_1242443764_C01_RelativeResponseUncertainty_FinalResults.txt
V1-psd=/home/katerina.chatziioannou/BayesWaveO3/GW190521/PSDFix/glitch_median_PSD_forLI_V1.dat
V1-spcal-envelope=/home/cbc/pe/O3/calibrationenvelopes/Virgo/V_earlyO3_calibrationUncertaintyEnvelope_magnitude5percent_phase2degrees10microseconds.txt
a_spin1-max=0.99
a_spin2-max=0.99
adapt-temps=
amporder=1
approx=NRSur7dq4pseudoFourPN
chirpmass-max=150.0
chirpmass-min=70.0
comp-max=300.0
comp-min=30.0
distance-max=10000
enable-spline-calibration=
fref=11
maxmcmc=10000
mtotal-min=200
neff=1000
nlive=2048
ntemps=8
progress=
q-max=1.0
q-min=0.17
resume=
seglen=8.0
spcal-nodes=10
srate=1024.0
tolerance=0.1


[input]
analyse-all-time=False
events=all
gps-end-time=1242442986
gps-start-time=1242442821
gps-time-file=/home/max.isi/projects/cbc_pe/GW190521g/gpstime.txt
ignore-gracedb-psd=True
ignore-state-vector=True
max-psd-length=10000
minimum_realizations_number=8
padding=16
threshold-snr=3
timeslides=False


[lalinference]
flow={'H1': 11.0, 'L1': 11.0, 'V1': 11.0}


[ligo-skymap-from-samples]
enable-multiresolution=


[ligo-skymap-plot]
annotate=
contour=50 90


[mpi]
machine-count=8
machine-memory=2000
mpi_task_count=8


[paths]
basedir=/home/max.isi/projects/cbc_pe/GW190521g/prod_nr7dq4_nest_ao1_8s_newpsd/
daglogdir=/home/max.isi/projects/cbc_pe/GW190521g/prod_nr7dq4_nest_ao1_8s_newpsd
webdir=/home/max.isi/public_html/cbc_pe/GW190521g/prod_nr7dq4_nest_ao1_8s_newpsd/


[ppanalysis]


[resultspage]
skyres=0.5


[singularity]


[skyarea]
maxpts=2000


[statevector]
bits=['Bit 0', 'Bit 1', 'Bit 2']
state-vector-channel={'H1': 'H1:DCS-CALIB_STATE_VECTOR', 'L1': 'L1:DCS-CALIB_STATE_VECTOR','V1': 'V1:DQ_ANALYSIS_STATE_VECTOR'}