From ae776f2017d4605b220319c1f3ef8cf216746d79 Mon Sep 17 00:00:00 2001
From: Sarah Caudill <sarah.caudill@ligo.org>
Date: Tue, 15 May 2018 06:14:50 -0700
Subject: [PATCH] gstlal_inspiral_plotsummary: adding support for Virgo SNR

---
 .../bin/gstlal_inspiral_plotsummary           | 79 ++++++++++++-------
 1 file changed, 50 insertions(+), 29 deletions(-)

diff --git a/gstlal-inspiral/bin/gstlal_inspiral_plotsummary b/gstlal-inspiral/bin/gstlal_inspiral_plotsummary
index c4b3855ba3..0d9e1560c6 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_plotsummary
+++ b/gstlal-inspiral/bin/gstlal_inspiral_plotsummary
@@ -604,17 +604,24 @@ SELECT distinct_ifos.ifos, count(*) FROM coinc_inspiral JOIN distinct_ifos ON (d
 			("L1 m<sub>1</sub>", "number"),
 			("L1 m<sub>2</sub>", "number"),
 			("L1 s<sub>1z</sub>", "number"),
-			("L1 s<sub>2z</sub>", "number")
+			("L1 s<sub>2z</sub>", "number"),
+			("V1 &rho;", "string"),
+			("V1 &chi;<sup>2</sup>", "number"),
+			("V1 m<sub>1</sub>", "number"),
+			("V1 m<sub>2</sub>", "number"),
+			("V1 s<sub>1z</sub>", "number"),
+			("V1 s<sub>2z</sub>", "number")
 		]
 		for rank, values in enumerate(candidates, 1):
 			row = [rank] + [float(v) for v in values[:7]] + list(values[7:9])
 			# values[9] is a string that is e.g., H1:4.8993754:1.0139208:2.061641:1.145543 L1:8.2582664:1.1890973:2.061641:1.145543
-			ifodict = {"H1": [-1,-1,-1,-1,-1,-1], "L1": [-1,-1,-1,-1,-1,-1]}
+			ifodict = {"H1": [-1,-1,-1,-1,-1,-1], "L1": [-1,-1,-1,-1,-1,-1], "V1": [-1,-1,-1,-1,-1,-1]}
 			for ifo_row in values[9].split():
 				ifodata = ifo_row.split(":")
 				ifodict[ifodata[0]] = [float(v) for v in ifodata[1:]]
 			row.extend(ifodict["H1"])
 			row.extend(ifodict["L1"])
+			row.extend(ifodict["V1"])
 			data.append(row)
 		to_google_json(fname, description, data)
 
@@ -765,35 +772,41 @@ FROM
 					return sim.get_chirp_eff_dist(list(instruments)[0])
 
 			def found_decisive_charsnr(sim, oninstruments, partinstruments):
-				if len(oninstruments) == 2:
-					if len(partinstruments) == 2:
-						return sorted([sim.alpha1, sim.alpha2])[0]
-					if len(partinstruments) == 1:
-						ifo, = partinstruments
-						if ifo == 'H1':
-							return sim.alpha1
-						if ifo == 'L1':
-							return sim.alpha2
-				elif len(oninstruments) == 1:
+				if len(oninstruments) == 1:
 					ifo, = oninstruments
 					if ifo == 'H1':
-						return sim.alpha1
+						return sim.alpha4
 					if ifo == 'L1':
-						return sim.alpha2
+						return sim.alpha5
+					if ifo == 'V1':
+						return sim.alpha6
+				if len(oninstruments) > 3:
+					raise ValueError("More than 3 instruments not supported for injection snr calculation at this time.")
 				else:
-					raise ValueError("More than 2 instruments not supported for injection snr calculation at this time.")
+					if len(partinstruments) == 1:
+						ifo, = partinstruments
+						if ifo == 'H1':
+							return sim.alpha4
+						if ifo == 'L1':
+							return sim.alpha5
+						if ifo == 'V1':
+							return sim.alpha6
+					else:
+						return sorted([sim.alpha4, sim.alpha5, sim.alpha6])[1]
 
 			def missed_decisive_charsnr(sim, oninstruments):
