#!/home/conda/feedstock_root/build_artifacts/bld/rattler-build_pycbc_1768936173/host_env_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehol/bin/python

# Copyright (C) 2021 Francesco Pannarale & Cameron Mills
#
# 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 3 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.

"""
Processes PyGRB triggers and injections to create html results tables.
"""

# =============================================================================
# Preamble
# =============================================================================
import sys
import os
import logging
import numpy as np

import pycbc.version
from pycbc.conversions import mchirp_from_mass1_mass2
from pycbc.detector import Detector
from pycbc.events.coherent import reweightedsnr_cut
from pycbc import init_logging
import pycbc.results
from pycbc.results import pygrb_postprocessing_utils as ppu
from pycbc.io.hdf import HFile

__author__ = "Francesco Pannarale <francesco.pannarale@ligo.org>"
__version__ = pycbc.version.git_verbose_msg
__date__ = pycbc.version.date
__program__ = "pycbc_pygrb_page_tables"


# =============================================================================
# Functions
# =============================================================================
def additional_injection_data(data, ifos):
    """Provides data with chirp masses and effective distances"""

    data['mchirp'] = mchirp_from_mass1_mass2(data['mass1'],
                                             data['mass2'])
    eff_dist = 0
    for ifo in ifos:
        antenna = Detector(ifo)
        data['eff_dist_'+ifo] = antenna.effective_distance(
                                    data['distance'],
                                    data['ra'],
                                    data['dec'],
                                    data['polarization'],
                                    data['tc'],
                                    data['inclination']
                                    )
        eff_dist += 1.0 / data['eff_dist_'+ifo]
    data['eff_dist'] = 1.0 / eff_dist

    return data


