Add details of closest triggers to the injection minifollowup (#4503)

* Add details of closest triggers (not just event) to the injection minifollowup

* Dont know how this was deleted

* Wrap titles with quotes

Co-authored-by: Ian Harry <[email protected]>

* Use existing functions to reduce the number of trigger times being stored

* Let the max window be set as a command line argument, deal with case where closest trigger is far away from the injection

* remove import I thought I would use but didnt

* Use the loudest SNR event in a tiny window instead of simply closest time

* missed re-setting the ifo in this loop

* Apply suggestions from code review

* TD comment

---------

Co-authored-by: Ian Harry <[email protected]>

bin/minifollowups/pycbc_injection_minifollowup

@@ -25,7 +25,7 @@ from pycbc.detector import Detector
 import pycbc.workflow.minifollowups as mini
 from pycbc.workflow.core import resolve_url_to_file
 import pycbc.version, pycbc.pnutils
-from pycbc.io.hdf import SingleDetTriggers
+from pycbc.io.hdf import SingleDetTriggers, HFile
 
 parser = argparse.ArgumentParser(description=__doc__[1:])
 parser.add_argument('--version', action='version', version=pycbc.version.git_verbose_msg)
@@ -53,6 +53,10 @@ parser.add_argument('--ifar-threshold', type=float, default=None,
 parser.add_argument('--maximum-decisive-snr', type=float, default=None,
                     help="If given, only followup injections where the "
                          "decisive SNR is smaller than this value.")
+parser.add_argument('--nearby-triggers-window', type=float, default=0.05,
+                    help="Maximum time difference between the missed "
+                         "injection and the loudest SNR nearby trigger to "
+                         "display, seconds. Default=0.05")
 wf.add_workflow_command_line_group(parser)
 wf.add_workflow_settings_cli(parser, include_subdax_opts=True)
 args = parser.parse_args()
@@ -140,6 +144,50 @@ except KeyError:
     dec_chirp_dist = pycbc.pnutils.chirp_distance(dec_dist, mchirp)
     sorting = dec_chirp_dist.argsort()  # ascending order of dec chirp distance
 
+# Get the trigger SNRs and times
+# But only ones which are within a small window of the missed injection
+missed_inj_times = numpy.sort(f['injections/tc'][:][missed])
+
+# Note: Adding Earth diameter in light seconds to the window here
+# to allow for different IFO's arrival times of the injection
+safe_window = args.nearby_triggers_window + 0.0425
+
+def nearby_missedinj(endtime, snr):
+    """
+    Convenience function to check if trigger times are within a small
+    window of the injections
+
+    Parameters
+    ----------
+    endtime: numpy array
+        Trigger times to be checked against the missed injection times
+
+    snr: numpy array
+        Required by design of the HFile.select method but not used,
+        SNR of the triggers
+
+    Returns
+    -------
+    boolean array
+        True for triggers which are close to any missed injection
+    """
+    left = numpy.searchsorted(missed_inj_times - safe_window, endtime)
+    right = numpy.searchsorted(missed_inj_times + safe_window, endtime)
+    return left != right
+
+trigger_idx = {}
+trigger_snrs = {}
+trigger_times = {}
+for trig in single_triggers:
+    ifo = trig.ifo
+    with HFile(trig.lfn, 'r') as trig_f:
+        trigger_idx[ifo], trigger_times[ifo], trigger_snrs[ifo] = \
+            trig_f.select(
+                nearby_missedinj,
+                f'{ifo}/end_time',
+                f'{ifo}/snr',
+                return_indices=True)
+
 if len(missed) < num_events:
     num_events = len(missed)
 
@@ -178,7 +226,29 @@ for num_event in range(num_events):
                      (workflow, single_triggers, tmpltbank_file,
                       injection_file, args.output_dir, trig_id=trig_id,
                       file_substring='found_after_vetoes',
+                      title="Details of closest event",
                       tags=args.tags + [str(num_event)])[0],)]
+
+    for sngl in single_triggers:
+        # Find the triggers close to this injection at this IFO
+        ifo = sngl.ifo
+        trig_tdiff = abs(inj_params[ifo + '_end_time'] - trigger_times[ifo])
+        nearby = trig_tdiff < args.nearby_triggers_window
+        if not any(nearby):
+            # If there are no triggers in the defined window,
+            # do not print any info
+            continue
+        # Find the loudest SNR in this window
+        loudest = numpy.argmax(trigger_snrs[ifo][nearby])
+        # Convert to the indexin the trigger file
+        nearby_trigger_idx = trigger_idx[ifo][nearby][loudest]
+        # Make the info snippet
+        sngl_info = mini.make_sngl_ifo(workflow, sngl, tmpltbank_file,
+            nearby_trigger_idx, args.output_dir, ifo,
+            title=f"Parameters of loudest SNR nearby trigger in {ifo}",
+            tags=args.tags + [str(num_event)])[0]
+        layouts += [(sngl_info,)]
+
     files += mini.make_trigger_timeseries(workflow, single_triggers,
                               ifo_times, args.output_dir,
                               tags=args.tags + [str(num_event)])

