SNR Threshold Finder Example

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

np.random.seed(1234)

With injection catalog as hf5 file

# Example if you have the injection catalog file locally
# you can download the file for O4 injections from zenodo:
# !wget https://zenodo.org/records/16740117/files/samples-rpo4a_v2_20250503133839UTC-1366933504-23846400.hdf?download=1

# file_name = 'samples-rpo4a_v2_20250503133839UTC-1366933504-23846400.hdf'

# test = SNRThresholdFinder(
#     catalog_file = file_name,
#     # 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=None, # you can provide parameter values (np.ndarray) here if needed, if you don't want to use the catalog
#             ),
#     projected_detection_statistic = dict(
#                 key_name='observed_snr_net',
#                 threshold=None, # to be determined
#                 threshold_search_bounds=(6, 12),
#                 parameter=None, # you can provide parameter values (np.ndarray) here if needed, if you don't want to use the catalog
#             ),
#     parameters_to_fit = dict(
#         key_name = 'z',
#         parameter=None, # you can provide parameter values (np.ndarray) here if needed, if you don't want to use the catalog
#     ),
#     sample_size=20000,
#     selection_range = dict(
#         key_name = 'mass1_source',
#         range = (30, 60), # in solar masses
#         parameter = None, # you can provide parameter values (np.ndarray) here if needed, if you don't want to use the catalog
#     ),
#     multiprocessing_verbose=True,
# )

# best_thr, del_H, H, H_true, snr_thrs = test.find_threshold(iteration=10, print_output=True, no_multiprocessing=True)

With injection catalog parameters as numpy arrays

# injection_data_bbh.json is the reduced data extracted from the above hdf file for testing purpose
! rm injection_data_bbh.json # remove if already exists
! wget https://raw.githubusercontent.com/hemantaph/gwsnr/refs/heads/main/tests/unit/injection_data_bbh.json

--2026-07-03 14:24:52--  https://raw.githubusercontent.com/hemantaph/gwsnr/refs/heads/main/tests/unit/injection_data_bbh.json
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.111.133, 185.199.108.133, 185.199.109.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.111.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 10458366 (10.0M) [text/plain]
Saving to: ‘injection_data_bbh.json’

injection_data_bbh. 100%[===================>]   9.97M  11.0MB/s    in 0.9s

2026-07-03 14:24:54 (11.0 MB/s) - ‘injection_data_bbh.json’ saved [10458366/10458366]

from gwsnr.utils import get_param_from_json
params = get_param_from_json('injection_data_bbh.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=(6, 12),
                parameter=observed_snr_net,
            ),
    parameters_to_fit = dict(
        key_name = 'z',
        parameter=z,
    ),
    sample_size=20000,
    selection_range = dict(
        key_name = 'mass1_source',
        range = (30, 60), # in solar masses
        parameter = mass1_source,
    ),
)

best_thr, del_H, H, H_true, snr_thrs = test.find_threshold(iteration=10, print_output=True, no_multiprocessing=True)

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

Best SNR threshold: 9.21

from scipy.interpolate import UnivariateSpline

spline = UnivariateSpline(snr_thrs, del_H, s=0)
snr_fine = np.linspace(snr_thrs.min(), snr_thrs.max(), 100)
del_H_fine = spline(snr_fine)

import matplotlib.pyplot as plt
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "serif",
    "font.serif": ["Computer Modern Roman"],
})
plt.figure(figsize=(5, 2))
plt.plot(snr_fine, del_H_fine, 'b-', alpha=0.7, label='Spline interpolation')
plt.plot(snr_thrs, del_H, 'ro', markersize=6, label='Data points')
plt.axvline(best_thr, color='gray', linestyle='--', alpha=0.7, label=f'Best thr: {best_thr:.2f}')
plt.xlabel(r'$\rho_{\rm obs, net}$ Threshold (L1+H1)')
plt.ylabel(r'Cross-Entropy Difference ($\Delta H$)')
plt.grid(alpha=0.4)
plt.legend()
plt.tight_layout()
plt.savefig('snr_threshold_cross_entropy_difference.png', dpi=300, bbox_inches='tight')
plt.show()

Visualizing how well the detectable population matches when using SNRth and FAR cutoffs as detection criteria

snr_th = best_thr
plt.figure(figsize=(5, 2))
plt.hist(z[observed_snr_net>=snr_th], bins=50, density=True, alpha=0.5, color='C0', histtype='step', label=f'SNR greater than {best_thr:.2f}')
far_th = 1
plt.hist(z[gstlal_far<far_th], bins=50, density=True, alpha=0.5, color='C1', histtype='step', label=f'FAR less than {far_th:.2f}')
plt.hist(z, bins=100, density=True, alpha=0.5, color='C2', histtype='step', label='All injections')
plt.xlabel('Redshift z')
plt.xlim(0, 1.5)
plt.ylabel('Probability Density')
plt.legend()
plt.grid(alpha=0.4)
plt.tight_layout()
plt.savefig('redshift_distribution.png', dpi=300, bbox_inches='tight')
plt.show()

# # check rates
# with h5py.File(file_name, 'r') as obj:
#     attrs = dict(obj.attrs.items())
#     events = obj['events'][:]

# observed_snr_net = events['observed_snr_net']
# gstlal_far = events['gstlal_far']

# fraction_above_snr_thresh = np.sum(observed_snr_net >= best_thr) / len(observed_snr_net)
# fraction_above_gstlal_far = np.sum(gstlal_far < 1) / len(observed_snr_net)

# print(f"Fraction above SNR threshold {best_thr:.2f}: {fraction_above_snr_thresh:.6f}")
# print(f"Fraction above GstLAL FAR threshold 1 per year: {fraction_above_gstlal_far:.6f}")