def load_missed_found_injections(hdf_file, ifos, bank_file, snr_threshold=None,
                                 background_bestnrs=None):
    """Loads found and missed injections from an hdf file as two dictionaries

    Parameters
    ----------
    hdf_file: str
        File path
    ifos: list
    bank_file: h5py.File object
    snr_threshold: float, optional [default: None]
        Reweighted SNR threshold
    background_bestnrs: numpy.array, optional [default: None]
        Used to compute FAP of quiet injections

    Returns
    -------
    data: tuple of dictionaries
        Found, missed, and missed after the cut in reweighted SNR injection
        parameter dictionaries.
    """

    logging.info('Loading injections...')
    inj_data = HFile(hdf_file, 'r')
    inj_params = ['mass1', 'mass2', 'distance', 'inclination', 'ra', 'dec',
                  'polarization', 'spin1x', 'spin1y', 'spin1z', 'spin2x',
                  'spin2y', 'spin2z', 'tc']
    found_data = {}
    # Missed injections (ones not recovered at all)
    missed_data = {}

    # Load injections parameters
    for param in inj_params:
        missed_data[param] = inj_data['missed/'+param][...]
        found_data[param] = inj_data['found/'+param][...]

    # Calculate effective distance for the ifos
    found_data = additional_injection_data(found_data, ifos)
    missed_data = additional_injection_data(missed_data, ifos)

    # Get recovered parameters and statistic values for the found injections
    # Recovered parameters
    for param in ['mass1', 'mass2', 'spin1z', 'spin2z']:
        found_data['rec_'+param] = \
            np.array(bank_file[param])[inj_data['network/template_id']]
    # If there are no found injections, simply move on
    # (there are no injections missed after the cut to return)
    if 'network/end_time_gc' not in inj_data.keys():
        return found_data, missed_data, {}
    # Otherwise carry on getting the recovered parameters and statistic values
    # of the found injections
    found_data['time_diff'] = \
        found_data['tc'] - inj_data['network/end_time_gc'][...]
    found_data['rec_mchirp'] = mchirp_from_mass1_mass2(
        found_data['rec_mass1'],
        found_data['rec_mass2'])
    # Recovered RA and Dec
    found_data['rec_ra'] = inj_data['network/ra'][...]
    found_data['rec_dec'] = inj_data['network/dec'][...]
    # Statistics values
    for param in ['coherent_snr', 'reweighted_snr', 'null_snr']:
        found_data[param] = inj_data['network/'+param][...]
    found_data['chisq'] = inj_data['network/my_network_chisq'][...]
    found_data['nifos'] = inj_data['network/nifo'][...].astype(int)
    for ifo in ifos:
        if np.all(inj_data['network/event_id'][...] ==
                  inj_data[ifo+'/event_id'][...]):
            found_data['sigmasq_'+ifo] = inj_data[ifo+'/sigmasq'][...]
            found_data['snr_'+ifo] = inj_data[ifo+'/snr'][...]
            found_data[ifo+'/end_time'] = inj_data[ifo+'/end_time'][...]
        else:
            # Sort the ifo event_id with respect to the network event_id
            ifo_sorted_indices = np.argsort(inj_data['network/event_id'][...][
                np.argsort(inj_data['network/event_id'])].searchsorted(
                    inj_data[ifo+'/event_id'][...]))
            found_data['sigmasq_'+ifo] = \
                inj_data[ifo+'/sigmasq'][...][ifo_sorted_indices]
            found_data['snr_'+ifo] = \
                inj_data[ifo+'/snr'][...][ifo_sorted_indices]
    # BestNRs
    found_data['bestnr'] = reweightedsnr_cut(found_data['reweighted_snr'][...],
                                             snr_threshold)
    # Apply reweighted SNR cut to the found injections
    cut_data = {}
    if snr_threshold:
        logging.info("%d found injections loaded.",
                     len(found_data[inj_params[0]]))
        logging.info("%d missed injections loaded.",
                     len(missed_data[inj_params[0]]))
        logging.info("Applying reweighted SNR cut at %s.", snr_threshold)
        rw_snr_cut = found_data['reweighted_snr'] < snr_threshold
        for k, v in found_data.items():
            cut_data[k] = v[rw_snr_cut]
            found_data[k] = v[~rw_snr_cut]
    del found_data['reweighted_snr']
    del cut_data['reweighted_snr']

    if background_bestnrs is not None:
        found_data['fap'] = np.array(
                [sum(background_bestnrs > bestnr) for bestnr in
                 found_data['bestnr']],
                dtype=float) / len(background_bestnrs)
    # Antenna responses
    f_resp = {}
    for ifo in ifos:
        if sum(found_data['sigmasq_'+ifo] == 0):
            logging.info("%s: sigmasq not set for at least one trigger.", ifo)
        if sum(found_data['sigmasq_'+ifo] != 0) == 0:
            logging.info("%s: sigmasq not set for any trigger.", ifo)
            if len(ifos) == 1:
                msg = "This is a single ifo analysis. "
                msg += "Setting sigmasq to unity for all triggers."
                logging.info(msg)
                found_data['sigmasq_'+ifo][:] = 1.0
        antenna = Detector(ifo)
        f_resp[ifo] = ppu.get_antenna_responses(antenna, found_data['ra'],
                                                found_data['dec'],
                                                found_data['tc'])

    inj_sigma_mult = \
        np.asarray([f_resp[ifo] *
                   found_data['sigmasq_'+ifo] for ifo in ifos])
    inj_sigma_tot = np.sum(inj_sigma_mult, axis=0)
    for ifo in ifos:
        found_data['inj_sigma_mean_'+ifo] = np.mean(
            found_data['sigmasq_'+ifo] * f_resp[ifo] / inj_sigma_tot)
    # Close the hdf file
    inj_data.close()

    logging.info("%d found injections.", len(found_data['mchirp']))
    logging.info("%d missed injections.", len(missed_data['mchirp']))
    logging.info("%d injections cut.", len(cut_data['mchirp']))

    return found_data, missed_data, cut_data


