Analyzing public parameter estimation samples up to GWTC-3

The primary use case for gwpopulation_pipe is analyzing the gravitational-wave transient catalog (GWTC).

This will generate results equivalent to a subset of the result in the LIGO-Virgo-Kagra astrophysical distributions paper accompanying GWTC-3.

Warning

This example was written using gwpopulation_pipe==0.3.2 and may not work with the latest version of the code.

First, we need to obtain various pieces of data.

Downloading single-event fiducial posterior samples

The samples from each version of GWTC are available individually. The following samples will download the data and unpack them into individual directories. They can also be manually downloaded for GWTC-1, GWTC-2 and GWTC-3.

DistancePrior="nocosmo"

# fetch O1/O2 samples using GWTC-2.1 release
mkdir gwtc1_samples
cd gwtc1_samples
declare -a GWTC1_EVENTS=(
    "150914_095045" "151012_095443" "151226_033853" "170104_101158" "170608_020116"
    "170729_185629" "170809_082821" "170814_103043" "170818_022509" "170823_131358"
)
for EVENT in "${GWTC1_EVENTS[@]}"
do
  wget -nc "https://zenodo.org/records/6513631/files/IGWN-GWTC2p1-v2-GW${EVENT}_PEDataRelease_mixed_${DistancePrior}.h5"
done
cd -

# fetch O3a samples
mkdir o3a_samples
cd o3a_samples
declare -a O3A_EVENTS=(
    "190403_051519" "190408_181802"
    "190412_053044" "190413_052954" "190413_134308" "190421_213856" "190426_190642"
    "190503_185404" "190512_180714" "190513_205428" "190514_065416" "190517_055101"
    "190519_153544" "190521_030229" "190521_074359" "190527_092055" "190602_175927"
    "190620_030421" "190630_185205" "190701_203306" "190706_222641" "190707_093326"
    "190708_232457" "190719_215514" "190720_000836" "190725_174728" "190727_060333"
    "190728_064510" "190731_140936" "190803_022701" "190805_211137" "190814_211039"
    "190828_063405" "190828_065509" "190910_112807" "190915_235702" "190916_200658"
    "190917_114630" "190924_021846" "190925_232845" "190926_050336" "190929_012149"
    "190930_133541"
)
for EVENT in "${O3A_EVENTS[@]}"
do
  wget -nc "https://zenodo.org/records/6513631/files/IGWN-GWTC2p1-v2-GW${EVENT}_PEDataRelease_mixed_${DistancePrior}.h5"
done
cd -

# fetch O3b samples
echo "Fetching O3b samples from zenodo.org/record/5546663"
mkdir o3b_samples
cd o3b_samples
declare -a O3B_EVENTS=(
    "191103_012549" "191105_143521" "191109_010717" "191113_071753" "191126_115259"
    "191127_050227" "191129_134029" "191204_110529" "191204_171526" "191215_223052"
    "191216_213338" "191219_163120" "191222_033537" "191230_180458" "200105_162426"
    "200112_155838" "200115_042309" "200128_022011" "200129_065458" "200202_154313"
    "200208_130117" "200208_222617" "200209_085452" "200210_092255" "200216_220804"
    "200219_094415" "200220_061928" "200220_124850" "200224_222234" "200225_060421"
    "200302_015811" "200306_093714" "200308_173609" "200311_115853" "200316_215756"
    "200322_091133"
)
for EVENT in "${O3B_EVENTS[@]}"
do
  wget -nc "https://zenodo.org/record/5546663/files/IGWN-GWTC3p0-v1-GW${EVENT}_PEDataRelease_mixed_${DistancePrior}.h5"
done
cd -

Warning

This will download > 20 GB of data.

Download sensitivity data products

The sensitivity of the gravitational-wave detector network is assessed using simulated signals that are injected into the data stream and recovered with the same search pipelines that identify real gravitational-wave signals. Since for this example we’re just looking at binary black hole systems we will only download the data relevant to those systems.

$ wget -nc https://zenodo.org/record/5636816/files/o1%2Bo2%2Bo3_bbhpop_real%2Bsemianalytic-LIGO-T2100377-v2.hdf5

Additional sensitivity products for other classes of system can also be downloaded from the same location.

Setup the configuration file

Now we can write our configuration file gwtc3-bbh.ini for gwpopulation_pipe.

Configuration file

################################################################################
## Generic arguments
################################################################################

run-dir=GWTC-3
log-dir=GWTC-3/logs
label=bbh
# This option is only relevant if you are running on a cluster that requires
# a user accounting tag.
user=colm.talbot
# use the vt-file from the docker image used for testing that already has the
# FAR and SNR cuts applied
vt-file=/test-data/sensitivity/o1+o2+o3_bbhpop_real+semianalytic-LIGO-T2100377-v2-far-1-snr-10.hdf5
vt-ifar-threshold=1.0
vt-snr-threshold=10.0
vt-function=injection_resampling_vt
prior-file=gwtc3-bbh.prior
# the test image does not have a GPU, so we set this to 0
request-gpu=0
backend=jax

################################################################################
## Analysis models
# A tensor product will be done with the models listed here, i.e., if you list
# two mass models, two spin magnitude models, one spin tilt model, and two
# redshift models, a total of 8 analysis jobs will be set up
################################################################################

all-models={mass: [two_component_primary_mass_ratio], redshift: [powerlaw]}

################################################################################
## Data collection arguments
################################################################################

