pe-configurator 0.1.0

Tools to identify waveform settings and prior bounds for gravitational wave parameter estimation.

Introduction

In this repo we provide scripts that automatically try to determine some settings that must be applied for time domain EOB waveforms to run.

Top-level description

The main script to use is proc_samples.py. Given a preliminary PE run, it will perform some basic checks and then generate recommended PE settings. It can also optionally generate a json format result for automatic ingestion by the robot and a more detailed standalone report (still under construction).

NB: the script is intended to be used for "vanilla" BBH events.

Usage

Here is a typical usage example:

peconfigurator --HM /home/zoheyr.doctor/public_html/O3/O3aCatalog/data_release/all_posterior_samples/GW190413_052954_comoving.h5  --dataset C01:SEOBNRv4P --json_file test.json --report_file test.pdf

For all the available options, please check

peconfigurator -h

If all goes right it will print something like

Processing....
Estimating seglen based on posterior samples
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10885/10885 [02:22<00:00, 76.48it/s]
Estimating srate based on posterior samples
Using ell = 3 for Nyquist check
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10885/10885 [00:00<00:00, 30557.61it/s]
Estimating f_start based on posterior samples
Done processing

############################## SUMMARY ##############################
Posterior parameters
----------------------------------------
         total mass             q            Mc         chi1z         chi2z            dL
count  10885.000000  10885.000000  10885.000000  10885.000000  10885.000000  10885.000000
mean      93.716478      0.721174     39.633770      0.001211     -0.027717   3628.546845
std        8.950171      0.171014      3.951967      0.261996      0.297189   1048.258790
min       66.155982      0.145312     25.004179     -0.967360     -0.953884   1026.785460
25%       87.516019      0.595916     37.017804     -0.131398     -0.196356   2904.252323
50%       93.068042      0.729647     39.459142      0.005121     -0.012037   3506.855087
75%       99.238665      0.863340     42.139057      0.152395      0.134719   4216.193654
max      124.959332      0.999992     53.957413      0.906870      0.973028   7975.515016
----------------------------------------
Sanity checks
Checking for railing in chirp mass: PASSED
Checking for railing in luminosity distance: PASSED
Checking that f_ref is not above f_MECO: PASSED

++++++++++++++++++++++++++++++ Suggested settings ++++++++++++++++++++++++++++++
Real required srate  753.7701187444659, as power of 2: 1024
Real seglen: 4.046571094240219, with padding+rollon: 6.44657109424022, as power of 2: 8.0
Check from samples:
Real seglen: 0.8090184427276651, with padding+rollon: 3.209018442727665, as power of 2: 4.0
f_start: 13.333333333333334, as amp_order: 1
Real bounds on chirp mass: [25.004179128030085,53.95741304487518], suggested bounds: [20.00334330242407,64.74889565385021]
Real distance max: 8000.0

The test.json file contains some machine-readable information as shown below:

{
    "srate": 1024,
    "f_start": 13.333333333333334,
    "f_ref": 20.0,
    "amp_order": 1,
    "seglen": 8.0,
    "chirpmass_min": 20.00334330242407,
    "chirpmass_max": 64.74889565385021,
    "metadata": {
        "date": "2020-11-05 04:38:55.172784",
        "command_line_args": "Namespace(HM=True, bounds_tol=0.2, dataset='C01:SEOBNRv4P', ell_max=3, f_ref=20.0, flow=20.0, json_file='test.json', nbins=50, q_min=None, report_file='test.pdf', samples_file='/home/zoheyr.doctor/public_html/O3/O3aCatalog/data_release/all_posterior_samples/GW190413_052954_comoving.h5', tolerance=2.0)"
    }

An example of the generated report can be seen here:

https://ldas-jobs.ligo.caltech.edu/~serguei.ossokine/LVC/O3b_settings/test.html

It is recommended that the report is reviewed carefully before the settings are approved.

Running on many events

To run on many events, a separate script, run_on_many_events.py is provided. The most important settings to it are the path to the repo containing the superevents (specified via --repo) and a json file containing the names of events to run on and the corresponding preferred data set inside the PEsummary metafile (specified via --run_map). An example of such a json file is below:

{
    "S191105e":"EXP1",
    "S191129u":"EXP4"
}

Here is an example of calling the script:

python  run_on_many_events.py --repo /home/gareth.davies/git_repos/o3b_catalog_events/ --run_map ../test_run_map.json --json_output --full_reports --output_dir example

With these options, a separate directory will be created for each superevent inside the example directory. Inside each directory there will be 4 files, corresponding to the raw log (ending in .log), the machine readable json file with recommended settings (ending in .json), and 2 files corresponding to a more verbose report (ending with .ipynb and .html). The html report is the exectued and rendered version of the ipynb file, provided for debugging purposes.

An example of the directory and output can be found here: https://ldas-jobs.ligo.caltech.edu/~serguei.ossokine/LVC/O3b_settings/full_test/

Detailed description of internal logic

This section describes the detailed logic used in deciding the various settings. Many of these settings depend on the masses of the binary and are thus affected by proposed changes to the chirp mass prior which is described below. For every setting we generate 2 estimates: one coming directly from the samples and a conservative one, which takes into account the wider chirp mass prior. In particular, the following strategy is adopted:

1. Estimate the quantity of interest from the samples.
2. Pick the most extreme value of the quantity of interest (e.g. max seglen or min $f_{\rm start}$) and find its parameters (mass ratio and spins)
3. Use either the lower or upper proposed chirp mass bounds (as appropriate) and the mass ratio from step 2 to estimate a new total mass
4. Use this total mass along with the other params (as appropriate) to produce a new estimate of the quantity of interest

Deciding the chirp mass bounds for the prior

The script suggests a chirp mass prior that is wider than the minimum and maximum chirp masses found in the samples, to guard against the moving around of the posterior due to e.g. using a different waveform model.

The script proposes chirp mass bounds in 2 steps:

1. Given the preliminary PE run, it checks of railing is present in the chirp mass posterior. This process involves a couple of free parameters that control how strict this check is. This is governed by the options --tolerance and nbins. See -h for more details.
2. If no railing is found, the proposed bounds are given by $[(1-f)\mathcal{M}_{\rm min},(1+f)\mathcal{M}_{\rm max}]$ where $f$ is a free parameter, given by the option --bounds_tol with the default value of 0.2.

Deciding the sampling rate

The sampling rate is determined by the requirement that it's high enough to resolve the highest frequency we are interested in. This is generally given by the frequency of the highest $(\ell,m)$ mode in the ringdown. This is estimated by using the SEOBNRv4HM routines. The user is given control which $\ell$ is used for this Nyquist frequency check via the option --ell-max. By default it is set to 3. To be precise, the sampling rate is computed as $2f_{\rm max}$ and then rounded up to the next power of 2.

Deciding the starting frequency

The starting frequency determines the length of the generated waveform. For SEOBNR time-domain family of waveforms there are two considerations that come into play:

1. The starting frequency (which corresponds to the frequency of the (2,2) mode) has to be such that the desired $(\ell,m)$ mode is present at the provided flow (which gives the lower bound of the likelihood integral and is an argument supplied by the user, set by default to 20 Hz). The script uses the PN relation (valid in the inspiral) that $f_{\ell m}=\frac{m}{2}f_{22}$ to dermine this. In principle, this consideration applies to any time-domain waveform.
2. The starting frequency is not too close to merger. This is a limitation of the SEOBNR time domain waveforms due to the need i) for a long enough inspiral to attach ringdown ii) to generate quasi-circular initial conditions. The script uses the same condition as the waveform model internally to check if the starting frequency is below this frequency limit.

Deciding the seglen

The raw segment length is determined by the following formula $t = T(1+s)+t_{\rm pad}+t_{\rm rollon}$, where $T$ is an estimate of the duration of the waveform, including inspiral, merger and ringdown; $s$ is a safety factor (default value of 0.03), $t_{\rm pad}$ is padding (default is 2 seconds) and $t_{\rm rollon}$ is the Tukey window roll on duration (default is 0.4 seconds). The duration $T$ is taken to be a function of the 2 component masses $m_{i}$, the z-components of the dimensionless spin $\chi_{i}$ and of course the low frequency cutoff flow and computed using the SEOBNRv4_ROM function. Notice that this means explicitly that we only care about the duration of the waveform from flow onwards and not from f_start. Finally the seglen is rounded to the next power of 2.

Checking the reference frequency

The reference frequency is an input parameter, encoded in the option --f_ref with a default value of 20 Hz. The script then computes the MECO frequency for all the samples, and checks that $f_{\rm ref}\leq 0.97f_{\rm MECO}$. If not, a message is printed out and even if specified a json report is not generated. Note that the script does not provide a recommendation for the reference frequency as this should be guided by physics.

Note that the failing of any railing checks and the reference frequency check can be ignored by passing the flag --override_safeties.