def format_pvalue_str(pvalue, n_trials):
    """Format p-value as a string."""
    return f'< {(1./n_trials):.3g}' if pvalue == 0 else f'{pvalue:.3g}'


# =============================================================================
# Main script starts here
# =============================================================================
parser = ppu.pygrb_initialize_plot_parser(description=__doc__)
parser.add_argument("-F", "--offsource-file", action="store", required=True,
                    help="Location of off-source trigger file")
parser.add_argument("--onsource-file", action="store",
                    help="Location of on-source trigger file.")
parser.add_argument("--found-missed-file", action="store",
                    help="HDF format file with injections to output " +
                    "details about.")
parser.add_argument("--num-loudest-off-trigs", action="store",
                    type=int, default=30, help="Number of loudest " +
                    "offsouce triggers to output details about.")
parser.add_argument("--bank-file", action="store", type=str, required=True,
                    help="Location of the full template bank used.")
parser.add_argument("--quiet-found-injs-output-file",
                    help="Quiet-found injections html output file.")
parser.add_argument("--missed-found-injs-output-file",
                    help="Missed-found injections html output file.")
parser.add_argument("--quiet-found-injs-h5-output-file",
                    help="Quiet-found injections h5 output file.")
parser.add_argument("--loudest-offsource-trigs-output-file",
                    help="Loudest offsource triggers html output file.")
parser.add_argument("--loudest-offsource-trigs-h5-output-file",
                    help="Loudest offsource triggers h5 output file.")
parser.add_argument("--loudest-onsource-trig-output-file",
                    help="Loudest onsource trigger html output file.")
parser.add_argument("--loudest-onsource-trig-h5-output-file",
                    help="Loudest onsource trigger h5 output file.")
parser.add_argument("-g", "--glitch-check-factor", action="store",
                    type=float, default=1.0, help="When deciding " +
                    "exclusion efficiencies this value is multiplied " +
                    "to the offsource around the injection trigger to " +
                    "determine if it is just a loud glitch.")
parser.add_argument("-C", "--cluster-window", action="store", type=float,
                    default=0.1, help="The cluster window used " +
                    "to cluster triggers in time.")
ppu.pygrb_add_bestnr_cut_opt(parser)
ppu.pygrb_add_slide_opts(parser)
opts = parser.parse_args()
ppu.slide_opts_helper(opts)

init_logging(opts.verbose, format="%(asctime)s: %(levelname)s: %(message)s")

# Store options used multiple times in local variables
offsource_file = opts.offsource_file
onsource_file = opts.onsource_file
found_missed_file = opts.found_missed_file
lofft_outfile = opts.loudest_offsource_trigs_output_file
lofft_h5_outfile = opts.loudest_offsource_trigs_h5_output_file
lont_outfile = opts.loudest_onsource_trig_output_file
lont_h5_outfile = opts.loudest_onsource_trig_h5_output_file
qf_outfile = opts.quiet_found_injs_output_file
mf_outfile = opts.missed_found_injs_output_file
qf_h5_outfile = opts.quiet_found_injs_h5_output_file

# Set output files and directories
output_files = []

# Check for correct input
if [found_missed_file, onsource_file].count(None) == 0:
    parser.error('Please provide --found-missed-file to process injections, ' +
                 '--onsource-file to process the on-source, or neither of ' +
                 'them to process the off-source triggers.')
# The user may process injections...
elif found_missed_file is not None:
    output_files = [qf_outfile, mf_outfile, qf_h5_outfile]
    if None in output_files:
        parser.error('Please provide all 3 injections output files when ' +
                     'using --found-missed-file')
# ...or triggers in the onsource...
elif onsource_file is not None:
    output_files = [lont_outfile, lont_h5_outfile]
    if None in output_files:
        parser.error('Please provide both on-source output files ' +
                     'when using --onsource-file.')
