#!/usr/bin/python3.13 -s

# 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_%s' % ifo] = antenna.effective_distance(
                                    data['distance'],
                                    data['ra'],
                                    data['dec'],
                                    data['polarization'],
                                    data['tc'],
                                    data['inclination']
                                    )
        eff_dist += 1.0 / data['eff_dist_%s' % ifo]
    data['eff_dist'] = 1.0 / eff_dist

    return data


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

    Parameters
    ----------
    hdf_file: str
        File path
    ifos: list
    snr_threshold: float
        NewSNR threshold
    bank_file: HFile object
    background_bestnrs: numpy.array, optional
        Used to compute FAP of quiet injections.

    Returns
    -------
    data: tuple of dictionaries
        Found and missed injection parameter dictionaries.
    """

    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 = {}
    logging.info('Loading injections...')

    # Load injections parameters
    for param in inj_params:
        missed_data[param] = inj_data['missed/%s' % param][...]
        found_data[param] = inj_data['found/%s' % 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_%s' % param] = \
            np.array(bank_file[param])[inj_data['network/template_id']]
    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/%s' % 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['%s/event_id' % ifo][...]):
            found_data['sigmasq_%s' % ifo] = inj_data['%s/sigmasq' % ifo][...]
            found_data['snr_%s' % ifo] = inj_data['%s/snr' % ifo][...]
        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['%s/event_id' % ifo][...]))
            found_data['sigmasq_%s' % ifo] = \
                inj_data['%s/sigmasq' % ifo][...][ifo_sorted_indices]
            found_data['snr_%s' % ifo] = \
                inj_data['%s/snr' % ifo][...][ifo_sorted_indices]
    # BestNRs
    found_data['bestnr'] = reweightedsnr_cut(found_data['reweighted_snr'][...],
                                             snr_threshold)
    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_%s' % ifo] == 0):
            logging.info("%s: sigmasq not set for at least one trigger.", ifo)
        if sum(found_data['sigmasq_%s' % 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_%s' % 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_%s' % ifo] for ifo in ifos])
    inj_sigma_tot = np.sum(inj_sigma_mult, axis=0)
    for ifo in ifos:
        found_data['inj_sigma_mean_%s' % ifo] = np.mean(
            found_data['sigmasq_%s' % ifo] * f_resp[ifo] / inj_sigma_tot)
    # Close the hdf file
    inj_data.close()

    return found_data, missed_data


# =============================================================================
# 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)
opts = parser.parse_args()

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 = []
# The user may want to process injections (possibly also against onsource)...
if 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 just triggers in offsource and onsource
elif onsource_file is not None:
    output_files = [lont_outfile, lont_h5_outfile, lofft_outfile,
                    lofft_h5_outfile]
    if None in output_files:
        parser.error('Please provide both on-source output files ' +
                     'and both off-source output files when using ' +
                     '--onsource-file.')
else:
    parser.error('Please provide --found-missed-file and/or --onsource-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 and vetoes
ifos, vetoes = ppu.extract_ifos_and_vetoes(offsource_file, opts.veto_files,
                                           opts.veto_category)

# Load triggers, time-slides, and segment dictionary
logging.info("Loading triggers.")
trig_data = ppu.load_triggers(offsource_file, ifos, None,
                              rw_snr_threshold=opts.newsnr_threshold)
logging.info("%d offsource triggers surviving reweighted SNR cut.",
             len(trig_data['network/event_id']))
logging.info("Loading timeslides.")
slide_dict = ppu.load_time_slides(offsource_file)
logging.info("Loading segments.")
segment_dict = ppu.load_segment_dict(offsource_file)

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

# Construct trials
logging.info("Constructing trials.")
trial_dict = ppu.construct_trials(opts.seg_files, segment_dict,
                                  ifos, slide_dict, vetoes)
total_trials = sum([len(trial_dict[slide_id]) for slide_id in slide_dict])
logging.info("%d trials generated.", total_trials)

# Extract basic trigger properties and store as dictionaries
trig_time, trig_snr, trig_bestnr = \
    ppu.extract_basic_trig_properties(trial_dict, trig_data, slide_dict,
                                      segment_dict, opts)
# Calculate SNR and BestNR values and maxima
time_veto_max_snr = {}
time_veto_max_bestnr = {}
for slide_id in slide_dict:
    num_slide_segs = len(trial_dict[slide_id])
    time_veto_max_snr[slide_id] = np.zeros(num_slide_segs)
    time_veto_max_bestnr[slide_id] = np.zeros(num_slide_segs)

for slide_id in slide_dict:
    for j, trial in enumerate(trial_dict[slide_id]):
        trial_cut = (trial[0] <= trig_time[slide_id])\
                          & (trig_time[slide_id] < trial[1])
        if not trial_cut.any():
            continue
        # Max SNR
        time_veto_max_snr[slide_id][j] = \
            max(trig_snr[slide_id][trial_cut])
        # Max BestNR
        time_veto_max_bestnr[slide_id][j] = \
            max(trig_bestnr[slide_id][trial_cut])
        # Max SNR for triggers passing SBVs
        sbv_cut = trig_bestnr[slide_id] != 0
        if not (trial_cut & sbv_cut).any():
            continue

logging.info("SNR and bestNR maxima calculated.")

# Output details of loudest offsouce triggers, sorted by BestNR
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_bestnr[slide_id],
                               sorted_trigs[slide_id]))

offsource_trigs.sort(key=lambda element: element[0])
offsource_trigs.reverse()

# Median and max values of SNR and BestNR
_, median_snr, _ = ppu.max_median_stat(slide_dict, time_veto_max_snr,
                                       trig_snr, total_trials)
max_bestnr, median_bestnr, full_time_veto_max_bestnr =\
    ppu.max_median_stat(slide_dict, time_veto_max_bestnr, trig_bestnr,
                        total_trials)

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['%s/event_id' % ifo] == 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
        num_trials_louder = 0
        for slide_id in slide_dict:
            for val in time_veto_max_bestnr[slide_id]:
                if val > bestnr:
                    num_trials_louder += 1
        fap = num_trials_louder/total_trials
        pval = '< %.3g' % (1./total_trials) if fap == 0 else '%.3g' % fap
        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]],
             mchirps[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['%s/snr' % ifo][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
    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(['%s SNR' % ifo for ifo in ifos])
    th.extend(['%s time shift (s)' % ifo for ifo in ifos])
    th.append('BestNR')

    # 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))
    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))

    # 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),
    #            full_time_veto_max_bestnr, delimiter='/t')


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

    # Get trigs
    on_trigs = ppu.load_triggers(onsource_file, ifos, None,
                                 rw_snr_threshold=opts.newsnr_threshold)

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

        # Gather bestNR index
        bestNR_event = np.argmax(on_trigs_bestnrs)
        loud_on_bestnr_trigs, loud_on_bestnr = \
            (on_trigs['network/event_id'][bestNR_event],
             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
    loud_on_fap = 1
    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['%s/event_id' % ifo] == trig_id)[0][0]
            for ifo in ifos
        }
        num_trials_louder = 0
        tot_off_snr = np.array([])
        for slide_id in slide_dict:
            num_trials_louder += sum(time_veto_max_bestnr[slide_id] >
                                     loud_on_bestnr)
            tot_off_snr = np.concatenate([tot_off_snr,
                                          time_veto_max_bestnr[slide_id]])
        fap = num_trials_louder/total_trials
        fap_test = sum(tot_off_snr > loud_on_bestnr)/total_trials
        pval = '< %.3g' % (1./total_trials) if fap == 0 else '%.3g' % fap
        loud_on_fap = fap
        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]],
             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['%s/snr' % ifo][ifo_trig_index[ifo]] for ifo in ifos] + \
            [loud_on_bestnr]
        td.append(d)
    else:
        td.append(["There are no events"] + [0 for number in range(11)] +
                  [0 for ifo in ifos] + [0])

    # 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'] + ['%s SNR' % ifo for ifo in ifos] + ['BestNR']

    td = list(zip(*td))

    # 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
    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)

else:
    tot_off_snr = np.array([])
    for slide_id in slide_dict:
        tot_off_snr = np.concatenate([tot_off_snr,
                                      time_veto_max_bestnr[slide_id]])
    med_snr = np.median(tot_off_snr)
    fap = sum(tot_off_snr > med_snr)/total_trials

# =======================
# Post-process injections
# =======================
if found_missed_file is not None:
    found_injs, missed_injs = load_missed_found_injections(
        found_missed_file, ifos, opts.newsnr_threshold, bank_data,
        background_bestnrs=full_time_veto_max_bestnr)
    logging.info("Missed/found injections/triggers loaded.")
    logging.info("%d found injections found.", len(found_injs['mchirp']))
    logging.info("%d missed injections found.", len(missed_injs['mchirp']))
    # Construct conditions for injection:
    # 1) found louder than background,
    zero_fap = found_injs['bestnr'] > max_bestnr

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

    # 3) missed after being recovered (i.e., vetoed)
    # -- > question: is there ever another way this happens other than veto?
    # vetoed_trigs = (~zero_fap) & (~nonzero_fap)
    vetoed_trigs = found_injs['bestnr'] == 0

    logging.info("%d found injections analysed.", len(found_injs['mchirp']))

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

    logging.info("%d missed injections analysed.", len(missed_injs['mchirp']))

    # Write quiet triggers to file
    sites = [ifo[0] for ifo in ifos]
    th = ['Dist'] + ['Eff. Dist. %s' % 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 %s' % ifo for ifo in ifos] +\
         ['BestNR', '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_%s' % ifo for ifo in ifos]
    keys = ['distance']
    keys += ['eff_dist_%s' % 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', 'spin1x', 'spin1y', 'spin1z', 'spin2x', 'spin2y',
             'spin2z', 'rec_spin1z', 'rec_spin2z']
    # The following parameters are available only for recovered injections
    na_keys = ['time_diff', 'rec_mass1', 'rec_mass2', 'rec_mchirp',
               'rec_spin1z', 'rec_spin2z', 'rec_ra', 'rec_dec', 'coherent_snr',
               'chisq', 'null_snr', 'bestnr']
    na_keys += sngl_snr_keys
    td = []
    for key in keys:
        if key in na_keys:
            td += [np.concatenate((found_injs[key][nonzero_fap],
                                   found_injs[key][vetoed_trigs],
                                   missed_na))]
        else:
            td += [np.concatenate((found_injs[key][nonzero_fap],
                                   found_injs[key][vetoed_trigs],
                                   missed_injs[key]))]
    td = list(zip(*td))
    td.sort(key=lambda elem: elem[0])
    td = list(zip(*td))

    # Write to h5 file
    logging.info("Writing %d quiet-found injections to h5 file.", len(td))
    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
    logging.info("Writing %d quiet-found injections to html file.",
                 len(td))
    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 to html file
    t_missed = []
    for key in keys:
        t_missed += [found_injs[key][vetoed_trigs]]
    t_missed = list(zip(*t_missed))
    t_missed.sort(key=lambda elem: elem[0])
    logging.info("Writing %d missed-found 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 found injections",
            'caption': "Recovered parameters and statistic values of \
            injections that are recovered, but downwieghted to BestNR = 0 \
            (i.e., vetoed).",
            '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
