#
# =============================================================================
#
# Preamble
#
# =============================================================================
#
""" This module contains functions for calculating expected rates of noise
and signal coincidences.
"""
import itertools
import logging
import numpy
import pycbc.detector
logger = logging.getLogger('pycbc.events.coinc_rate')
[docs]def multiifo_noise_lograte(log_rates, slop):
"""
Calculate the expected rate of noise coincidences for multiple
combinations of detectors
Parameters
----------
log_rates: dict
Key: ifo string, Value: sequence of log single-detector trigger rates,
units assumed to be Hz
slop: float
time added to maximum time-of-flight between detectors to account
for timing error
Returns
-------
expected_log_rates: dict
Key: ifo combination string
Value: expected log coincidence rate in the combination, units log Hz
"""
expected_log_rates = {}
# Order of ifos must be stable in output dict keys, so sort them
ifos = sorted(list(log_rates.keys()))
ifostring = ' '.join(ifos)
# Calculate coincidence for all-ifo combination
expected_log_rates[ifostring] = \
combination_noise_lograte(log_rates, slop)
# If more than one possible coincidence type exists,
# calculate coincidence for subsets through recursion
if len(ifos) > 2:
# Calculate rate for each 'miss-one-out' detector combination
subsets = itertools.combinations(ifos, len(ifos) - 1)
for subset in subsets:
rates_subset = {}
for ifo in subset:
rates_subset[ifo] = log_rates[ifo]
sub_coinc_rates = multiifo_noise_lograte(rates_subset, slop)
# add these sub-coincidences to the overall dictionary
for sub_coinc in sub_coinc_rates:
expected_log_rates[sub_coinc] = sub_coinc_rates[sub_coinc]
return expected_log_rates
[docs]def combination_noise_rate(rates, slop):
"""
Calculate the expected rate of noise coincidences for a combination of
detectors
WARNING: for high stat values, numerical underflow can occur
Parameters
----------
rates: dict
Key: ifo string, Value: sequence of single-detector trigger rates,
units assumed to be Hz
slop: float
time added to maximum time-of-flight between detectors to account
for timing error
Returns
-------
numpy array
Expected coincidence rate in the combination, units Hz
"""
logger.warning('combination_noise_rate() is liable to numerical '
'underflows, use combination_noise_lograte '
'instead')
log_rates = {k: numpy.log(r) for (k, r) in rates.items()}
# exp may underflow
return numpy.exp(combination_noise_lograte(log_rates, slop))
[docs]def combination_noise_lograte(log_rates, slop, dets=None):
"""
Calculate the expected rate of noise coincidences for a combination of
detectors given log of single detector noise rates
Parameters
----------
log_rates: dict
Key: ifo string, Value: sequence of log single-detector trigger rates,
units assumed to be Hz
slop: float
time added to maximum time-of-flight between detectors to account
for timing error
dets: dict
Key: ifo string, Value: pycbc.detector.Detector object
If not provided, will be created from ifo strings
Returns
-------
numpy array
Expected log coincidence rate in the combination, units Hz
"""
# multiply product of trigger rates by the overlap time
allowed_area = multiifo_noise_coincident_area(list(log_rates),
slop,
dets=dets)
# list(dict.values()) is python-3-proof
rateprod = numpy.sum(list(log_rates.values()), axis=0)
return numpy.log(allowed_area) + rateprod
[docs]def multiifo_noise_coincident_area(ifos, slop, dets=None):
"""
Calculate the total extent of time offset between 2 detectors,
or area of the 2d space of time offsets for 3 detectors, for
which a coincidence can be generated
Cannot yet handle more than 3 detectors.
Parameters
----------
ifos: list of strings
list of interferometers
slop: float
extra time to add to maximum time-of-flight for timing error
dets: dict
Key: ifo string, Value: pycbc.detector.Detector object
If not provided, will be created from ifo strings
Returns
-------
allowed_area: float
area in units of seconds^(n_ifos-1) that coincident values can fall in
"""
# set up detector objects
if dets is None:
dets = {ifo: pycbc.detector.Detector(ifo) for ifo in ifos}
n_ifos = len(ifos)
if n_ifos == 2:
allowed_area = 2. * \
(dets[ifos[0]].light_travel_time_to_detector(dets[ifos[1]]) + slop)
elif n_ifos == 3:
tofs = numpy.zeros(n_ifos)
ifo2_num = []
# calculate travel time between detectors (plus extra for timing error)
# TO DO: allow for different timing errors between different detectors
for i, ifo in enumerate(ifos):
ifo2_num.append(int(numpy.mod(i + 1, n_ifos)))
det0 = dets[ifo]
det1 = dets[ifos[ifo2_num[i]]]
tofs[i] = det0.light_travel_time_to_detector(det1) + slop
# combine these to calculate allowed area
allowed_area = 0
for i, _ in enumerate(ifos):
allowed_area += 2 * tofs[i] * tofs[ifo2_num[i]] - tofs[i]**2
else:
raise NotImplementedError("Not able to deal with more than 3 ifos")
return allowed_area
[docs]def multiifo_signal_coincident_area(ifos):
"""
Calculate the area in which signal time differences are physically allowed
Parameters
----------
ifos: list of strings
list of interferometers
Returns
-------
allowed_area: float
area in units of seconds^(n_ifos-1) that coincident signals will occupy
"""
n_ifos = len(ifos)
if n_ifos == 2:
det0 = pycbc.detector.Detector(ifos[0])
det1 = pycbc.detector.Detector(ifos[1])
allowed_area = 2 * det0.light_travel_time_to_detector(det1)
elif n_ifos == 3:
dets = {}
tofs = numpy.zeros(n_ifos)
ifo2_num = []
# set up detector objects
for ifo in ifos:
dets[ifo] = pycbc.detector.Detector(ifo)
# calculate travel time between detectors
for i, ifo in enumerate(ifos):
ifo2_num.append(int(numpy.mod(i + 1, n_ifos)))
det0 = dets[ifo]
det1 = dets[ifos[ifo2_num[i]]]
tofs[i] = det0.light_travel_time_to_detector(det1)
# calculate allowed area
phi_12 = numpy.arccos((tofs[0]**2 + tofs[1]**2 - tofs[2]**2)
/ (2 * tofs[0] * tofs[1]))
allowed_area = numpy.pi * tofs[0] * tofs[1] * numpy.sin(phi_12)
else:
raise NotImplementedError("Not able to deal with more than 3 ifos")
return allowed_area