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.

[7]:

## original catalogue data: https://zenodo.org/records/16740117/files/samples-rpo4a_v2_20250503133839UTC-1366933504-23846400.hdf?download=1

## Uncomment the following lines to download the original catalogue data for testing purpose
# injection_data_bns.json is the reduced data extracted from the above hdf file for testing purpose
# it contains L1 and H1 injection data
! 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-07-08 23:23:52--  https://raw.githubusercontent.com/hemantaph/gwsnr/refs/heads/main/tests/unit/injection_data_bns.json
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.110.133, 185.199.108.133, 185.199.111.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.110.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  8.12MB/s    in 0.8s

2026-07-08 23:23:54 (8.12 MB/s) - ‘injection_data_bns.json’ saved [6444606/6444606]

[16]:

import h5py
import numpy as np
np.random.seed(42)
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:14<00:00,  1.43s/it]

Best SNR threshold: 11.62

[17]:

snr_th, snr_th_net = best_thr, best_thr
print(f"Best threshold for SNR: {snr_th:.2f}, Best threshold for SNR_net: {snr_th_net:.2f}")

Best threshold for SNR: 11.62, Best threshold for SNR_net: 11.62

[18]:

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,
    ifos=['H1', 'L1'],
    snr_method="interpolation_no_spins",
    mtot_min=mtot_min,
    mtot_max=mtot_max,
)


Initializing GWSNR class...

Interpolator will be loaded for H1 detector from ./interpolator_json/H1/partialSNR_dict_3.json
Interpolator will be loaded for L1 detector from ./interpolator_json/L1/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:  ['H1', 'L1']
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 0x30cf05bc0>], [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 0x30fd9fb50>]]

Horizon Distance - Analytic

[19]:

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

[19]:

{'horizon_distance_H1': array([286.83332778]),
 'horizon_distance_L1': array([286.83332778])}

Horizon Distance - Numerical

[20]:

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

hd_dict_ # in Mpc

[20]:

{'horizon_distance_H1': 286.8333285966073,
 'horizon_distance_L1': 286.8333285966073,
 'horizon_distance_net': 457.9805338138831}

[21]:

# optimal orientation
optimal_orientation

[21]:

{'optimal_orientation_H1': {'ra': 2.251317689687112,
  'dec': -0.8112775059498398,
  'theta_jn': 3.141592653589793},
 'optimal_orientation_L1': {'ra': 5.893058202133175,
  'dec': 0.5327335489575311,
  'theta_jn': 0.0},
 'optimal_orientation_net': {'ra': 2.5299508320516666,
  'dec': -0.6873371132879433,
  'theta_jn': 3.435784008143372e-05}}

[ ]:

[ ]: