GitLab

cc: @sylvia.biscoveanu

When running a matched-filter search on an event multiple times with Gaussian noise, we expect to see each interferometer's SNR distributed as a Gaussian about the optimal SNR with unit variance (see fig. 3 in [1]). For a network of N (in this case, four) detectors, the resulting network SNR^2 should be distributed as a noncentral chi-square distribution with N degrees of freedom and the non-centrality parameter equal to the network optimal SNR^2. (see [2])

However, this isn't the distribution that I'm seeing, for 100 matched-filter searches of the same event using bayestar-realize-coincs:

Here is code to replicate this example. The exact files can also be found on the CIT cluster at /home/geoffrey.mo/bayestar-test/noise_snr. My ligo.skymap version is 0.1.16.

lalapps_inspinj -o inj.xml --m-distr fixMasses --fixed-mass1 1.3 --fixed-mass2 1.3 --t-distr uniform --time-step 100 --gps-start-time 1000000000 --gps-end-time 1000000099 --d-distr volume --min-distance 150000 --max-distance 160000 --l-distr random --i-distr uniform --f-lower 20 --disable-spin --waveform IMRPhenom_NRTidal

# Make a copy of this, call it inj_many.xml.
cp inj.xml inj_many.xml

# DO THIS MANUALLY # 
# Then copy the line in the sim_inspiral inj_many.xml table many times (in this case, 100).
# Make sure to add the delimiter at the end of the line before copying.

bayestar-sample-model-psd -o psd.xml --H1=aLIGOAdVO4T1800545 --L1=aLIGOAdVO4T1800545 --V1=AdVMidLowSensitivityP1200087 --K1=KAGRAEarlySensitivityT1600593

# Get the zero-noise (optimal) SNR:
bayestar-realize-coincs -o coinc.xml inj.xml --reference-psd psd.xml --detector H1 L1 V1 K1 --snr-threshold 0.0 --waveform IMRPhenomD_NRTidal --min-triggers 1 --net-snr-threshold 1 -j 64

# Now localize the many injs:
bayestar-realize-coincs -o coinc_many.xml inj_many.xml --reference-psd psd.xml --detector H1 L1 V1 K1 --snr-threshold 0.0 --waveform IMRPhenomD_NRTidal --min-triggers 1 --net-snr-threshold 1 -j 64 --measurement-error gaussian-noise

# Run this python script to parse the xmls for the SNRs and plot the distribution
python parse_plot.py

And this is parse_plot.py:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import xml.etree.ElementTree as ET

# These don't seem to work due to some ligolw incompatabilities with the
# coinc_event_id column validity. I really don't want to debug this right now
# so I'm just going to do some ugly parsing of the xml to get the SNRs.
# coinc_one = EventTable.read('coinc.xml', tablename='coinc_inspiral')
# coinc_many = EventTable.read('coinc_many.xml', tablename='coinc_inspiral')

single_tree = ET.parse('coinc.xml')
single_root = single_tree.getroot()
optimal_snr = float(single_root[1][-1].text.split(',')[-1].split('\n')[0])

many_tree = ET.parse('coinc_many.xml')
many_root = many_tree.getroot()
events = [event.split(',') for event in many_root[1][-1].text.split('\n')[:-1] if
          len(event) > 1]
events[-1].append('') # last row needs fixing to look like the other rows
snrs = np.array([float(event[-2]) for event in events])

xarray = np.linspace(5, 15, 200)

plt.figure(figsize=(12,8))
# this is the distribution we expect for the quadrature sum
# of four Gaussians (one for each detector)
plt.hist(np.random.noncentral_chisquare(4, optimal_snr**2, 100000), 
    bins=200, histtype='step', density=True, label='Noncentral chi squared')
plt.hist(snrs**2, bins=20, histtype='step', density=True, label='Observed');
plt.axvline(optimal_snr**2, ls='--', color='black', label=r'$\rho_\mathrm{opt}$');
plt.legend(loc='upper left');
plt.xlabel('SNR$^2$');
plt.ylabel('Probability density');
plt.savefig('snr2_dist.png', dpi=300)

[1] https://arxiv.org/pdf/1809.02293.pdf

[2] Mathematica code:

In[1] = TransformedDistribution[
 s^2 + t^2 + u^2 + v^2, {s \[Distributed] NormalDistribution[w, 1], 
  t \[Distributed] NormalDistribution[x, 1], 
  u \[Distributed] NormalDistribution[y, 1], 
  v \[Distributed] NormalDistribution[z, 1]}]

Out[1] = NoncentralChiSquareDistribution[4, w^2 + x^2 + y^2 + z^2]