# ...or triggers in the offsource
# (both onsource_file and found_missed_file are None)
else:
    output_files = [lofft_outfile, lofft_h5_outfile]
    if None in output_files:
        parser.error('Please provide both off-source output files ' +
                     'when using --offsource-file.')
logging.info("Setting output directory.")
for output_file in output_files:
    if output_file:
        outdir = os.path.split(os.path.abspath(output_file))[0]
        if not os.path.isdir(outdir):
            os.makedirs(outdir)

# Extract IFOs
ifos = ppu.extract_ifos(offsource_file)

# Generate time-slides dictionary
slide_dict = ppu.load_time_slides(offsource_file)

# Generate segments dictionary
segment_dict = ppu.load_segment_dict(offsource_file)

# Construct trials removing vetoed times
trial_dict, total_trials = ppu.construct_trials(opts.seg_files, segment_dict,
                                                ifos, slide_dict,
                                                opts.veto_file)

# Load triggers (apply reweighted SNR cut, not vetoes)
trig_data = ppu.load_data(offsource_file, ifos, data_tag='offsource',
                          rw_snr_threshold=opts.newsnr_threshold,
                          slide_id=opts.slide_id)

# Extract needed trigger properties and store them as dictionaries
# Based on trial_dict: if vetoes were applied, trig_* are the veto survivors
# _av stands for after vetoes
keys = ['network/end_time_gc', 'network/coherent_snr',
        'network/reweighted_snr']
trig_data_av = ppu.extract_trig_properties(
    trial_dict,
    trig_data,
    slide_dict,
    segment_dict,
    keys
)

# Max SNR and BestNR values in each trial: these are stored in dictionaries
# keyed by slide_id, as arrays indexed by trial number
background_snr = {k: np.zeros(len(v)) for k, v in trial_dict.items()}
background = {k: np.zeros(len(v)) for k, v in trial_dict.items()}
for slide_id in slide_dict:
    trig_times = trig_data_av[keys[0]][slide_id]
    for j, trial in enumerate(trial_dict[slide_id]):
        # True whenever the trigger is in the trial
        trial_cut = (trial[0] <= trig_times) & (trig_times < trial[1])
        if not trial_cut.any():
            continue
        # Max SNR
        background_snr[slide_id][j] = \
            max(trig_data_av[keys[1]][slide_id][trial_cut])
        # Max BestNR
        background[slide_id][j] = \
            max(trig_data_av[keys[2]][slide_id][trial_cut])

# Max and median values of reweighted SNR,
# and sorted (loudest in trial) reweighted SNR values
max_bestnr, median_bestnr, sorted_bkgd =\
    ppu.max_median_stat(slide_dict, background,
                        trig_data_av[keys[2]], total_trials)
assert total_trials == len(sorted_bkgd)

# Median value of SNR
_, median_snr, _ = ppu.max_median_stat(slide_dict, background_snr,
                                       trig_data_av[keys[1]], total_trials)

logging.info("Background SNR and bestNR of trials calculated.")

# Output details of loudest offsouce triggers: only triggers compatible
# with the trial_dict are considered
offsource_trigs = []
sorted_trigs = ppu.sort_trigs(trial_dict, trig_data, slide_dict, segment_dict)
for slide_id in slide_dict:
    offsource_trigs.extend(
        zip(trig_data_av[keys[2]][slide_id], sorted_trigs[slide_id])
    )
offsource_trigs.sort(key=lambda element: element[0])
offsource_trigs.reverse()

# Calculate chirp masses of templates
logging.info('Loading triggers template masses')
bank_data = HFile(opts.bank_file, 'r')
template_mchirps = mchirp_from_mass1_mass2(
        bank_data['mass1'][...],
        bank_data['mass2'][...]
    )