-				if len(oninstruments) == 2:
-					return sorted([sim.alpha1, sim.alpha2])[0]
-				elif len(oninstruments) == 1:
+				if len(oninstruments) == 1:
 					ifo, = oninstruments
 					if ifo == 'H1':
-						return sim.alpha1
+						return sim.alpha4
 					if ifo == 'L1':
-						return sim.alpha2
+						return sim.alpha5
+					if ifo == 'V1':
+						return sim.alpha6
+				if len(oninstruments) > 3:
+					raise ValueError("More than 3 instruments not supported for injection snr calculation at this time.")
 				else:
-					raise ValueError("More than 2 instruments not supported for injection snr calculation at this time.")
+					return sorted([sim.alpha4, sim.alpha5, sim.alpha6])[1]
 
 			for cnt, (title, x_label, x_func, y_label, y_func, filename_fragment) in enumerate((
 				(r"$\textrm{Distance vs.\ Chirp Mass (With %s Operating)}$" % ", ".join(sorted(self.on_instruments)), r"$M_{\mathrm{chirp}}$ ($\mathrm{M}_{\odot}$)", lambda sim: sim.mchirp, r"$D$ ($\mathrm{Mpc}$)", lambda sim, instruments: sim.distance, "d_vs_mchirp"),
@@ -816,7 +829,7 @@ FROM
 					if cnt == 0:
 						self.missed_found_plots.injection_summary_data.append(["Missed", "".join(sorted(self.on_instruments)), "---", len(missed)])
 						for rank, sim in enumerate(missed):
-							self.missed_found_plots.missed_summary_data.append(["".join(sorted(self.on_instruments)), sim.waveform, float(sim.time_at_instrument("H1", {"H1": 0.0})), float(sim.time_at_instrument("L1", {"L1": 0.0})), sim.mass1, sim.mass2, sim.spin1x, sim.spin1y, sim.spin1z, sim.spin2x, sim.spin2y, sim.spin2z, sim.distance, decisive_chirp_distance(sim, self.on_instruments), sim.inclination, sim.alpha1, sim.alpha2, missed_decisive_charsnr(sim, self.on_instruments)])
+							self.missed_found_plots.missed_summary_data.append(["".join(sorted(self.on_instruments)), sim.waveform, float(sim.time_at_instrument("H1", {"H1": 0.0})), float(sim.time_at_instrument("L1", {"L1": 0.0})), float(sim.time_at_instrument("V1", {"V1": 0.0})), sim.mass1, sim.mass2, sim.spin1x, sim.spin1y, sim.spin1z, sim.spin2x, sim.spin2y, sim.spin2z, sim.distance, decisive_chirp_distance(sim, self.on_instruments), sim.inclination, sim.alpha4, sim.alpha5, sim.alpha6, missed_decisive_charsnr(sim, self.on_instruments)])
 					legend.append("Missed")
 					axes.semilogy([x_func(sim) for sim in missed], [y_func(sim, self.on_instruments) for sim in missed], "k.")
 
@@ -874,7 +887,7 @@ FROM
 				yield fig, "%s_%s" % (filename_fragment, "".join(sorted(on_instruments))), is_open_box
 
 		to_google_json(self.base + "_injection_summary.json", [("Found||Missed", "string"), ("On Instruments", "string"), ("Participating Instruments", "string"), ("count", "number")], self.injection_summary_data)
-		to_google_json(self.base + "_missed_summary.json", [("On Instruments", "string"), ("Waveform", "string"), ("H1 <i>t</i>", "number"), ("L1 <i>t</i>", "number"), ("<i>m</i><sub>1</sub>", "number"), ("<i>m</i><sub>2</sub>", "number"), ("<i>s</i><sub>1x</sub>", "number"), ("<i>s</i><sub>1y</sub>", "number"), ("<i>s</i><sub>1z</sub>", "number"), ("<i>s</i><sub>2x</sub>", "number"), ("<i>s</i><sub>2y</sub>", "number"), ("<i>s</i><sub>2z</sub>", "number"), ("D (Mpc)", "number"), ("Decisive D<sub>chirp,eff</sub> (Mpc)" , "number"), ("Inclination", "number"), ("H1 &rho;", "number"), ("L1 &rho;", "number"), ("Decisive &rho;", "number")], self.missed_summary_data)
+		to_google_json(self.base + "_missed_summary.json", [("On Instruments", "string"), ("Waveform", "string"), ("H1 <i>t</i>", "number"), ("L1 <i>t</i>", "number"), ("V1 <i>t</i>", "number"), ("<i>m</i><sub>1</sub>", "number"), ("<i>m</i><sub>2</sub>", "number"), ("<i>s</i><sub>1x</sub>", "number"), ("<i>s</i><sub>1y</sub>", "number"), ("<i>s</i><sub>1z</sub>", "number"), ("<i>s</i><sub>2x</sub>", "number"), ("<i>s</i><sub>2y</sub>", "number"), ("<i>s</i><sub>2z</sub>", "number"), ("D (Mpc)", "number"), ("Decisive D<sub>chirp,eff</sub> (Mpc)" , "number"), ("Inclination", "number"), ("H1 &rho;", "number"), ("L1 &rho;", "number"), ("V1 &rho;", "number"), ("Decisive &rho;", "number")], self.missed_summary_data)
 
 
 #
@@ -1004,19 +1017,27 @@ class ParameterAccuracyPlots(object):
 			if instrument == 'H1':
 				fig, axes = create_plot(r"Inj. H1 SNR", r"Rec. H1 SNR")
 				axes.set_title(r"SNR Recovery in %s (%s Injections)" % (instrument, waveform if "_" not in waveform else waveform.replace("_", "\_")))
-				if max([sim.alpha1 for sim, sngl, far in pairs]) > 20:
-					axes.loglog([sim.alpha1 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
+				if max([sim.alpha4 for sim, sngl, far in pairs]) > 20:
+					axes.loglog([sim.alpha4 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
 				else:
-					axes.plot([sim.alpha1 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
+					axes.plot([sim.alpha4 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
 				yield fig, "snr_rec_scatter_%s_%s" % (waveform, instrument), False
 
 			if instrument == 'L1':
 				fig, axes = create_plot(r"Inj. L1 SNR", r"Rec. L1 SNR")
 				axes.set_title(r"SNR Recovery in %s (%s Injections)" % (instrument, waveform if "_" not in waveform else waveform.replace("_", "\_")))
-				if max([sim.alpha2 for sim, sngl, far in pairs]) > 20:
-					axes.loglog([sim.alpha2 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
+				if max([sim.alpha5 for sim, sngl, far in pairs]) > 20:
+					axes.loglog([sim.alpha5 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
+				else:
+					axes.plot([sim.alpha5 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
+				yield fig, "snr_rec_scatter_%s_%s" % (waveform, instrument), False
+			if instrument == 'V1':
+				fig, axes = create_plot(r"Inj. V1 SNR", r"Rec. V1 SNR")
+				axes.set_title(r"SNR Recovery in %s (%s Injections)" % (instrument, waveform if "_" not in waveform else waveform.replace("_", "\_")))
+				if max([sim.alpha6 for sim, sngl, far in pairs]) > 20:
+					axes.loglog([sim.alpha6 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
 				else:
-					axes.plot([sim.alpha2 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
+					axes.plot([sim.alpha6 for sim, sngl, far in pairs], [sngl.snr for sim, sngl, far in pairs], "kx")
 				yield fig, "snr_rec_scatter_%s_%s" % (waveform, instrument), False
 
 #
-- 
GitLab