Change cov matrix calculation to use actual waveform amplitude/parameters...

Change cov matrix calculation to use actual waveform amplitude/parameters instead of stationary phase approx

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 = {}

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 = []

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