parameters=[mass_1, mass_ratio, a_1, a_2, cos_tilt_1, cos_tilt_2, redshift]
ignore=[GW170817, S190425z, GW200115, 190403, 190426, 190514, S190814bv, 190916, GW190917, 190926, 191113, 191126, GW191219, 200105, 200115, 200210, 200220, 200306, 200308, 200322]
# the sample files point to the location in the docker image used for testing
# we downsample the events by month to make the tests faster but still cover
# each observing run
sample-regex={"GWTC1": "/test-data/gwtc1-downsampled/*GW1708*.h5", "O3a": "/test-data/o3a-downsampled/*GW1908*.h5", "O3b": "/test-data/o3b-downsampled/*GW1911*.h5"}
preferred-labels=[PrecessingSpinIMRHM, C01:IMRPhenomXPHM, IMRPhenomXPHM]
plot=True
data-label=posteriors
distance-prior=comoving
mass-prior=flat-detector-components

################################################################################
## Arguments describing analysis jobs
################################################################################

max-redshift=1.9
minimum-mass=2
maximum-mass=100
# sampler settings chosen to make the tests faster
# these are likely not good settings for a real analysis
sampler=nestle
sampler-kwargs={nlive:50}
vt-parameters=[mass, redshift]
enforce-minimum-neffective-per-event=False

################################################################################
## Post processing arguments
################################################################################

post-plots=True
make-summary=False
n-post-samples=5000

This references a prior file which should contain the following

Prior file

alpha = Uniform(minimum=-4, maximum=12, name='alpha', latex_label='$\\alpha$')
alpha_1 = Uniform(minimum=-4, maximum=12, name='alpha_1', latex_label='$\\alpha_{1}$')
alpha_2 = Uniform(minimum=-4, maximum=12, name='alpha_2', latex_label='$\\alpha_{2}$')
beta = Uniform(minimum=-2, maximum=7, name='beta', latex_label='$\\beta_{q}$')
mmax = Uniform(minimum=62, maximum=100, name='mmax', latex_label='$m_{\\max}$')
mmin = Uniform(minimum=2, maximum=3.5, name='mmin', latex_label='$m_{\\min}$')
break_fraction = Uniform(minimum=0, maximum=1, name='break_fraction', latex_label='$b$')
lam = Uniform(minimum=0, maximum=1, name='lambda', latex_label='$\\lambda_{m}$')
lam_1 = Uniform(minimum=0, maximum=1, name='lambda_1', latex_label='$\\lambda_{m1}$')
mpp = Uniform(minimum=20, maximum=50, name='mpp', latex_label='$\\mu_{m}$')
sigpp = Uniform(minimum=1, maximum=10, name='sigpp', latex_label='$\\sigma_{m}$')
mpp_1 = Uniform(minimum=5, maximum=20, name='mpp_1', latex_label='$\\mu_{m1}$')
mpp_2 = Uniform(minimum=20, maximum=100, name='mpp_2', latex_label='$\\mu_{m2}$')
sigpp_1 = Uniform(minimum=0.1, maximum=5, name='sigpp_1', latex_label='$\\sigma_{m1}$')
sigpp_2 = Uniform(minimum=1, maximum=25, name='sigpp_2', latex_label='$\\sigma_{m2}$')
delta_m = Uniform(minimum=0, maximum=10, name='delta_m', latex_label='$\\delta_{m}$')
amax = 1
mu_chi = Uniform(minimum=0, maximum=1, name='mu_chi', latex_label='$\\mu_{\\chi}$')
sigma_chi = Uniform(minimum=0.005, maximum=0.25, name='sigma_chi', latex_label='$\\sigma^{2}_{\\chi}$')
alpha_chi = Constraint(minimum=1, maximum=1e5)
beta_chi = Constraint(minimum=1, maximum=1e5)
xi_spin = Uniform(minimum=0, maximum=1, name='xi_spin', latex_label='$\\zeta$')
sigma_spin = Uniform(minimum=1e-1, maximum=4e0, name='sigma_t', latex_label='$\\sigma_t$')
lamb = Uniform(minimum=-2, maximum=8, name='lamb', latex_label="$\\lambda_z$")
gaussian_mass_maximum = 100.0

Arbitrary models can be passed as python paths, e.g., my_population_inference_package.my_population_inference_model. The only limitation is that this function must follow the standard format for input arguments an return values, e.g.,

def my_population_inference_model(dataset, parameter_1, parameter_2, parameter_3):
    return (dataset["mass_1"] / parameter_1) ** parameter_2 + parameter_3

Using custom population models may make some of the output filenames a little clunky.

Run `gwpopulation_pipe`

Running the following will set up everything needed to perform the analysis.

$ gwpopulation_pipe gwtc3-bbh.ini

49 bilby_pipe INFO    : dag file written to GWTC-3/submit/o1o2o3.dag
49 bilby_pipe INFO    : shell script written to GWTC-3/submit/o1o2o3.sh
49 bilby_pipe INFO    : Now run condor_submit_dag GWTC-3/submit/o1o2o3.dag

The three lines of output will direct you to the submit directory that contains a bash script with all of the required commands and also files required to submit to a cluster running the HTCondor schedule manager. If you are running the jobs locally you can just call

$ bash GWTC-3/submit/o1o2o3.sh

to run all of the required stages in serial.

If these jobs run to completion you should see a set of results in the result and summary directories.