From e8477cb855dd6c17576eca1e643c4c0b991bcb18 Mon Sep 17 00:00:00 2001 From: Gareth S Cabourn Davies Date: Mon, 18 Sep 2023 16:05:17 +0100 Subject: [PATCH] Trigger timeseries ranking (#4491) * allow all rankings in pycbc_plot_trigger_timeseries * remove uinrelated chnage * Update bin/minifollowups/pycbc_plot_trigger_timeseries Co-authored-by: Tito Dal Canton * help message fix * Some more safety about unwanted events going away from the plot * more logging --------- Co-authored-by: Tito Dal Canton --- .../pycbc_plot_trigger_timeseries | 101 ++++++++++-------- 1 file changed, 55 insertions(+), 46 deletions(-) diff --git a/bin/minifollowups/pycbc_plot_trigger_timeseries b/bin/minifollowups/pycbc_plot_trigger_timeseries index 617cd491251..14ef0e64e51 100644 --- a/bin/minifollowups/pycbc_plot_trigger_timeseries +++ b/bin/minifollowups/pycbc_plot_trigger_timeseries @@ -19,6 +19,7 @@ import argparse, logging, pycbc.version, pycbc.results, sys from pycbc.types import MultiDetOptionAction from pycbc.events import ranking from pycbc.io import HFile +import h5py import matplotlib; matplotlib.use('Agg'); import pylab import numpy @@ -39,8 +40,9 @@ parser.add_argument('--times', nargs='+', type=float, parser.add_argument('--special-trigger-ids', nargs='+', type=int, action=MultiDetOptionAction, metavar="IFO:GPS_TIME", help="The set of special trigger ids to plot a star at") -parser.add_argument('--plot-type', choices=['snr', 'newsnr'], default='snr', - help="Which plot to make; an 'snr' or a newsnr' plot.") +parser.add_argument('--plot-type', + choices=ranking.sngls_ranking_function_dict, default='snr', + help="Which single-detector ranking statistic to plot.") parser.add_argument('--output-file') parser.add_argument('--log-y-axis', action='store_true') @@ -49,73 +51,80 @@ pycbc.init_logging(args.verbose) any_data = False -# organize the single detector triggers files by ifo in a dict fig = pylab.figure() + +min_rank = numpy.inf for ifo in args.single_trigger_files.keys(): + logging.info("Getting %s triggers", ifo) t = args.times[ifo] - if args.special_trigger_ids: - id_loud = args.special_trigger_ids[ifo] - data = HFile(args.single_trigger_files[ifo], 'r') # Identify trigger indices within window - select_func = lambda *inputs: ((inputs[0] < (t + args.window)) & - (inputs[0] > (t - args.window))) - times, snr, chisq, chisq_dof = data.select(select_func, ifo + '/end_time', - ifo + '/snr', ifo + '/chisq', - ifo + '/chisq_dof') - - # center times on the trigger/chosen time - times = times - t - - if args.plot_type == 'snr': - val = snr - - if args.special_trigger_ids: - top = data[ifo]['snr'][id_loud] - - if args.plot_type == 'newsnr': - rchisq = chisq / (2 * chisq_dof - 2) - if len(rchisq) > 0: - val = ranking.newsnr(snr, rchisq) - else: - val = numpy.array([]) - - # Get the newnsr of the loudest trigger so we can plot a star there - if args.special_trigger_ids: - rchisq = data[ifo]['chisq'][id_loud] / \ - (data[ifo]['chisq_dof'][id_loud] * 2 - 2) - top = ranking.newsnr(data[ifo]['snr'][id_loud], rchisq) - - if type(val) in [numpy.float32, numpy.float64] or len(val) > 0: - any_data = True - - pylab.scatter(times, val, color=pycbc.results.ifo_color(ifo), marker='x', + keep = abs(data[ifo + '/end_time'][:] - t) < args.window + if not any(keep): + # No triggers in this window, add to the legend and continue + # Make sure it isnt on the plot + pylab.scatter(-2 * args.window, 0, + color=pycbc.results.ifo_color(ifo), + marker='x', + label=ifo) + continue + + any_data = True + logging.info("Keeping %d triggers in the window", sum(keep)) + + + data_dict = {k: data[ifo][k][keep] for k in data[ifo] + if isinstance(data[ifo][k], h5py.Dataset) + and data[ifo][k].size == keep.size} + + logging.info("Getting %s", args.plot_type) + rank = ranking.get_sngls_ranking_from_trigs(data_dict, args.plot_type) + + if all(rank == 1): + # This is the default value to say "downranked beyond consideration" + # so skip these events + pylab.scatter(-2 * args.window, 0, + color=pycbc.results.ifo_color(ifo), + marker='x', + label=ifo) + continue + + pylab.scatter(data_dict['end_time'] - t, rank, + color=pycbc.results.ifo_color(ifo), marker='x', label=ifo) + # Minimum rank which hasn't been set to the default of 1 + min_rank = min(min_rank, rank[rank > 1].min()) + if args.special_trigger_ids: - any_data = True - pylab.scatter([0], [top], marker='*', s=50, color='yellow') + special_idx = numpy.flatnonzero(keep) == args.special_trigger_ids[ifo] + if not any(special_idx): + logging.info("IDX %d not in kept list", args.special_trigger_ids[ifo]) + continue + + pylab.scatter((data_dict['end_time'] - t)[special_idx], + rank[special_idx], marker='*', s=50, color='yellow') if args.log_y_axis and any_data: pylab.yscale('log') +if not numpy.isinf(min_rank): + pylab.ylim(ymin=min_rank) + pylab.xlabel('time (s)') pylab.ylabel(args.plot_type) -try: - if len(val) > 0: - pylab.ylim(ymin=val.min()) -except TypeError: - pylab.ylim(ymin=val) - pylab.xlim(xmin=-args.window, xmax=args.window) pylab.legend() pylab.grid() + +logging.info("Saving figure") pycbc.results.save_fig_with_metadata(fig, args.output_file, cmd = ' '.join(sys.argv), title = 'Single Detector Trigger Timeseries (%s)' % args.plot_type, caption = 'Time series showing the single detector triggers' ' centered around the time of the trigger of interest.', ) +logging.info("Done!")