# =========================================
# Output of loudest offsource triggers data
# =========================================
if lofft_outfile:
    # td: table data
    td = []

    # Gather properties of the loudest offsource triggers
    for i in range(min(len(offsource_trigs), opts.num_loudest_off_trigs)):
        bestnr = offsource_trigs[i][0]
        trig_id = offsource_trigs[i][1]
        trig_index = \
            np.where(trig_data['network/event_id'] == trig_id)[0][0]
        ifo_trig_index = {
            ifo: np.where(trig_data[ifo+'/event_id'] == trig_id)[0][0]
            for ifo in ifos
        }
        trig_slide_id = int(trig_data['network/slide_id'][trig_index])

        # Get trial of trigger, triggers with 'No trial' should have
        # already been removed!
        for j, trial in enumerate(trial_dict[trig_slide_id]):
            if trig_data['network/end_time_gc'][trig_index] in trial:
                chunk_num = j
                break
        else:
            chunk_num = 'No trial'

        # Get FAP of trigger
        pval = sum(sorted_bkgd > bestnr) / total_trials
        pval = format_pvalue_str(pval, total_trials)
        d = [chunk_num, trig_slide_id, pval,
             trig_data['network/end_time_gc'][trig_index],
             bank_data['mass1'][trig_data['network/template_id'][trig_index]],
             bank_data['mass2'][trig_data['network/template_id'][trig_index]],
             template_mchirps[trig_data['network/template_id'][trig_index]],
             bank_data['spin1z'][trig_data['network/template_id'][trig_index]],
             bank_data['spin2z'][trig_data['network/template_id'][trig_index]],
             trig_data['network/ra'][trig_index],
             trig_data['network/dec'][trig_index],
             trig_data['network/coherent_snr'][trig_index],
             trig_data['network/my_network_chisq'][trig_index],
             trig_data['network/null_snr'][trig_index]]
        d.extend([trig_data[ifo+'/snr'][ifo_trig_index[ifo]]
                  for ifo in ifos])
        d.extend([slide_dict[trig_slide_id][ifo] for ifo in ifos])
        d.append(bestnr)
        td.append(d)

    # th: table header [pycbc_pygrb_minifollowups looks for 'Slide Num']
    th = ['Trial', 'Slide Num', 'p-value', 'GPS time',
          'Rec. m1', 'Rec. m2', 'Rec. Mc', 'Rec. spin1z', 'Rec. spin2z',
          'Rec. RA', 'Rec. Dec', 'SNR', 'Chi^2', 'Null SNR']
    th.extend([ifo+' SNR' for ifo in ifos])
    th.extend([ifo+' time shift (s)' for ifo in ifos])
    th.append('BestNR')

    # When len(offsource_trigs) == 0, the loop above leaves td = [] unchanged
    # and this case needs to be handled adequately prior to moving on
    if not td:
        td = [[]] * len(th)

    # To ensure desired formatting in the h5 file and html table:
    # 1) "transpose" the data preserving its dtype
    td = list(zip(*td))

    # Write to h5 file
    logging.info("Writing %d loudest offsource triggers to h5 file.",
                 len(td[0]))
    lofft_h5_fp = HFile(lofft_h5_outfile, 'w')
    for i, key in enumerate(th):
        lofft_h5_fp.create_dataset(key, data=td[i])
    lofft_h5_fp.close()

    # Write to html file
    logging.info("Writing %d loudest triggers to html file.", len(td[0]))

    # To ensure desired formatting in the html table:
    # 2) convert the columns to numpy arrays
    # This is necessary as the p-values need to be treated as strings,
    # because they may contain a '<'
    td = [np.asarray(d) for d in td]

    # Format of table data
    format_strings = ['##.##', '##.##', None, '##.#####',
                      '##.##', '##.##', '##.##', '##.##', '##.##',
                      '##.##', '##.##', '##.##', '##.##', '##.##']
    format_strings.extend(['##.##' for ifo in ifos])
    format_strings.extend(['##.##' for ifo in ifos])
    format_strings.extend(['##.##'])
    html_table = pycbc.results.html_table(td, th,
                                          format_strings=format_strings,
                                          page_size=30)
    kwds = {'title': "Parameters of loudest offsource triggers",
            'caption': "Parameters of the " +
                       str(min(len(offsource_trigs),
                           opts.num_loudest_off_trigs)) +
                       " loudest offsource triggers.  " +
                       "The median reweighted SNR value is " +
                       str(median_bestnr) +
                       ".  The median SNR value is " +
                       str(median_snr),
            'cmd': ' '.join(sys.argv), }
    pycbc.results.save_fig_with_metadata(str(html_table),
                                         lofft_outfile, **kwds)

    # Store BestNR and FAP values: for collective FAP value studies at the
    # end of an observing run collectively
    # TODO: Needs a final place in the results webpage
    # np.savetxt('%s/bestnr_vs_fap_numbers.txt' %(outdir),
    #            sorted_bkgd, delimiter='/t')


