gstlal_inspiral_compute_dtdphideff_cov_matrix : added KAGRA

gstlal-inspiral/bin/gstlal_inspiral_compute_dtdphideff_cov_matrix

@@ -19,7 +19,6 @@ import argparse
 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
@@ -46,16 +45,12 @@ if args.flow >= args.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")
 psd = series.read_psd_xmldoc(ligolw_utils.load_filename(refpsd, verbose = True, contenthandler = series.PSDContentHandler))
 psddict = {}
 for ifo in rho:
-	ifo_psd = ifo
-	# FIXME use Virgo PSD for KAGRA temporarily. This has to be fixed once KAGRA PSD is measured.
-	if ifo == "K1":
-		ifo = "V1"
 	f = numpy.arange(len(psd[ifo].data.data)) * psd[ifo].deltaF
-	psddict[ifo_psd] = scipy.interpolate.interp1d(f, psd[ifo].data.data)
+	psddict[ifo] = scipy.interpolate.interp1d(f, psd[ifo].data.data)
 
 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)
@@ -83,7 +78,7 @@ transdd = {}
 # transdd = {frozenset(("H1", "L1")): 1. / numpy.sqrt(sigsqdd["H1"] + sigsqdd["L1"])  , frozenset(("H1", "V1")): 1. / numpy.sqrt(sigsqdd["H1"] + sigsqdd["V1"]), frozenset(("L1", "V1")): 1. / numpy.sqrt(sigsqdd["L1"] + sigsqdd["V1"])}
 # transdd = {frozenset(instrument_pair) : 1. / numpy.sqrt(sigsqdd[instrument_pair[0]] + sigsqdd[instrument_pair[1]]) for instrument_pair in TimePhaseSNR.instument_pairs(("H1", "L1", "V1", "K1"))}
 
-pairs = [tuple(sorted(pair)) for pair in itertools.combinations(("H1", "L1", "V1", "K1"), 2)]
+pairs = [tuple(sorted(pair)) for pair in itertools.combinations(instruments, 2)]
 for pair in pairs:
 	a,b = pair
 	m11 = sigsqtt[a] + sigsqtt[b]
@@ -108,13 +103,17 @@ for group, mat in zip((h5_transtt, h5_transtp, h5_transpt, h5_transpp, h5_transd
 	for k,v in mat.items():
 		group.create_dataset(",".join(sorted(k)), data = float(v))
 
-combos = TimePhaseSNR.instrument_combos(("H1","L1","V1", "K1"), min_instruments=1)
+combos = set()
+for i in range(1, len(instruments) + 1):
+	for choice in itertools.combinations(instruments, i):
+		combos.add(tuple(sorted(choice)))
+combos = tuple(sorted(list(combos)))
 norm = dict((frozenset(k), 0.) for k in combos)
 norm = {}
 h5_norm = f.create_group("norm")
-for k in combos:
-	h5_norm.create_dataset(",".join(sorted(k)), data = 1.0)
-	norm[k] = 1.0
+for combo in combos:
+	h5_norm.create_dataset(",".join(sorted(combo)), data = 1.0)
+	norm[combo] = 1.0
 f.close()
 
 print "transtt =", transtt