gstlal_cherenkov_plot_rankingstat: initial commit for plotting program

gstlal-burst/bin/gstlal_cherenkov_plot_rankingstat

+#!/usr/bin/env python3
+#
+# Copyright (C) 2021  Soichiro Kuwahara
+#
+# 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.
+
+"""Cherenkov burst rankingstat plotting tool"""
+
+#
+# =============================================================================
+#
+#                                   Preamble
+#
+# =============================================================================
+#
+
+
+import math
+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
+})
+import numpy
+from optparse import OptionParser
+import sys
+
+from lal import rate
+from ligo.lw import utils as ligolw_utils
+from gstlal.cherenkov import rankingstat as cherenkov_rankingstat
+
+
+__author__ = "Soichrio Kuwahara <soichiro.kuwahara@ligo.org>"
+
+program_name = "gstlal_cherenkov_plot_rankingstat"
+
+
+#
+# =============================================================================
+#
+#                                 Command Line
+#
+# =============================================================================
+#
+
+
+def parse_command_line():
+	parser = OptionParser(
+		description = "GstLAL-based cherenkov burst ranking statistic plotting tool."
+	)
+	parser.add_option("--verbose", action = "store_true", help = "Be verbose.")
+	options, filenames = parser.parse_args()
+
+	# save for the process_params table
+	options.options_dict = dict(options.__dict__)
+
+	return options, filenames
+
+
+#
+# =============================================================================
+#
+#                                     Main
+#
+# =============================================================================
+#
+
+
+options, filenames = parse_command_line()
+
+
+#
+# load the ranking stat file
+#
+
+
+rankingstat = cherenkov_rankingstat.RankingStat.from_xml(ligolw_utils.load_filename(filenames[0], verbose = options.verbose, contenthandler = cherenkov_rankingstat.LIGOLWContentHandler), "rankingstat")
+
+
+#
+# plot denominator snr, chisq PDFs
+#
+
+
+for instrument, pdf in rankingstat.denominator.lnpdfs.items():
+	# FIXME:  this won't live here permanently, just doing it temporarily for now
+	pdf.normalize()
+
+	x,y = pdf.bins.centres()
+	lnpdf = pdf.at_centres()
+
+	fig = figure.Figure()
+	FigureCanvas(fig)
+	axes = fig.gca()
+
+	norm = matplotlib.colors.Normalize(vmin = -30., vmax = lnpdf.max())
+	mesh = axes.pcolormesh(x, y, lnpdf.T, norm = None, cmap = "afmhot", shading = "gouraud")
+	axes.contour(x, y, lnpdf.T, colors = "k", linestyles = "-", linewidths = .5, alpha = .3)
+	axes.loglog()
+	axes.grid(which = "both")
+	axes.set_xlim((4, 500))
+	axes.set_ylim((.0001, 5))
+	fig.colorbar(mesh, ax = axes)
+	axes.set_xlabel("$\mathrm{SNR}$")
+	axes.set_ylabel("$\chi^{2}/ \mathrm{SNR}^{2}$")
+	fig.savefig("test_pdf_%s_25.png" % instrument)
+
+
+	rate.filter_array(pdf.array, rate.gaussian_window(4., 4.))
+	pdf.normalize()
+	lnpdf = pdf.at_centres()
+
+	fig = figure.Figure()
+	FigureCanvas(fig)
+	axes = fig.gca()
+
+	norm = matplotlib.colors.Normalize(vmin = -30, vmax = lnpdf.max())
+	cmap = matplotlib.cm.afmhot
+	cmap.set_bad("k")
+	mesh = axes.pcolormesh(x, y, lnpdf.T, norm = norm, cmap = cmap, shading = "gouraud")
+	axes.contour(x, y, lnpdf.T, norm = norm, colors = "k", linestyles = "-", linewidths = .5, alpha = .3)
+	axes.loglog()
+	axes.grid(which = "both")
+	axes.set_xlim((4, 500))
+	axes.set_ylim((.0001, 5))
+	colorbar = fig.colorbar(mesh, ax = axes)
+	colorbar.set_clim(-30, lnpdf.max())
+	axes.set_xlabel("$\mathrm{SNR}$")
+	axes.set_ylabel("$\chi^{2}/ \mathrm{SNR}^{2}$")
+	fig.savefig("test_kde_pdf_%s_25.png" % instrument)