# =======================
# Load on source triggers
# =======================
if onsource_file:

    # Get trigs
    on_trigs = ppu.load_data(onsource_file, ifos, data_tag=None,
                             rw_snr_threshold=opts.newsnr_threshold,
                             slide_id=0)

    # Record loudest trig by BestNR
    loud_on_bestnr = 0
    if on_trigs:
        on_trigs_bestnrs = on_trigs['network/reweighted_snr'][...]

        # Gather bestNR index
        if on_trigs_bestnrs.size > 0:
            bestNR_event = np.argmax(on_trigs_bestnrs)
            loud_on_bestnr_trigs = on_trigs['network/event_id'][bestNR_event]
            loud_on_bestnr = on_trigs_bestnrs[bestNR_event]
    # If the loudest event has bestnr = 0, there is no event at all!
    if loud_on_bestnr == 0:
        loud_on_bestnr_trigs = None

    logging.info("Onsource analysed.")

    # Table data
    td = []

    # Gather data
    if loud_on_bestnr_trigs:
        trig_id = loud_on_bestnr_trigs
        trig_index = np.where(on_trigs['network/event_id'] == trig_id)[0][0]
        ifo_trig_index = {
            ifo: np.where(on_trigs[ifo+'/event_id'] == trig_id)[0][0]
            for ifo in ifos
        }
        num_trials_louder = 0
        pval = sum(sorted_bkgd > loud_on_bestnr)/total_trials
        pval = format_pvalue_str(pval, total_trials)
        d = [pval,
             on_trigs['network/end_time_gc'][trig_index],
             bank_data['mass1'][on_trigs['network/template_id'][trig_index]],
             bank_data['mass2'][on_trigs['network/template_id'][trig_index]],
             template_mchirps[on_trigs['network/template_id'][trig_index]],
             bank_data['spin1z'][on_trigs['network/template_id'][trig_index]],
             bank_data['spin2z'][on_trigs['network/template_id'][trig_index]],
             on_trigs['network/ra'][trig_index],
             on_trigs['network/dec'][trig_index],
             on_trigs['network/coherent_snr'][trig_index],
             on_trigs['network/my_network_chisq'][trig_index],
             on_trigs['network/null_snr'][trig_index]] + \
            [on_trigs[ifo+'/snr'][ifo_trig_index[ifo]] for ifo in ifos] + \
            [loud_on_bestnr]
        td.append(d)

    # Table header
    th = ['p-value', 'GPS time', 'Rec. m1', 'Rec. m2', 'Rec. Mc',
          'Rec. spin1z', 'Rec. spin2z', 'Rec. RA', 'Rec. Dec', 'SNR', 'Chi^2',
          'Null SNR'] + [ifo+' SNR' for ifo in ifos] + ['BestNR']

    td = list(zip(*td))

    # Handle the case in which there is no data to be placed in the table
    if not td:
        td = [[]] * len(th)

    # Write to h5 file
    logging.info("Writing loudest onsource trigger to h5 file.")
    with HFile(lont_h5_outfile, 'w') as lont_h5_fp:
        for i, key in enumerate(th):
            lont_h5_fp.create_dataset(key, data=td[i])

    # Write to html file
    logging.info("Writing loudest onsource trigger to html file.")

    # Format of table data
    format_strings = [None, '##.#####', '##.##', '##.##', '##.##', '##.##',
                      '##.##', '##.##', '##.##', '##.##', '##.##', '##.##']
    format_strings.extend(['##.##' for ifo in ifos])
    format_strings.extend(['##.##'])

    # Table data: assemble human readable message when no trigger is recovered
    if not loud_on_bestnr_trigs:
        td = [list("-" * len(format_strings))]
        td[0][0] = "There are no events"
    td = [np.asarray(d) for d in td]
    html_table = pycbc.results.html_table(td, th,
                                          format_strings=format_strings,
                                          page_size=1)
    kwds = {'title': "Loudest event",
            'caption': "Recovered parameters and statistic values of the \
            loudest trigger.",
            'cmd': ' '.join(sys.argv), }
    pycbc.results.save_fig_with_metadata(str(html_table), lont_outfile,
                                         **kwds)

