The DCC document https://dcc.ligo.org/T2100512 contains the maximum posterior draw from the Power Law + Dip + Break (ind) and Power Law + Dip + Break (pair) models highlighted in “The population of merging compact binaries inferred using gravitational waves through GWTC-3”.
The models are implemented in this fork of gwpopulation.
The first file, O1O2O3all_mass_g_iid_mag_iid_tilt_powerlaw_redshift_maxP_sample.json is a human-readable file that contains the maximum-posterior draw for the inferred hyperparameters of Power Law + Dip + Break (ind) model when fit to data from GWTC-3. This includes values for the rate and spin hyperparameters. The following examples use the independent-pairing model, but you can swap it out for the pairing function model if you use the files that contain mass_h instead of mass_g and use the matter_matters_pairing function from the above-linked fork of gwpopulation.
You can evaluate the models by plugging the json files into the above linked function, like so:
from gwpopulation.models.mass import matter_matters_primary_secondary_independent # change to matter_matters_pairing if you use model h
from bilby.hyper.model import Model
import pandas as pd
import numpy as np
model = Model([matter_matters_primary_secondary_independent]) # again, change to matter_matters_pairing if you use model h
# load it in your favorite way, I prefer pandas
post = pd.io.json.read_json("O1O2O3all_mass_g_iid_mag_iid_tilt_powerlaw_redshift_maxP_sample.json", orient='index', typ='series')
model.parameters.update(post)
and then calling the model on whatever data you like. For example,
ms = np.logspace(0,np.log10(100), 2000)
mass_1_grid, mass_2_grid = np.meshgrid(ms, ms)
data = dict(mass_1=mass_1_grid, mass_2=mass_2_grid)
unnormed_prob = model.prob(data)
prob = unnormed_prob * maxp['normalization']
rate = prob * maxp['rate']
If you don’t want to clone the above-linked repo and install that version of gwpopulation, you can also find an h5 file of dR/dm1dm2 evaluated at this maximum posterior draw using essentially the code above. It has 4 entries:
mass : mass grid on which the model is evaluated. It is log-spaced from 1 to 100 solar masses. This is the same for m1 and m2 so for simplicity I just provide one 1D array here.p_mass_1_mass_2 : 2D array of the differential merger rate evaluated at each m1 and m2 (dR/dm1dm2). m1 is on axis=1 and m2 is on axis 0. For convenience, I’ve also separately marginalized over each component mass for the next two entries.p_mass_1 : This is actually dR/dm1. it is p_mass_1_mass_2 marginalized over m2p_mass_2 : This is actually dR/dm2. it is p_mass_1_mass_2 marginalized over m1pdf : dR/dm where “m” is generic component mass. It is obtained by taking the average of p_mass_1 and p_mass_2 at each mass point.You should therefore be able to plot the model using something like the following example:
from astropy.table import Table
# use mass_h if you want the powerlaw pairing model
filename = "O1O2O3all_mass_g_iid_mag_iid_tilt_powerlaw_redshift_mass_maxP_pdf.h5"
lines = Table.read(filename)
mass = lines['mass']
p_mass_1 = lines['p_mass_1']
plt.figure()
plt.semilogy(mass, p_mass_1)
plt.xlabel('primary mass [Msun]')
plt.ylabel('differential merger rate [Msun^-1 Gpc^-3 yr^-1]')
plt.savefig('p_m1_maxP.png')
plt.close()
Or, if you want to use h5py:
import h5py
# use mass_h if you want the powerlaw pairing model
filename = "O1O2O3all_mass_g_iid_mag_iid_tilt_powerlaw_redshift_mass_maxP_pdf.h5"
f = h5py.File(filename, 'r')
mass = f['lines']['mass']
p_mass_1 = f['lines']['p_mass_1']
plt.figure()
plt.semilogy(mass, p_mass_1)
plt.xlabel('primary mass [Msun]')
plt.ylabel('differential merger rate [Msun^-1 Gpc^-3 yr^-1]')
plt.savefig('p_m1_maxP.png')
plt.close()
You can also integrate the differential merger rate over the specific component mass regions to get rates for certain source types. For example, a rate of BNS-like events can be obtained by assuming a maximum NS mass (I would suggest using the value of NSmaxin the json file) and integrating p(m1) up to that value.