gstlal-inspiral/python/plots/dtdphi.py, Makefile.am : add a routine to create...

gstlal-inspiral/python/plots/dtdphi.py, Makefile.am : add a routine to create dtdphi plot for a given coinc

gstlal-inspiral/python/plots/Makefile.am

@@ -3,6 +3,7 @@ plotsdir = $(pkgpythondir)/plots
 
 plots_PYTHON = \
 	bank.py \
+	dtdphi.py \
 	far.py \
 	horizon.py \
 	lloid.py \

gstlal-inspiral/python/plots/dtdphi.py

+# Copyright (C) 2023  Leo Tsukada
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+## @file
+
+## @package plotdtdphi
+
+#
+# =============================================================================
+#
+#                                   Preamble
+#
+# =============================================================================
+#
+
+import matplotlib
+from matplotlib import figure
+from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
+matplotlib.rcParams.update({
+	"font.size": 10.0,
+	"axes.titlesize": 10.0,
+	"axes.labelsize": 10.0,
+	"xtick.labelsize": 8.0,
+	"ytick.labelsize": 8.0,
+	"legend.fontsize": 8.0,
+	"figure.dpi": 300,
+	"savefig.dpi": 300,
+	"text.usetex": True
+})
+from matplotlib import pyplot
+import numpy
+import lal
+from lalburst.snglcoinc import light_travel_time
+from gstlal.stats import inspiral_extrinsics
+
+TPDPDF = inspiral_extrinsics.InspiralExtrinsics()
+R = {"H1":lal.CachedDetectors[lal.LHO_4K_DETECTOR].response,"L1":lal.CachedDetectors[lal.LLO_4K_DETECTOR].response,"V1":lal.CachedDetectors[lal.VIRGO_DETECTOR].response,"K1":lal.CachedDetectors[lal.KAGRA_DETECTOR].response}
+X = {"H1":lal.CachedDetectors[lal.LHO_4K_DETECTOR].location,"L1":lal.CachedDetectors[lal.LLO_4K_DETECTOR].location,"V1":lal.CachedDetectors[lal.VIRGO_DETECTOR].location,"K1":lal.CachedDetectors[lal.KAGRA_DETECTOR].location}
+
+def dtdphi2prob(dt, dphi, ifo1, ifo2, refsnr, refhorizon):
+	"""function that takes dt,dphi samples and return the probability (density) values derived from the new pdf
+
+	:dt: the time delay between detector pair
+	:dphi: the difference of the coalescence phase between detector pair
+	:returns: dtdpipdf(dt, dphi) * snr^4 (assuming snr = {"H1": 5., "V1": 2.25} and horizon = {"H1": 110., "V1": 45.})
+
+	"""
+
+	t = {ifo1: 0, ifo2: dt}
+	phi = {ifo1: 0, ifo2: dphi}
+	snr = {ifo1: refsnr[ifo1], ifo2: refsnr[ifo2]}
+	horizon = {ifo1: refhorizon[ifo1], ifo2: refhorizon[ifo2]}
+	# signature is (time, phase, snr, horizon)
+	p = TPDPDF.time_phase_snr(t, phi, snr, horizon) * (sum(s**2 for s in snr.values())**.5)**4
+	return float(p)
+
+
+def init_plot(figsize):
+	fig = figure.Figure()
+	FigureCanvas(fig)
+	fig.set_size_inches(figsize)
+	axes = fig.gca()
+
+	return fig, axes
+
+
+
+def plots_dtdphi(ifo1, ifo2, snrs, horizons, sngl=None):
+	len_arr = 100
+	dt_max = light_travel_time(ifo1, ifo2)
+	dtvec = numpy.linspace(-dt_max, dt_max, len_arr)
+	dphivec = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, len_arr)
+	DTProb = numpy.zeros(len_arr)
+	DPProb = numpy.zeros(len_arr)
+	Prob = numpy.zeros((len_arr,len_arr))
+
+	for j, dt in enumerate(dtvec):
+		for i, dphi in enumerate(dphivec):
+			p = dtdphi2prob(dt, dphi, ifo1, ifo2, snrs, horizons)
+			Prob[j,i] = p
+			DPProb[i] += p
+			DTProb[j] += p
+
+	fig = pyplot.figure(figsize=(7.5,7.5))
+	ax1 = fig.add_subplot(221)
+	ax1.plot(dphivec, DPProb / numpy.sum(DPProb) / (dphivec[1] - dphivec[0]), label ="Direct Calculation")
+	ax1.set_ylabel(r"""$P(\Delta\phi | s, \{D_{%s}, D_{%s}\}, \{\rho_{%s}, \rho_{%s}\})$""" % (ifo1, ifo2, ifo1, ifo2))
+	ax1.legend(bbox_to_anchor=(1, 1), loc='upper left', ncol=1)
+
+	ax2 = fig.add_subplot(223)
+	ax2.pcolor(dphivec, dtvec, Prob)
+	ax2.set_xlabel(r"""$\phi_{%s} - \phi_{%s}$""" % (ifo2, ifo1))
+	ax2.set_ylabel(r"""$t_{%s} - t_{%s}$""" % (ifo2, ifo1))
+
+	ax3 = fig.add_subplot(224)
+	ax3.plot(DTProb / numpy.sum(DTProb) / (dtvec[1] - dtvec[0]), dtvec)
+	ax3.set_xlabel(r"""$P(\Delta t | s, \{D_{%s}, D_{%s}\}, \{\rho_{%s}, \rho_{%s}\})$""" % (ifo1, ifo2, ifo1, ifo2))
+
+	if sngl is not None:
+		dt_ref = sngl["dt"]
+		dphi_ref = sngl["dphi"]
+		ax1.axvline(dphi_ref, ls="--", color = "r", lw=4)
+		ax2.plot(dphi_ref, dt_ref, "ko", mfc = "None", mec = "r", ms = 14, mew=4)
+		ax3.axhline(dt_ref, ls="--", color = "r", lw=4)
+	fig.tight_layout(pad = .8)
+
+	return fig