# =======================
# Post-process injections
# =======================
if found_missed_file is not None:
    # Load injections applying reweighted SNR cut
    found_injs, missed_injs, cut_injs = load_missed_found_injections(
        found_missed_file, ifos, bank_data,
        snr_threshold=opts.newsnr_threshold,
        background_bestnrs=sorted_bkgd
    )

    # Split in injections found surviving vetoes and ones found but vetoed
    found_after_vetoes, vetoed, *_ = ppu.apply_vetoes_to_found_injs(
        found_missed_file,
        found_injs,
        ifos,
        veto_file=opts.veto_file
    )

    if 'bestnr' not in found_after_vetoes:
        found_after_vetoes['bestnr'] = np.array([])

    # Construct conditions for injection:
    # 1) found louder than background,
    zero_fap = found_after_vetoes['bestnr'] > max_bestnr

    # 2) found (bestnr > 0) but not louder than background (non-zero FAP)
    nonzero_fap = ~zero_fap & (found_after_vetoes['bestnr'] != 0)

    # 3) missed after being recovered: vetoed (these have bestnr = 0)

    # Avoids a problem with formatting in the non-static html output file
    # missed_na = [-0] * len(missed_injs['mchirp'])

    # Write quiet triggers to file
    sites = [ifo[0] for ifo in ifos]
    th = ['Dist'] + ['Eff. Dist. '+site for site in sites] +\
         ['GPS time', 'GPS time - Rec. Time'] +\
         ['Inj. m1', 'Inj. m2', 'Inj. Mc', 'Rec. m1', 'Rec. m2', 'Rec. Mc',
          'Inj. inc', 'Inj. RA', 'Inj. Dec', 'Rec. RA', 'Rec. Dec', 'SNR',
          'Chi^2', 'Null SNR'] +\
         ['SNR '+ifo for ifo in ifos] +\
         ['BestNR', 'p-value',
          'Inj S1x', 'Inj S1y', 'Inj S1z',
          'Inj S2x', 'Inj S2y', 'Inj S2z',
          'Rec S1z', 'Rec S2z']
    # Format of table data
    format_strings = ['##.##']
    format_strings.extend(['##.##' for ifo in ifos])
    format_strings.extend(['##.#####', '##.#####',
                           '##.##', '##.##', '##.##',
                           '##.##', '##.##', '##.##',
                           '##.##', '##.##', '##.##',
                           '##.##', '##.##', '##.##',
                           '##.##', '##.##'])
    format_strings.extend(['##.##' for ifo in ifos])
    format_strings.extend(['##.##', '#.######',
                           '##.##', '##.##', '##.##',
                           '##.##', '##.##', '##.##',
                           '##.##', '##.##'])
    sngl_snr_keys = ['snr_'+ifo for ifo in ifos]
    keys = ['distance']
    keys += ['eff_dist_'+ifo for ifo in ifos]
    keys += ['tc', 'time_diff', 'mass1', 'mass2', 'mchirp', 'rec_mass1',
             'rec_mass2', 'rec_mchirp', 'inclination', 'ra', 'dec', 'rec_ra',
             'rec_dec', 'coherent_snr', 'chisq', 'null_snr']
    keys += sngl_snr_keys
    keys += ['bestnr', 'fap', 'spin1x', 'spin1y', 'spin1z', 'spin2x', 'spin2y',
             'spin2z', 'rec_spin1z', 'rec_spin2z']
    for key in keys:
        if key not in found_after_vetoes:
            found_after_vetoes[key] = np.array([])
    td = [found_after_vetoes[key][nonzero_fap] for key in keys]
    td = list(zip(*td))
    td.sort(key=lambda elem: elem[0])
    logging.info("Writing %d quiet-found injections to h5 and html files.",
                 len(td))
    td = list(zip(*td))

    # Handle the case in which there is no data to be placed in the table
    if not td:
        td = [[]] * len(th)

    # Write to h5 file
    with HFile(qf_h5_outfile, 'w') as qf_h5_fp:
        for i, key in enumerate(th):
            qf_h5_fp.create_dataset(key, data=td[i])

    # Write to html file
    td = [np.asarray(d) for d in td]
    html_table = pycbc.results.html_table(td, th,
                                          format_strings=format_strings,
                                          page_size=20)
    kwds = {'title': "Quiet found injections",
            'caption': "Recovered parameters and statistic values of \
            injections that are recovered, but not louder than \
            background.", 'cmd': ' '.join(sys.argv), }
    pycbc.results.save_fig_with_metadata(str(html_table), qf_outfile,
                                         **kwds)

    # Write quiet triggers to html file
    if len(cut_injs) == 0:
        cut_injs = dict.fromkeys(keys, np.array([]))
    for key in keys:
        if key not in vetoed:
            vetoed[key] = np.full(len(vetoed['mass1']), -1)
        if key not in cut_injs:
            cut_injs[key] = np.full(len(cut_injs['mass1']), -1)
        if key not in missed_injs:
            missed_injs[key] = np.full(len(missed_injs['mass1']), -1)
    # Update the data to include a label stating whether the injection was
    # vetoed, found and then removed due to the reweighted SNR cut, or missed
    vetoed['category'] = np.full(len(vetoed['mass1']), 'Vetoed')
    cut_injs['category'] = np.full(len(cut_injs['mass1']), 'Cut')
    missed_injs['category'] = np.full(len(missed_injs['mass1']), 'Missed')
    # Update the header and table formatter
    th += ['Category']
    format_strings.extend(['##'])
    t_missed = [np.concatenate((missed_injs[key], vetoed[key], cut_injs[key]))
                for key in keys+['category']]
    t_missed = list(zip(*t_missed))
    t_missed.sort(key=lambda elem: elem[0])
    logging.info("Writing %d missed, vetoed, and cut injections to html file.",
                 len(t_missed))
    t_missed = zip(*t_missed)
    t_missed = [np.asarray(d) for d in t_missed]
    html_table = pycbc.results.html_table(t_missed, th,
                                          format_strings=format_strings,
                                          page_size=20)
    kwds = {'title': "Missed, vetoed, and cut injections",
            'caption': "Parameters of missed injections and \
            recovered parameters and statistic values of \
            injections that are recovered but then vetoed or cut because \
            their reweighted SNR is below threshold.",
            'cmd': ' '.join(sys.argv), }
    pycbc.results.save_fig_with_metadata(str(html_table), mf_outfile,
                                         **kwds)

# Close the bank file
bank_data.close()

# Post-processing of injections ends here
