From 1392f6696b9e7ea777e97620a75ee4178e2f5f34 Mon Sep 17 00:00:00 2001
From: Ryan Michael Magee <ryan.magee@ldas-pcdev6.ligo.caltech.edu>
Date: Fri, 21 Jun 2019 10:24:57 -0700
Subject: [PATCH] Change cov matrix calculation to use actual waveform
amplitude/parameters instead of stationary phase approx
---
...lal_inspiral_compute_dtdphideff_cov_matrix | 22 +++++++++--------
...inspiral_create_dt_dphi_snr_ratio_pdfs_dag | 6 ++++-
gstlal-inspiral/python/templates.py | 24 +++++++++++--------
3 files changed, 31 insertions(+), 21 deletions(-)
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix b/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix
index 57e18a23e1..57f0cbcb9e 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix
+++ b/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix
@@ -20,6 +20,7 @@ import itertools
import numpy, scipy.interpolate
from lal import series
from gstlal.stats.inspiral_extrinsics import TimePhaseSNR
+from gstlal import templates
from ligo.lw import utils as ligolw_utils
import sys
import h5py
@@ -29,6 +30,10 @@ parser.add_argument('--psd-xml', help = 'XML containing HLV psd')
parser.add_argument('--H-snr', type = float, help = 'H characteristic SNR')
parser.add_argument('--L-snr', type = float, help = 'L characteristic SNR')
parser.add_argument('--V-snr', type = float, help = 'V characteristic SNR')
+parser.add_argument('--m1', type = float, default = 1.4, help = 'primary component mass')
+parser.add_argument('--m2', type = float, default = 1.4, help = 'secondary component mass')
+parser.add_argument('--s1', type = float, default = 0., help = 'primary (z) spin')
+parser.add_argument('--s2', type = float, default = 0., help = 'secondary (z) spin')
parser.add_argument('--flow', type = float, default = 10., help = 'Low frequency cut-off. Default 10 Hz')
parser.add_argument('--fhigh', type = float, default = 1024., help = 'High frequency cut-off. Default 1024 Hz')
parser.add_argument("--output", help = "set the output h5 file, e.g., covmat.h5")
@@ -47,13 +52,10 @@ for ifo in psd:
f = numpy.arange(len(psd[ifo].data.data)) * psd[ifo].deltaF
psddict[ifo] = scipy.interpolate.interp1d(f, psd[ifo].data.data)
-def fn(fmin, fmax, n, ifo, psddict = psddict):
- farr = numpy.linspace(fmin, fmax, 100000)
- df = farr[1] - farr[0]
- return 4 * numpy.sum(farr**(-7/3.) * farr**n * df / psddict[ifo](farr))
-
-def moment(fmin, fmax, n, ifo, psddict = psddict):
- return fn(fmin, fmax, n, ifo, psddict) / fn(fmin, fmax, 0, ifo, psddict)
+def moment(fmin, fmax, n, m1, m2, s1, s2, ifo, psddict = psddict, delta_f = 0.25):
+ farr = numpy.linspace(0, fmax, fmax / delta_f + delta_f)
+ h = templates.hplus_of_f(m1, m2, s1, s2, fmin, fmax, delta_f)
+ return templates.moment(farr, h, n, psddict[ifo])
sigsqtt = {}
sigsqpp = {}
@@ -61,10 +63,10 @@ sigsqtp = {}
sigsqdd = {}
for ifo in rho:
- sigsqf = moment(args.flow, args.fhigh, 2, ifo) - moment(args.flow, args.fhigh, 1, ifo)**2
+ sigsqf = moment(args.flow, args.fhigh, 2, args.m1, args.m2, args.s1, args.s2, ifo) - moment(args.flow, args.fhigh, 1, args.m1, args.m2, args.s1, args.s2, ifo)**2
sigsqtt[ifo] = (1. / (2 * 3.14 * rho[ifo] * sigsqf**.5)**2)
- sigsqpp[ifo] = moment(args.flow, args.fhigh, 2, ifo) / (rho[ifo]**2 * sigsqf)
- sigsqtp[ifo] = moment(args.flow, args.fhigh, 1, ifo) / (2 * 3.14 * rho[ifo]**2 * sigsqf)
+ sigsqpp[ifo] = moment(args.flow, args.fhigh, 2, args.m1, args.m2, args.s1, args.s2, ifo) / (rho[ifo]**2 * sigsqf)
+ sigsqtp[ifo] = moment(args.flow, args.fhigh, 1, args.m1, args.m2, args.s1, args.s2, ifo) / (2 * 3.14 * rho[ifo]**2 * sigsqf)
sigsqdd[ifo] = 1. / rho[ifo]**2
transtt = {}
diff --git a/gstlal-inspiral/bin/gstlal_inspiral_create_dt_dphi_snr_ratio_pdfs_dag b/gstlal-inspiral/bin/gstlal_inspiral_create_dt_dphi_snr_ratio_pdfs_dag
index 3dcf5ef1c8..e8e83a2bf4 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_create_dt_dphi_snr_ratio_pdfs_dag
+++ b/gstlal-inspiral/bin/gstlal_inspiral_create_dt_dphi_snr_ratio_pdfs_dag
@@ -25,6 +25,10 @@ parser.add_argument('--psd-xml', help = 'XML containing HLV psd')
parser.add_argument('--H-snr', type = float, help = 'H characteristic SNR')
parser.add_argument('--L-snr', type = float, help = 'L characteristic SNR')
parser.add_argument('--V-snr', type = float, help = 'V characteristic SNR')
+parser.add_argument('--m1', type = float, default = 1.4, help = 'primary component mass')
+parser.add_argument('--m2', type = float, default = 1.4, help = 'secondary component mass')
+parser.add_argument('--s1', type = float, default = 0., help = 'primary (z) spin')
+parser.add_argument('--s2', type = float, default = 0., help = 'secondary (z) spin')
parser.add_argument('--flow', type = float, default = 10., help = 'Low frequency cut-off. Default 10 Hz')
parser.add_argument('--fhigh', type = float, default = 1024., help = 'High frequency cut-off. Default 1024 Hz')
args = parser.parse_args()
@@ -42,7 +46,7 @@ covJob = dagparts.DAGJob("gstlal_inspiral_compute_dtdphideff_cov_matrix", condor
margJob = dagparts.DAGJob("gstlal_inspiral_create_dt_dphi_snr_ratio_pdfs", condor_commands = {"request_memory":"7GB", "want_graceful_removal":"True", "kill_sig":"15", "accounting_group":"ligo.prod.o3.cbc.uber.gstlaloffline"})
addJob = dagparts.DAGJob("gstlal_inspiral_add_dt_dphi_snr_ratio_pdfs", condor_commands = {"request_memory":"4GB", "want_graceful_removal":"True", "kill_sig":"15", "accounting_group":"ligo.prod.o3.cbc.uber.gstlaloffline"})
-covnode = dagparts.DAGNode(covJob, dag, parent_nodes = [], opts = {"H-snr": args.H_snr, "L-snr": args.L_snr, "V-snr": args.V_snr, "flow": args.flow, "fhigh": args.fhigh}, input_files = {"psd-xml": args.psd_xml}, output_files = {"output":"covmat.h5"})
+covnode = dagparts.DAGNode(covJob, dag, parent_nodes = [], opts = {"H-snr": args.H_snr, "L-snr": args.L_snr, "V-snr": args.V_snr, "flow": args.flow, "fhigh": args.fhigh, "m1": args.m1, "m2": args.m2, "s1": args.s1, "s2": args.s2}, input_files = {"psd-xml": args.psd_xml}, output_files = {"output":"covmat.h5"})
num = 1000
margnodes = []
diff --git a/gstlal-inspiral/python/templates.py b/gstlal-inspiral/python/templates.py
index c522c4409a..6da52485e3 100644
--- a/gstlal-inspiral/python/templates.py
+++ b/gstlal-inspiral/python/templates.py
@@ -365,16 +365,9 @@ def time_slices(
return numpy.array(time_freq_boundaries,dtype=[('rate','int'),('begin','float'),('end','float')])
-def bandwidth(m1, m2, s1, s2, f_min, f_max, delta_f, psd):
- def fn(farr, h, n, psd):
- df = farr[1] - farr[0]
- return 4 * numpy.sum(numpy.abs(h)**2 * farr**n * df / psd(farr))
-
- def moment(farr, h, n, psd):
- return fn(farr, h, n, psd) / fn(farr, h, 0, psd)
-
- spin1 = [0., 0., s1]; spin2 = [0., 0., s2]
+def hplus_of_f(m1, m2, s1, s2, f_min, f_max, delta_f):
farr = numpy.linspace(0, f_max, f_max / delta_f + delta_f)
+ spin1 = [0., 0., s1]; spin2 = [0., 0., s2]
hp, hc = lalsim.SimInspiralFD(
m1 * lal.MSUN_SI, m2 * lal.MSUN_SI,
spin1[0], spin1[1], spin1[2],
@@ -394,5 +387,16 @@ def bandwidth(m1, m2, s1, s2, f_min, f_max, delta_f, psd):
)
assert hp.data.length > 0, "huh!? h+ has zero length!"
- h = hp.data.data[:len(farr)]
+ return hp.data.data[:len(farr)]
+
+def fn(farr, h, n, psd):
+ df = farr[1] - farr[0]
+ return 4 * numpy.sum(numpy.abs(h)**2 * farr**n * df / psd(farr))
+
+def moment(farr, h, n, psd):
+ return fn(farr, h, n, psd) / fn(farr, h, 0, psd)
+
+def bandwidth(m1, m2, s1, s2, f_min, f_max, delta_f, psd):
+ h = hplus_of_f(m1, m2, s1, s2, f_min, f_max, delta_f)
+ farr = numpy.linspace(0, f_max, f_max / delta_f + delta_f)
return (moment(farr, h, 2, psd) - moment(farr, h, 1, psd)**2)**.5
--
GitLab