From 7a7b68f99a0df82e6dacb92310a8a0a963723f97 Mon Sep 17 00:00:00 2001
From: Leo Tsukada <leo.tsukada@ligo.org>
Date: Thu, 1 Jun 2023 18:49:17 -0700
Subject: [PATCH] gstlal-inspiral/bin/gstlal_ll_inspiral_event_plotter : add
 DTDPHI_PLOTS mode to upload dtdphi plot

---
 .../bin/gstlal_ll_inspiral_event_plotter      | 47 +++++++++++++++++--
 1 file changed, 44 insertions(+), 3 deletions(-)

diff --git a/gstlal-inspiral/bin/gstlal_ll_inspiral_event_plotter b/gstlal-inspiral/bin/gstlal_ll_inspiral_event_plotter
index adacc1a781..32aef85a7f 100755
--- a/gstlal-inspiral/bin/gstlal_ll_inspiral_event_plotter
+++ b/gstlal-inspiral/bin/gstlal_ll_inspiral_event_plotter
@@ -51,6 +51,7 @@ from gstlal import events
 from gstlal import far
 from gstlal import inspiral
 from gstlal import lvalert_helper
+from gstlal.plots import dtdphi as plotdtdphi
 from gstlal.plots import far as plotfar
 from gstlal.plots import psd as plotpsd
 
@@ -112,6 +113,7 @@ def parse_command_line():
 		'2. RANKING_PLOTS: upload background plots from the ranking statistic data. '
 		'3. SNR_PLOTS: upload SNR time series plots. '
 		'4. PSD_PLOTS: upload PSD plots. '
+		'5. DTDPHI_PLOTS: upload DTDPHI plots. '
 		' Can be given multiple times. If not specified, all plots are made.'))
 
 	options, args = parser.parse_args()
@@ -130,6 +132,7 @@ class Plots(Enum):
 	RANKING_PLOTS = 2
 	SNR_PLOTS = 3
 	PSD_PLOTS = 4
+	DTDPHI_PLOTS = 5
 
 class EventPlotter(events.EventProcessor):
 	"""
@@ -210,7 +213,7 @@ class EventPlotter(events.EventProcessor):
 			logging.info('found new event at {} from bin {}'.format(time, bank_bin))
 			self.events[event_key] = self.new_event(time, bank_bin)
 
-		### ranking stat 
+		### ranking stat
 		if message.topic() == self.ranking_stat_topic:
 			self.events[event_key]['ranking_data_path'] = payload['ranking_data_path']
 			# we'll just take the xmldoc from the preferred event, which will be identical
@@ -254,11 +257,14 @@ class EventPlotter(events.EventProcessor):
 					self.upload_ranking_plots(event)
 					uploaded['RANKING_PLOTS'] = True
 				if 'PSD_PLOTS' in uploaded and not uploaded['PSD_PLOTS'] and event['psd']:
-					self.upload_psd_plots(event) 
+					self.upload_psd_plots(event)
 					uploaded['PSD_PLOTS'] = True
 				if 'SNR_PLOTS' in uploaded and not uploaded['SNR_PLOTS']:
 					self.upload_snr_plots(event)
 					uploaded['SNR_PLOTS'] = True
+				if 'DTDPHI_PLOTS' in uploaded and not uploaded['DTDPHI_PLOTS'] and event['ranking_data_path']:
+					self.upload_dtdphi_plots(event)
+					uploaded['DTDPHI_PLOTS'] = True
 
 		# clean out events once all plots are uploaded
 		# and clean out old events
@@ -483,7 +489,7 @@ class EventPlotter(events.EventProcessor):
 			ax.set_xlabel(r'$\mathrm{{Time\,from\,{}}}$'.format(peaktime))
 			ax.legend(loc='best')
 			ax.grid()
-			
+
 		fig.tight_layout()
 		filename = '{}_snrtimeseries.{}'.format(event['gid'], self.format)
 
@@ -498,6 +504,41 @@ class EventPlotter(events.EventProcessor):
 
 		logging.info('finished processing SNR time series plot for {}'.format(event['gid']))
 
+
+	def upload_dtdphi_plots(self, event):
+		sngl_inspiral_table = lsctables.SnglInspiralTable.get_table(event['coinc'])
+		offset_vectors = lsctables.TimeSlideTable.get_table(event['coinc']).as_dict()
+		assert len(offset_vectors.keys()) == 1, "the time slide table has multiple time-slide entries"
+
+		offset_vector = offset_vectors[list(offset_vectors.keys())[0]]
+		rankingstat, _ = far.parse_likelihood_control_doc(ligolw_utils.load_filename(event['ranking_data_path'], contenthandler = far.RankingStat.LIGOLWContentHandler))
+
+		event_kwargs = rankingstat.kwargs_from_triggers(sngl_inspiral_table, offset_vector)
+		ifos = sorted(event_kwargs["snrs"].keys())
+		ifo1 = ifos.pop(0)
+		snrs = event_kwargs["snrs"]
+		# remove segments for ifos not having horizon distance after
+		# compression, which would cause an error inside
+		# local_mean_horizon_distance()
+		segs_nonzero = {ifo:seg for ifo, seg in event_kwargs["segments"].items() if len(rankingstat.numerator.horizon_history[ifo]) > 0}
+		horizons = rankingstat.numerator.local_mean_horizon_distance(segs_nonzero, template_id = event_kwargs["template_id"])
+
+		### generate and upload plots
+		for ifo2 in ifos:
+			ifo_pair = ifo1[0] + ifo2[0]
+			dt_ref = event_kwargs["dt"][ifo2] - event_kwargs["dt"][ifo1]
+			dphi_ref = event_kwargs["phase"][ifo2] - event_kwargs["phase"][ifo1]
+
+			fig = plotdtdphi.plots_dtdphi(ifo1, ifo2, snrs, horizons, sngl={"dt":dt_ref, "dphi":dphi_ref})
+			filename = '{}_dtdphi_{}.{}'.format(event['gid'], ifo_pair, self.format)
+			if not self.no_upload:
+				lvalert_helper.upload_fig(fig, self.client, event['gid'], filename = filename, log_message = '{} dtdphi 2D pdf plot'.format(ifo_pair), tagname = ('dtdphi', 'background'))
+			if self.output_path is not None:
+				filename = os.path.join(self.output_path, filename)
+				logging.info('writing {} ...'.format(filename))
+				fig.savefig(filename)
+			logging.info('finished processing {} dtdphi pdf plot for {}'.format(ifo_pair, event['gid']))
+
 	def finish(self):
 		"""
 		upload remaining files/plots before shutting down
-- 
GitLab