bin/minifollowups/pycbc_page_coincinfo

@@ -41,16 +41,20 @@ parser.add_argument('--statmap-file-subspace-name', default='background_exc',
 trig_input = parser.add_mutually_exclusive_group(required=True)
 trig_input.add_argument('--n-loudest', type=int,
     help="Examine the n'th loudest trigger, use with statmap file")
-parser.add_argument('--sort-variable', default='ifar',
-                    help='Which subgroup of --analysis-category to use for '
-                         'sorting. Default=ifar')
-parser.add_argument('--sort-order', default='descending',
-                    choices=['ascending','descending'],
-                    help='Which direction to use when sorting on '
-                         '--sort-variable. Default=descending')
 trig_input.add_argument('--trigger-id', type=int,
     help="Examine the trigger with specified ID, use with statmap file. An "
          "alternative to --n-loudest. Cannot be used together")
+parser.add_argument('--sort-variable', default='ifar',
+    help='Which subgroup of --analysis-category to use for '
+         'sorting if using --n-loudest. Default=ifar')
+parser.add_argument('--sort-order', default='descending',
+    choices=['ascending','descending'],
+    help='Which direction to use when sorting on '
+         '--sort-variable with --n-loudest. Default=descending')
+parser.add_argument('--title',
+    help="Supply a title for the event details. Defaults are "
+         "'Parameters of event ranked N' if --n-loudest is given, or "
+         "'Details of trigger' for --trigger-id.")
 parser.add_argument('--include-summary-page-link', action='store_true',
     help="If given, will include a link to the DQ summary page on the "
          "single detector trigger tables.")
@@ -211,5 +215,5 @@ html += str(pycbc.results.static_table(data, headers))
 
 pycbc.results.save_fig_with_metadata(html, args.output_file, {},
                         cmd=' '.join(sys.argv),
-                        title=title,
+                        title=args.title if args.title else title,
                         caption=caption)

bin/minifollowups/pycbc_page_snglinfo

@@ -51,6 +51,9 @@ trig_args.add_argument("--n-loudest", type=int, default=None,
 trig_args.add_argument("--trigger-id", type=int, default=None,
     help="Use the trigger with this ID in the HDF file. Must supply either "
          "this option or --n-loudest.")
+parser.add_argument('--title',
+    help="Title for the produced snippet. If not given, will default to "
+         "pre-sets according to the options given")
 parser.add_argument('--include-summary-page-link', action='store_true',
     help="If given, will include a link to the DQ summary page")
 parser.add_argument('--include-gracedb-link', action='store_true',
@@ -194,6 +197,7 @@ html = pycbc.results.dq.redirect_javascript + \
         str(pycbc.results.static_table(data, headers))
 ###############################################################################
 
+# Set up default titles and the captions for the file
 if args.n_loudest:
     title = 'Parameters of single-detector event ranked %s' \
         % (args.n_loudest + 1)
@@ -204,5 +208,5 @@ else:
 
 pycbc.results.save_fig_with_metadata(html, args.output_file, {},
                         cmd = ' '.join(sys.argv),
-                        title = title,
+                        title = args.title if args.title else title,
                         caption = caption)

pycbc/workflow/minifollowups.py

@@ -527,7 +527,7 @@ def make_inj_info(workflow, injection_file, injection_index, num, out_dir,
 
 def make_coinc_info(workflow, singles, bank, coinc, out_dir,
                     n_loudest=None, trig_id=None, file_substring=None,
-                    sort_order=None, sort_var=None, tags=None):
+                    sort_order=None, sort_var=None, title=None, tags=None):
     tags = [] if tags is None else tags
     makedir(out_dir)
     name = 'page_coincinfo'
@@ -545,6 +545,8 @@ def make_coinc_info(workflow, singles, bank, coinc, out_dir,
         node.add_opt('--n-loudest', str(n_loudest))
     if trig_id is not None:
         node.add_opt('--trigger-id', str(trig_id))
+    if title is not None:
+        node.add_opt('--title', f'"{title}"')
     if file_substring is not None:
         node.add_opt('--statmap-file-subspace-name', file_substring)
     node.new_output_file_opt(workflow.analysis_time, '.html', '--output-file')
@@ -553,7 +555,7 @@ def make_coinc_info(workflow, singles, bank, coinc, out_dir,
     return files
 
 def make_sngl_ifo(workflow, sngl_file, bank_file, trigger_id, out_dir, ifo,
-                  tags=None):
+                  title=None, tags=None):
     """Setup a job to create sngl detector sngl ifo html summary snippet.
     """
     tags = [] if tags is None else tags
@@ -566,6 +568,8 @@ def make_sngl_ifo(workflow, sngl_file, bank_file, trigger_id, out_dir, ifo,
     node.add_input_opt('--bank-file', bank_file)
     node.add_opt('--trigger-id', str(trigger_id))
     node.add_opt('--instrument', ifo)
+    if title is not None:
+        node.add_opt('--title', f'"{title}"')
     node.new_output_file_opt(workflow.analysis_time, '.html', '--output-file')
     workflow += node
     files += node.output_files