Skip to content
Snippets Groups Projects
Commit b5f55795 authored by Ryan Magee's avatar Ryan Magee
Browse files

Revert "gstlal_inspiral_compute_dtdphideff_cov_matrix: dont approximate pi,...

Revert "gstlal_inspiral_compute_dtdphideff_cov_matrix: dont approximate pi, change spaces to tabs, fix snr not being stored"

This reverts commit 687a10b3.
parent 9ba7faab
No related branches found
No related tags found
No related merge requests found
......@@ -18,7 +18,6 @@
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
......@@ -41,8 +40,10 @@ 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 = {}
......@@ -62,9 +63,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 * numpy.pi * rho[ifo] * sigsqf**.5)**2)
sigsqtt[ifo] = (1. / (2 * 3.14 * 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 * numpy.pi * 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 = {}
......@@ -92,10 +93,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", shape = (1,), data=refpsd)
f.create_dataset("psd", data=args.psd_xml)
h5_snr = f.create_group("SNR")
for ifo, snr in rho.items():
h5_snr.create_dataset(ifo, data=numpy.array([snr]))
h5_snr.create_dataset(ifo, data=snr)
h5_transtt = f.create_group("transtt")
h5_transtp = f.create_group("transtp")
h5_transpt = f.create_group("transpt")
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment