Skip to content
Snippets Groups Projects
Commit 687a10b3 authored by Ryan Magee's avatar Ryan Magee Committed by Patrick Godwin
Browse files

gstlal_inspiral_compute_dtdphideff_cov_matrix: dont approximate pi, change...

gstlal_inspiral_compute_dtdphideff_cov_matrix: dont approximate pi, change spaces to tabs, fix snr not being stored
parent e30b5358
No related branches found
No related tags found
1 merge request!58dtdphi dag fixes
......@@ -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")
......
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