# Module with utilities for estimating candidate events source probabilities
# Initial code by A. Curiel Barroso, August 2019
# Modified by V. Villa-Ortega, January 2020
"""Functions to compute the area corresponding to different CBC on the m1 & m2
plane when given a central mchirp value and uncertainty.
It also includes a function that calculates the source frame when given the
detector frame mass and redshift.
"""
import math
from pycbc.conversions import mass2_from_mchirp_mass1 as m2mcm1
from scipy.integrate import quad
from pycbc.cosmology import _redshift
[docs]def insert_args(parser):
mchirp_group = parser.add_argument_group("Arguments for estimating the "
"source probabilities of a "
"candidate event using the snr, "
"mchirp, and effective distance.")
mchirp_group.add_argument('--src-class-mass-range', type=float, nargs=2,
metavar=('MIN_M2', 'MAX_M1'),
default=[1.0, 45.0],
help="Minimum and maximum values for the mass "
"of the binary components, used as limits "
"of the mass plane when computing the area "
"corresponding to different CBC sources.")
mchirp_group.add_argument('--src-class-mass-gap', type=float, nargs=2,
metavar=('MAX_NS', 'MIN_BH'), default=[3.0, 5.0],
help="Limits of the mass gap, that correspond "
"to the maximum mass of a neutron star "
"and the minimum mass of a black hole. "
"Used as limits of integration of the "
"different CBC regions.")
mchirp_group.add_argument('--src-class-mchirp-to-delta', type=float,
metavar='m0', default=0.01,
help='Coefficient to estimate the value of the '
'mchirp uncertainty by mchirp_delta = '
'm0 * mchirp.')
mchirp_group.add_argument('--src-class-eff-to-lum-distance', type=float,
metavar='a0', default=0.759,
help='Coefficient to estimate the value of the '
'luminosity distance from the minimum '
'eff distance by D_lum = a0 * min(D_eff).')
mchirp_group.add_argument('--src-class-lum-distance-to-delta', type=float,
nargs=2, metavar=('b0', 'b1'),
default=[-0.449, -0.342],
help='Coefficients to estimate the value of the '
'uncertainty on the luminosity distance '
'from the estimated luminosity distance and'
' the coinc snr by delta_lum = D_lum * '
'exp(b0) * coinc_snr ** b1.')
mchirp_group.add_argument('--src-class-mass-gap-separate',
action='store_true',
help='Gives separate probabilities for each kind'
' of mass gap CBC sources: GNS, GG, BHG.')
[docs]def from_cli(args):
return {'mass_limits': {'max_m1': args.src_class_mass_range[1],
'min_m2': args.src_class_mass_range[0]},
'mass_bdary': {'ns_max': args.src_class_mass_gap[0],
'gap_max': args.src_class_mass_gap[1]},
'estimation_coeff': {'a0': args.src_class_eff_to_lum_distance,
'b0': args.src_class_lum_distance_to_delta[0],
'b1': args.src_class_lum_distance_to_delta[1],
'm0': args.src_class_mchirp_to_delta},
'mass_gap': args.src_class_mass_gap_separate}
[docs]def src_mass_from_z_det_mass(z, del_z, mdet, del_mdet):
"""Takes values of redshift, redshift uncertainty, detector mass and its
uncertainty and computes the source mass and its uncertainty.
"""
msrc = mdet / (1. + z)
del_msrc = msrc * ((del_mdet / mdet) ** 2.
+ (del_z / (1. + z)) ** 2.) ** 0.5
return (msrc, del_msrc)
[docs]def intmc(mc, x_min, x_max):
"""Returns the integral of a component mass as a function of the mass of
the other component, taking mchirp as an argument.
"""
integral = quad(lambda x, mc: m2mcm1(mc, x), x_min, x_max, args=mc)
return integral[0]
[docs]def calc_areas(trig_mc_det, mass_limits, mass_bdary, z, mass_gap):
"""Computes the area inside the lines of the second component mass as a
function of the first component mass for the two extreme values
of mchirp: mchirp +/- mchirp_uncertainty, for each region of the source
classifying diagram.
"""
trig_mc = src_mass_from_z_det_mass(z["central"], z["delta"],
trig_mc_det["central"],
trig_mc_det["delta"])
mcb = trig_mc[0] + trig_mc[1]
mcs = trig_mc[0] - trig_mc[1]
m2_min = mass_limits["min_m2"]
m1_max = mass_limits["max_m1"]
ns_max = mass_bdary["ns_max"]
gap_max = mass_bdary["gap_max"]
# The points where the equal mass line and a chirp mass
# curve intersect is m1 = m2 = 2**0.2 * mchirp
mib = (2.**0.2) * mcb
mis = (2.**0.2) * mcs
# AREA FOR BBH
if mib < gap_max:
abbh = 0.0
else:
limb_bbh = min(m1_max, m2mcm1(mcb, gap_max))
intb_bbh = intmc(mcb, mib, limb_bbh)
if mis < gap_max:
lims1_bbh = gap_max
lims2_bbh = lims1_bbh
else:
lims1_bbh = mis
lims2_bbh = min(m1_max, m2mcm1(mcs, gap_max))
ints_bbh = intmc(mcs, lims1_bbh, lims2_bbh)
limdiag_bbh = max(m2mcm1(mcs, lims1_bbh), gap_max)
intline_sup_bbh = 0.5 * (limdiag_bbh + mib) * (mib - lims1_bbh)
intline_inf_bbh = (limb_bbh - lims2_bbh) * gap_max
int_sup_bbh = intb_bbh + intline_sup_bbh
int_inf_bbh = ints_bbh + intline_inf_bbh
abbh = int_sup_bbh - int_inf_bbh
# AREA FOR BHG
if m2mcm1(mcb, gap_max) < ns_max or m2mcm1(mcs, m1_max) > gap_max:
abhg = 0.0
else:
if m2mcm1(mcb, m1_max) > gap_max:
limb2_bhg = m1_max
limb1_bhg = limb2_bhg
else:
limb2_bhg = min(m1_max, m2mcm1(mcb, ns_max))
limb1_bhg = max(gap_max, m2mcm1(mcb, gap_max))
intb_bhg = intmc(mcb, limb1_bhg, limb2_bhg)
if m2mcm1(mcs, gap_max) < ns_max:
lims2_bhg = gap_max
lims1_bhg = lims2_bhg
else:
lims1_bhg = max(gap_max, m2mcm1(mcs, gap_max))
lims2_bhg = min(m1_max, m2mcm1(mcs, ns_max))
intline_inf_bhg = (limb2_bhg - lims2_bhg) * ns_max
intline_sup_bhg = (limb1_bhg - lims1_bhg) * gap_max
ints_bhg = intmc(mcs, lims1_bhg, lims2_bhg)
int_sup_bhg = intb_bhg + intline_sup_bhg
int_inf_bhg = ints_bhg + intline_inf_bhg
abhg = int_sup_bhg - int_inf_bhg
# AREA FOR GG
if m2mcm1(mcs, gap_max) > gap_max or m2mcm1(mcb, ns_max) < ns_max:
agg = 0.0
else:
if m2mcm1(mcb, gap_max) > gap_max:
limb2_gg = gap_max
limb1_gg = limb2_gg
else:
limb1_gg = mib
limb2_gg = min(gap_max, m2mcm1(mcb, ns_max))
intb_gg = intmc(mcb, limb1_gg, limb2_gg)
if m2mcm1(mcs, ns_max) < ns_max:
lims2_gg = ns_max
lims1_gg = lims2_gg
else:
lims1_gg = mis
lims2_gg = min(gap_max, m2mcm1(mcs, ns_max))
ints_gg = intmc(mcs, lims1_gg, lims2_gg)
limdiag1_gg = max(m2mcm1(mcs, lims1_gg), ns_max)
limdiag2_gg = min(m2mcm1(mcb, limb1_gg), gap_max)
intline_sup_gg = (0.5 * (limb1_gg - lims1_gg)
* (limdiag1_gg + limdiag2_gg))
intline_inf_gg = (limb2_gg - lims2_gg) * ns_max
int_sup_gg = intb_gg + intline_sup_gg
int_inf_gg = ints_gg + intline_inf_gg
agg = int_sup_gg - int_inf_gg
# AREA FOR BNS
if m2mcm1(mcs, ns_max) > ns_max:
abns = 0.0
else:
if m2mcm1(mcb, ns_max) > ns_max:
limb2_bns = ns_max
limb1_bns = limb2_bns
else:
limb2_bns = min(ns_max, m2mcm1(mcb, m2_min))
limb1_bns = mib
intb_bns = intmc(mcb, limb1_bns, limb2_bns)
if mis < m2_min:
lims2_bns = m2_min
lims1_bns = lims2_bns
else:
lims2_bns = min(ns_max, m2mcm1(mcs, m2_min))
lims1_bns = mis
ints_bns = intmc(mcs, lims1_bns, lims2_bns)
intline_inf_bns = (limb2_bns - lims2_bns) * m2_min
limdiag1_bns = max(m2mcm1(mcs, lims1_bns), m2_min)
limdiag2_bns = min(m2mcm1(mcb, limb1_bns), ns_max)
intline_sup_bns = (0.5 * (limdiag1_bns + limdiag2_bns)
* (limb1_bns - lims1_bns))
int_sup_bns = intb_bns + intline_sup_bns
int_inf_bns = ints_bns + intline_inf_bns
abns = int_sup_bns - int_inf_bns
# AREA FOR GNS
if m2mcm1(mcs, gap_max) > ns_max or m2mcm1(mcb, ns_max) < m2_min:
agns = 0.0
else:
if m2mcm1(mcb, gap_max) > ns_max:
limb2_gns = gap_max
limb1_gns = limb2_gns
else:
limb2_gns = min(gap_max, m2mcm1(mcb, m2_min))
limb1_gns = max(ns_max, m2mcm1(mcb, ns_max))
intb_gns = intmc(mcb, limb1_gns, limb2_gns)
if m2mcm1(mcs, ns_max) < m2_min:
lims2_gns = ns_max
lims1_gns = lims2_gns
else:
lims1_gns = max(ns_max, m2mcm1(mcs, ns_max))
lims2_gns = min(gap_max, m2mcm1(mcs, m2_min))
intline_inf_gns = (limb2_gns - lims2_gns) * m2_min
intline_sup_gns = (limb1_gns - lims1_gns) * ns_max
ints_gns = intmc(mcs, lims1_gns, lims2_gns)
int_sup_gns = intb_gns + intline_sup_gns
int_inf_gns = ints_gns + intline_inf_gns
agns = int_sup_gns - int_inf_gns
# AREA FOR NSBH
if m2mcm1(mcs, m1_max) > ns_max or m2mcm1(mcb, gap_max) < m2_min:
ansbh = 0.0
else:
if m2mcm1(mcb, m1_max) > ns_max:
limb2_nsbh = m1_max
limb1_nsbh = limb2_nsbh
else:
limb1_nsbh = max(gap_max, m2mcm1(mcb, ns_max))
limb2_nsbh = min(m1_max, m2mcm1(mcb, m2_min))
intb_nsbh = intmc(mcb, limb1_nsbh, limb2_nsbh)
if m2mcm1(mcs, gap_max) < m2_min:
lims1_nsbh = gap_max
lims2_nsbh = lims1_nsbh
else:
lims1_nsbh = max(gap_max, m2mcm1(mcs, ns_max))
lims2_nsbh = min(m1_max, m2mcm1(mcs, m2_min))
intline_inf_nsbh = (limb2_nsbh - lims2_nsbh) * m2_min
intline_sup_nsbh = (limb1_nsbh - lims1_nsbh) * ns_max
ints_nsbh = intmc(mcs, lims1_nsbh, lims2_nsbh)
int_sup_nsbh = intb_nsbh + intline_sup_nsbh
int_inf_nsbh = ints_nsbh + intline_inf_nsbh
ansbh = int_sup_nsbh - int_inf_nsbh
if mass_gap:
return {
"BNS": abns,
"GNS": agns,
"NSBH": ansbh,
"GG": agg,
"BHG": abhg,
"BBH": abbh
}
return {
"BNS": abns,
"NSBH": ansbh,
"BBH": abbh,
"Mass Gap": agns + agg + abhg
}
[docs]def calc_probabilities(mchirp, snr, eff_distance, src_args):
"""Computes the different probabilities that a candidate event belongs to
each CBC source category taking as arguments the chirp mass, the
coincident SNR and the effective distance, and estimating the
chirp mass uncertainty, the luminosity distance (and its uncertainty)
and the redshift (and its uncertainty). Probability estimation is done
assuming it is directly proportional to the area laying in the
correspondent CBC region.
"""
mass_limits = src_args['mass_limits']
mass_bdary = src_args['mass_bdary']
coeff = src_args['estimation_coeff']
trig_mc_det = {'central': mchirp, 'delta': mchirp * coeff['m0']}
dist_estimation = coeff['a0'] * eff_distance
dist_std_estimation = (dist_estimation * math.exp(coeff['b0']) *
snr ** coeff['b1'])
z_estimation = _redshift(dist_estimation)
z_est_max = _redshift(dist_estimation + dist_std_estimation)
z_est_min = _redshift(dist_estimation - dist_std_estimation)
z_std_estimation = 0.5 * (z_est_max - z_est_min)
z = {'central': z_estimation, 'delta': z_std_estimation}
mass_gap = src_args['mass_gap']
# If the mchirp is greater than the mchirp corresponding to two masses
# equal to the maximum mass, the probability for BBH is 100%
mc_max = mass_limits['max_m1'] / (2 ** 0.2)
if trig_mc_det['central'] > mc_max * (1 + z['central']):
if mass_gap:
probabilities = {"BNS": 0.0, "GNS": 0.0, "NSBH": 0.0, "GG": 0.0,
"BHG": 0.0, "BBH": 1.0}
else:
probabilities = {"BNS": 0.0, "NSBH": 0.0, "BBH": 1.0,
"Mass Gap": 0.0}
else:
areas = calc_areas(trig_mc_det, mass_limits, mass_bdary, z, mass_gap)
total_area = sum(areas.values())
probabilities = {key: areas[key]/total_area for key in areas}
return probabilities