diff --git a/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix b/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix index 301984b77aecad4a431c4b9476c5b2f6b3dbdd63..10f837ba4e18807e197f671a45ddea91a960ec7b 100755 --- a/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix +++ b/gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix @@ -18,6 +18,7 @@ import argparse import itertools import numpy, scipy.interpolate +import numpy from lal import series from gstlal import templates from ligo.lw import utils as ligolw_utils @@ -40,10 +41,8 @@ args = parser.parse_args() if args.flow >= args.fhigh: raise ValueError("flow cannot be greater than fhigh") - refpsd = args.psd_xml rho = {"L1": args.L_snr, "H1": args.H_snr, "V1": args.V_snr, "K1": args.K_snr} - instruments = ("H1", "L1", "V1", "K1") snr = series.read_psd_xmldoc(ligolw_utils.load_filename(refpsd, verbose = True, contenthandler = series.PSDContentHandler)) psddict = {} @@ -63,9 +62,9 @@ sigsqdd = {} for ifo in rho: 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) + sigsqtt[ifo] = (1. / (2 * numpy.pi * rho[ifo] * sigsqf**.5)**2) 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) + sigsqtp[ifo] = moment(args.flow, args.fhigh, 1, args.m1, args.m2, args.s1, args.s2, ifo) / (2 * numpy.pi * rho[ifo]**2 * sigsqf) sigsqdd[ifo] = 1. / rho[ifo]**2 transtt = {} @@ -93,10 +92,10 @@ for pair in pairs: transdd[frozenset(pair)] = 1. / numpy.sqrt(sigsqdd[a] + sigsqdd[b]) f = h5py.File(args.output, "w") -f.create_dataset("psd", data=args.psd_xml) +f.create_dataset("psd", shape = (1,), data=refpsd) h5_snr = f.create_group("SNR") for ifo, snr in rho.items(): - h5_snr.create_dataset(ifo, data=snr) + h5_snr.create_dataset(ifo, data=numpy.array([snr])) h5_transtt = f.create_group("transtt") h5_transtp = f.create_group("transtp") h5_transpt = f.create_group("transpt")