From 24f8ea05e57993e618ed417b67bbbc3d426d7664 Mon Sep 17 00:00:00 2001
From: Kipp Cannon <kipp.cannon@ligo.org>
Date: Thu, 24 Aug 2017 17:53:00 +0900
Subject: [PATCH] plotpsd: add cumulative SNRs plot
---
gstlal/python/plotpsd.py | 67 ++++++++++++++++++++++++++++++++++++++++
1 file changed, 67 insertions(+)
diff --git a/gstlal/python/plotpsd.py b/gstlal/python/plotpsd.py
index 96a93bedc0..9dba8f8591 100644
--- a/gstlal/python/plotpsd.py
+++ b/gstlal/python/plotpsd.py
@@ -61,6 +61,52 @@ def summarize_coinc_xmldoc(coinc_xmldoc):
return sngl_inspirals, mass1, mass2, end_time, on_instruments
+def axes_plot_cummulative_snr(axes, psds, coinc_xmldoc):
+ sngl_inspirals, mass1, mass2, end_time, on_instruments = summarize_coinc_xmldoc(coinc_xmldoc)
+
+ axes.grid(which = "both", linestyle = "-", linewidth = 0.2)
+ axes.minorticks_on()
+
+ for instrument, sngl_inspiral in sngl_inspirals.items():
+ logging.info("found %s event with SNR %g" % (instrument, sngl_inspirals[instrument].snr))
+
+ if instrument not in psds:
+ logging.info("no PSD for %s" % instrument)
+ continue
+ psd = psds[instrument]
+ if psd is None:
+ logging.info("no PSD for %s" % instrument)
+ continue
+ psd_data = psd.data.data
+ f = psd.f0 + numpy.arange(len(psd_data)) * psd.deltaF
+ logging.info("found PSD for %s spanning [%g Hz, %g Hz]" % (instrument, f[0], f[-1]))
+
+ # FIXME: horizon distance stopped at 0.9 max frequency due
+ # to low pass filter messing up the end of the PSD. if we
+ # could figure out the frequency bounds and delta F we
+ # could move this out of the loop for some speed
+ horizon_distance = reference_psd.HorizonDistance(10., 0.9 * f[-1], psd.deltaF, mass1, mass2)
+
+ # generate inspiral spectrum and clip PSD to its frequency
+ # range
+ inspiral_spectrum_x, inspiral_spectrum_y = horizon_distance(psd, sngl_inspiral.snr)[1]
+ lo = int(round((inspiral_spectrum_x[0] - psd.f0) / psd.deltaF))
+ hi = int(round((inspiral_spectrum_x[-1] - psd.f0) / psd.deltaF)) + 1
+ f = f[lo:hi]
+ psd_data = psd_data[lo:hi]
+
+ # plot
+ snr2 = (inspiral_spectrum_y / psd_data).cumsum() * psd.deltaF
+ axes.semilogx(f, snr2**.5, color = plotutil.colour_from_instruments([instrument]), alpha = 0.8, linestyle = "-", label = "%s SNR = %.3g" % (instrument, sngl_inspiral.snr))
+
+ axes.set_ylim([0., axes.get_ylim()[1]])
+
+ axes.set_title(r"Cumulative SNRs for $%.3g\,\mathrm{M}_{\odot}$--$%.3g\,\mathrm{M}_{\odot}$ Merger Candidate at %.2f GPS" % (mass1, mass2, float(end_time)))
+ axes.set_xlabel(r"Frequency (Hz)")
+ axes.set_ylabel(r"Cumulative SNR")
+ axes.legend(loc = "upper left")
+
+
def axes_plot_psds(axes, psds, coinc_xmldoc = None):
"""!
Places a PSD plot into a matplotlib Axes object.
@@ -158,3 +204,24 @@ def plot_psds(psds, coinc_xmldoc = None, plot_width = 640):
axes_plot_psds(fig.gca(), psds, coinc_xmldoc = coinc_xmldoc)
fig.tight_layout(pad = .8)
return fig
+
+
+def plot_cumulative_snrs(psds, coinc_xmldoc, plot_width = 640):
+ """!
+ Produces a matplotlib figure of cumulative SNRs.
+
+ @param psds A dictionary of PSDs as REAL8FrequencySeries keyed by
+ instrument
+
+ @param coinc_xmldoc An XML document containing a single event with all
+ of the metadata as would be uploaded to gracedb.
+
+ @param plot_width How wide to make the figure object in pixels
+ (ignored if axes is provided).
+ """
+ fig = figure.Figure()
+ FigureCanvas(fig)
+ fig.set_size_inches(plot_width / float(fig.get_dpi()), int(round(plot_width / plotutil.golden_ratio)) / float(fig.get_dpi()))
+ axes_plot_cummulative_snr(fig.gca(), psds, coinc_xmldoc)
+ fig.tight_layout(pad = .8)
+ return fig
--
GitLab