Horizon Distance Calculation with Mass-Dependent SNR Thresholds

Find the threshold for BNS systems

• choose ‘mass1_source’ \(\in [1,3]\)

• we will use pre-selected injection catalogue data.

• I will take threshold for individual detector same as that of detector network, for simplicity.

• default sensitivity is consider with O4 psds.

# original catalogue data: https://zenodo.org/records/16740117/files/samples-rpo4a_v2_20250503133839UTC-1366933504-23846400.hdf?download=1
# injection_data_bns.json is the reduced data extracted from the above hdf file for testing purpose
! rm injection_data_bns.json # remove if exists
! wget https://raw.githubusercontent.com/hemantaph/gwsnr/refs/heads/main/tests/unit/injection_data_bns.json

--2026-01-02 15:35:22--  https://raw.githubusercontent.com/hemantaph/gwsnr/refs/heads/main/tests/unit/injection_data_bns.json
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.111.133, 185.199.109.133, 185.199.110.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.111.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 6444606 (6.1M) [text/plain]
Saving to: ‘injection_data_bns.json’

injection_data_bns. 100%[===================>]   6.15M  7.97MB/s    in 0.8s

2026-01-02 15:35:23 (7.97 MB/s) - ‘injection_data_bns.json’ saved [6444606/6444606]

import h5py
import numpy as np
import matplotlib.pyplot as plt
from gwsnr.threshold import SNRThresholdFinder

from gwsnr.utils import get_param_from_json
params = get_param_from_json('injection_data_bns.json')

gstlal_far = params['gstlal_far']
observed_snr_net = params['observed_snr_net']
z = params['z']
mass1_source = params['mass1_source']

test = SNRThresholdFinder(
    catalog_file = None,
    # below are all default values. You can omit them if you want.
    npool=4,
    original_detection_statistic = dict(
                key_name='gstlal_far',
                threshold=1,  # 1 per year
                parameter=gstlal_far,
            ),
    projected_detection_statistic = dict(
                key_name='observed_snr_net',
                threshold=None, # to be determined
                threshold_search_bounds=(8, 14),
                parameter=observed_snr_net,
            ),
    parameters_to_fit = dict(
        key_name = 'z',
        parameter=z,
    ),
    sample_size=20000,
    selection_range = dict(
        key_name = 'mass1_source',
        range = (1, 3), # in solar masses
        parameter = mass1_source,
    ),
)

best_thr, _, _, _, _ = test.find_threshold(iteration=10, print_output=True, no_multiprocessing=True)

100%|███████████████████████████████████████████████████████████████| 10/10 [00:10<00:00,  1.01s/it]

Best SNR threshold: 11.71

snr_th, snr_th_net = best_thr, best_thr

from gwsnr import GWSNR

mtot_min = (1+1)*(1+0) # (mass_1_min + mass_2_min)*(1 + z_min)
mtot_max = (2.3+2.3)*(1+5) # (mass_1_max + mass_2_max)*(1 + z_max)

gwsnr = GWSNR(npool=6, snr_method="interpolation_no_spins", mtot_min=mtot_min, mtot_max=mtot_max)

Initializing GWSNR class...

psds not given. Choosing bilby's default psds
Interpolator will be loaded for L1 detector from ./interpolator_json/L1/partialSNR_dict_4.json
Interpolator will be loaded for H1 detector from ./interpolator_json/H1/partialSNR_dict_4.json
Interpolator will be loaded for V1 detector from ./interpolator_json/V1/partialSNR_dict_4.json

Chosen GWSNR initialization parameters:

npool:  6
snr type:  interpolation_no_spins
waveform approximant:  IMRPhenomD
sampling frequency:  2048.0
minimum frequency (fmin):  20.0
reference frequency (f_ref):  20.0
mtot=mass1+mass2
min(mtot):  2
max(mtot) (with the given fmin=20.0): 27.599999999999998
detectors:  ['L1', 'H1', 'V1']
psds:  [[array([  10.21659,   10.23975,   10.26296, ..., 4972.81   ,
       4984.081  , 4995.378  ]), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,
       6.51153524e-46, 6.43165104e-46, 6.55252996e-46]), <scipy.interpolate._interpolate.interp1d object at 0x154cdf740>], [array([  10.21659,   10.23975,   10.26296, ..., 4972.81   ,
       4984.081  , 4995.378  ]), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,
       6.51153524e-46, 6.43165104e-46, 6.55252996e-46]), <scipy.interpolate._interpolate.interp1d object at 0x154cdfd30>], [array([   10.      ,    10.02306 ,    10.046173, ...,
        9954.0389  ,  9976.993   , 10000.      ]), array([1.22674387e-42, 1.20400299e-42, 1.18169466e-42, ...,
       1.51304203e-43, 1.52010157e-43, 1.52719372e-43]), <scipy.interpolate._interpolate.interp1d object at 0x154ed0130>]]

Horizon Distance - Analytic

hd_dict_ =  gwsnr.horizon_distance_analytical(mass_1=1.4, mass_2=1.4, snr_th=snr_th)
hd_dict_ # in Mpc

{'horizon_distance_L1': array([284.52248436]),
 'horizon_distance_H1': array([284.52248436]),
 'horizon_distance_V1': array([217.19546115])}

Horizon Distance - Numerical

hd_dict_, optimal_sky_dict =  gwsnr.horizon_distance_numerical(mass_1=1.4, mass_2=1.4, snr_th=snr_th)

hd_dict_ # in Mpc

{'horizon_distance_L1': 333.1959886579425,
 'horizon_distance_H1': 333.1959886579425,
 'horizon_distance_V1': 254.35127570125042,
 'horizon_distance_net': 469.8682819388923}

# optimal sky position
optimal_sky_dict

{'optimal_sky_location_L1': (5.893069985246483, 0.5327305795875431),
 'optimal_sky_location_H1': (2.2513256071396226, -0.81130311322348),
 'optimal_sky_location_V1': (4.51922945526775, -0.7614903798992666),
 'optimal_sky_location_net': (2.6487926231001, -0.6781071251407678)}