path: root/src/dc.c

/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
 *
 * 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, see <https://www.gnu.org/licenses/>.
 */

#include "dc.h"
#include "event.h"
#include "sort.h"
#include <stddef.h> /* for size_t */
#include "misc.h"
#include "pmt.h"
#include "vector.h"
#include <gsl/gsl_sort.h>
#include "zdab_utils.h"
#include "zebra.h"
#include <gsl/gsl_statistics_double.h>
#include "Record_Info.h"

#define MAX_PAIRS 10000

/* Returns 1 if the event is tagged as failing the ITC cut. The definition for
 * the ITC cut comes from the SNOMAN companion:
 *
 *     This tags events having fewer than 60% of the tubes within 93 ns
 *     coincidence ("in-time"). The number of in-time hits is the maximum
 *     number in any position of a 93 ns sliding window applied to the PMT time
 *     spectrum of the event.
 *
 */
int is_itc(event *ev)
{
    size_t i, j;
    double t_array[MAX_PMTS];
    int nhit, n, maxn;
    double start;

    nhit = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        /* Only loop over hit PMTs with good calibrations. */
        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].flags) continue;

        t_array[nhit++] = ev->pmt_hits[i].t;
    }

    /* Sort the times. */
    gsl_sort(t_array,1,nhit);

    maxn = 0;
    for (i = 0; i < nhit; i++) {
        start = t_array[i];
        n = 1;
        for (j = i+1; j < nhit; j++) {
            if (t_array[j] < start + ITC_TIME_WINDOW) n += 1;
            else break;
        }
        if (n > maxn) maxn = n;
    }

    return maxn/(double) nhit < ITC_TIME_FRACTION;
}

/* Returns 1 if the event is tagged as failing the FTS cut. The definition for
 * the FTS cut comes from the SNOMAN companion:
 *
 *     This tags events where the median time difference of hit PMT pairs
 *     within 3 m of each other is greater than 6.8 ns. Pairs with large time
 *     differences (> 25 ns) are not included in the median and at least 15
 *     pairs satisfying the distance and time difference criteria are required
 *     before the event can be tagged.
 *
 * I used the SNOMAN implementation in flt_fts_cut.for for reference.
 *
 */
int is_fts(event *ev)
{
    size_t i, j;
    double pmt_dir[3], distance, dt;
    double deltat_array[MAX_PAIRS];
    size_t count;
    double dtmedian;

    count = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        /* Only loop over hit PMTs with good calibrations. */
        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].flags) continue;

        for (j = i+1; j < MAX_PMTS; j++) {
            /* Only loop over hit PMTs with good calibrations. */
            if (i == j || !ev->pmt_hits[j].hit || ev->pmt_hits[j].flags) continue;

            /* Calculate distance between PMTs. */
            SUB(pmt_dir,pmts[i].pos,pmts[j].pos);
            distance = NORM(pmt_dir);

            /* Get the absolute value of the time difference. */
            dt = fabs(ev->pmt_hits[i].t - ev->pmt_hits[j].t);

            if (dt <= FTS_DT_THRESH && distance < FTS_DIST_THRESH) {
                deltat_array[count++] = dt;
                if (count >= MAX_PAIRS) goto end;
            }
        }
    }

end:

    if (count > FTS_COUNT_THRESH) {
        gsl_sort(deltat_array,1,count);
        dtmedian = gsl_stats_median_from_sorted_data(deltat_array,1,count);
        return dtmedian > FTS_MEDIAN_CUT;
    }

    return 0;
}

/* Returns 1 if the event is tagged as an OWL trigger event. The definition for
 * the OWL trigger cut comes from the SNOMAN companion:
 *
 *     This tags events where the OWLESUMHI trigger fires.
 *
 */
int is_owl_trigger(event *ev)
{
    return (ev->trigger_type & TRIG_OWLE_HI) != 0;
}

/* Returns 1 if the event is tagged as an OWL event. The definition for the OWL
 * cut comes from the SNOMAN companion:
 *
 *     This tags events where the number of OWLs + BUTTs is 3 or more.
 *
 */
int is_owl(event *ev)
{
    size_t i;
    size_t nhit;

    nhit = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (ev->pmt_hits[i].hit && (pmts[i].pmt_type == PMT_OWL || pmts[i].pmt_type == PMT_BUTT))
            nhit += 1;
    }

    return nhit >= 3;
}

/* Returns 1 if the event fails the ESUM cut. The definition for the ESUM cut
 * comes from the SNOMAN companion:
 *
 *     Tags an event if the ESUMLO AND ESUMHI trigger bits are set, AND NONE of
 *     the following trigger bits are set: (NHIT100LO, NHIT100MED,
 *     NHIT100HI,NHIT20, NHIT20LB, OWLN, OWLELO, OWLEHI, PULSE_GT, PRESCALE,
 *     MISS).
 *
 */
int is_esum(event *ev)
{
    return (ev->trigger_type & (TRIG_ESUM_LO | TRIG_ESUM_HI | TRIG_NHIT_100_LO | TRIG_NHIT_100_MED | TRIG_NHIT_100_HI | TRIG_NHIT_20 | TRIG_NHIT_20_LB | TRIG_OWLN | TRIG_OWLE_LO | TRIG_OWLE_HI | TRIG_PULSE_GT | TRIG_PRESCALE | TRIG_MISS)) == (TRIG_ESUM_LO | TRIG_ESUM_HI);
}

/* Returns the data cleaning bitmask for the event `ev`. */
uint32_t get_dc_word(event *ev, zebraFile *f, zebraBank *bmast, zebraBank *bev)
{
    uint32_t status = 0x0;

    if (is_muon(ev))
        status |= DC_MUON;
    if (junk_cut(f, bmast, bev))
        status |= DC_JUNK;
    if (crate_isotropy(ev))
        status |= DC_CRATE_ISOTROPY;
    if (qvnhit(ev))
        status |= DC_QVNHIT;
    if (is_neck_event(ev))
        status |= DC_NECK;
    if (is_flasher(ev))
        status |= DC_FLASHER;
    if (is_esum(ev))
        status |= DC_ESUM;
    if (is_owl(ev))
        status |= DC_OWL;
    if (is_owl_trigger(ev))
        status |= DC_OWL_TRIGGER;
    if (is_fts(ev))
        status |= DC_FTS;
    if (is_itc(ev))
        status |= DC_ITC;
    if (is_breakdown(ev))
        status |= DC_BREAKDOWN;

    return status;
}

/* Returns 1 if the event is tagged as an incoming muon. The goal of this cut
 * is to tag external muons *without* tagging muons which start within the PSUP
 * and then exit. In order to achieve that we use the ECA time of the OWL
 * tubes.
 *
 * Ideally we would have the full calibrated time of the OWL tubes, but there
 * was never any PCA for the OWL tubes. I also discussed this with Josh and he
 * said that the cables were never changed for the OWL tubes so we should at
 * least be able to rely on them being consistent.
 *
 * Since the average PCA correction is of order +/- 5 ns and the round trip
 * time for light from the edge of the PSUP to the cavity and back should be
 * longer than 15 ns (4 meters), the ECA time should be enough to distinguish
 * incoming vs outgoing muons.
 *
 * This function tags an event as an incoming muon if the following criteria
 * are met:
 *
 * - calibrated nhit >= 150
 * - at least 5 calibrated OWL hits
 * - at least 1 OWL hit with a QHS >= 20.0
 * - at least 2 OWL hits with ECA times earlier than the 10th percentile of the
 *   PSUP PMT times
 *
 * Idea: Should we also require at least one of the early tubes to have a high
 * charge? */
int is_muon(event *ev)
{
    size_t i;
    double ept[MAX_PMTS];
    size_t nhit, nhit_owl, nhit_owl_early;
    double early_time, max_owl_ehs;

    nhit = 0;
    nhit_owl = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        /* Require good calibrations.
         *
         * Note: Typically I would do this with the following check:
         *
         *     ev->pmt_hits[i].pf & (KPF_NO_CAL | KPF_BAD_CAL)
         *
         * However, Javi mentioned in an email that for very early hits which
         * are in the TAC curl region, the calibration processor will flag
         * them, set the time to -9999.0, and set the KPF_BAD_CAL bit. However,
         * in this case for some reason the ECA time still looks reasonable.
         * Therefore, since we only use the ECA time here we instead just make
         * sure that the ECA time is something reasonable. */
        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue;

        switch (pmts[i].pmt_type) {
        case PMT_OWL:
            if (nhit_owl == 0) {
                max_owl_ehs = ev->pmt_hits[i].ehs;
            } else if (ev->pmt_hits[i].ehs > max_owl_ehs) {
                max_owl_ehs = ev->pmt_hits[i].ehs;
            }
            nhit_owl += 1;
            break;
        case PMT_NORMAL:
            if (!ev->pmt_hits[i].flags) ept[nhit++] = ev->pmt_hits[i].ept;
            break;
        }
    }

    /* Require at least 150 normal hits. */
    if (nhit < MUON_MIN_NHIT) return 0;

    /* Require at least 5 OWL hits. */
    if (nhit_owl < MUON_MIN_OWL_NHIT) return 0;

    /* Require at least 1 OWL tube to have a QHS >= 20.0. */
    if (max_owl_ehs < MUON_OWL_QHS) return 0;

    /* Sort the times. */
    gsl_sort(ept,1,nhit);

    /* Calculate the time to define early OWL hits. */
    early_time = gsl_stats_quantile_from_sorted_data(ept,1,nhit,MUON_OWL_TIME_PERCENTILE);

    nhit_owl_early = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        /* Require good calibrations. */
        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue;

        if (pmts[i].pmt_type == PMT_OWL && ev->pmt_hits[i].ept < early_time) nhit_owl_early += 1;
    }

    return nhit_owl_early >= 2;
}

/* Returns 1 if the event fails the junk cut. The definition for the junk cut
 * comes from the SNOMAN companion:
 *
 *     This will remove orphans, ECA data and events containing the same PMT
 *     more than once.
 *
 * I also added a check that will flag the event if it contains more than
 * MAX_PMTS PMT hits (which shouldn't be possible).
 *
 * Note: This function may return -1 if there is an error traversing the PMT
 * banks. */
int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev)
{
    int i;
    int ids[MAX_PMTS];
    int nhit;
    static int pmt_links[] = {KEV_PMT,KEV_OWL,KEV_LG,KEV_FECD,KEV_BUTT,KEV_NECK};
    static char *pmt_names[] = {"PMT","OWL","LG","FECD","BUTT","NECK"};
    size_t index[MAX_PMTS];
    PMTBank bpmt;
    zebraBank b;
    int rv;

    /* Fail orphans. */
    if (bev->status & KEV_ORPHAN) return 1;

    /* Fail any events with the ECA bit set. */
    if (bmast->status & KMAST_ECA) return 1;

    /* Next we check to see if any PMTs got hit more than once. To do this we
     * construct an array of the PMT PIN numbers, sort it, and then check for
     * duplicates. */

    nhit = 0;
    for (i = 0; i < LEN(pmt_links); i++) {
        if (bev->links[pmt_links[i]-1] == 0) continue;

        rv = zebra_get_bank(f,&b,bev->links[pmt_links[i]-1]);

        if (rv) {
            fprintf(stderr, "error getting %s bank: %s\n", pmt_names[i], zebra_err);
            return -1;
        }

        while (1) {
            unpack_pmt(b.data, &bpmt);

            ids[nhit++] = bpmt.pin;

            if (nhit >= MAX_PMTS) {
                fprintf(stderr, "event with more than %i PMT hits\n", MAX_PMTS);
                return 1;
            }

            if (!b.next) break;

            rv = zebra_get_bank(f,&b,b.next);

            if (rv) {
                fprintf(stderr, "error getting %s bank: %s\n", pmt_names[i], zebra_err);
                return -1;
            }
        }
    }

    argsort(ids,nhit,index);

    for (i = 0; i < nhit - 1; i++) {
        /* Check if we have duplicate PMT hits. */
        if (ids[i] == ids[i+1]) return 1;
    }

    return 0;
}

/* Returns 1 if the event fails the crate isotropy test. The definition of
 * this cut comes from the SNOMAN companion:
 *
 *     This tags events where more than 70% of the hits occur on one crate and
 *     more than 80% of those hits occur on two adjacent cards, including a
 *     wrap around to count cards 0 and 15 together. 
 *
 */
int crate_isotropy(event *ev)
{
    size_t i;
    int crate_nhit[20];
    int card_nhit[16];
    int nhit;
    int crate, card;
    int max_crate;

    for (i = 0; i < LEN(crate_nhit); i++) crate_nhit[i] = 0;

    /* Loop over all the hits and add one to the crate_nhit array for each hit
     * in a crate. */

    nhit = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (!ev->pmt_hits[i].hit) continue;

        crate = i/512;

        crate_nhit[crate] += 1;

        nhit += 1;
    }

    /* Check if any crate has more than 70% of the hits. */
    max_crate = -1;
    for (i = 0; i < LEN(crate_nhit); i++) {
        if (crate_nhit[i] > 7*nhit/10) max_crate = i;
    }

    /* If no crates had more than 70% of the hits, return 0. */
    if (max_crate == -1) return 0;

    for (i = 0; i < LEN(card_nhit); i++) card_nhit[i] = 0;

    /* Now we count up the number of hits in adjacent cards. We store the
     * number of hits in the card_nhit array as follows: Hits from card 15 or
     * card 0 will be added to card_nhit[0], hits from cards 0 or 1 will go
     * into card_nhit[1], etc. */

    nhit = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (!ev->pmt_hits[i].hit) continue;

        crate = i/512;

        if (crate != max_crate) continue;

        card = (i % 512)/32;

        card_nhit[card] += 1;
        card_nhit[(card + 1) % LEN(card_nhit)] += 1;

        nhit += 1;
    }

    /* Check to see if 80% of the hits in this crate came from two adjacent
     * cards. */
    for (i = 0; i < LEN(card_nhit); i++) {
        if (card_nhit[i] > 8*nhit/10) return 1;
    }

    return 0;
}

/* Returns 1 if the event is tagged by the QvNHIT cut. The definition of
 * this cut comes from the SNOMAN companion:
 *
 *     Using pedestal subtracted charge and a hardwired gain of 32.3 first
 *     throw away the 10% of the tubes that have the largest charge and then
 *     divide this by 0.9* NHIT. Tag the event if this ratio is below 0.25.
 *
 * I copied the logic from the SNOMAN file flt_q_nhit_cut.for. */
int qvnhit(event *ev)
{
    size_t i;
    double ehl[MAX_PMTS];
    int nhit, max;
    double sum, qhl_ratio;

    nhit = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;

        /* Require good calibrations. */
        if (ev->pmt_hits[i].pf & (KPF_NO_CAL | KPF_BAD_CAL)) continue;

        /* Remove unphysical charge hits. */
        if (ev->pmt_hits[i].qhl <= -100) continue;

        /* FIXME: SNOMAN code checks for pmt terminator. I don't know where
         * that info is stored. */

        /* Note: We use QHL instead of QHS here because that's how SNOMAN did
         * it. In the SNOMAN code there is a comment which says:
         *
         *     we use QHL since it is a better measure of pickup (long integrate)
         *
         */
        if (ev->pmt_hits[i].hit) {
            ehl[nhit++] = ev->pmt_hits[i].ehl/32.3;
        }
    }

    gsl_sort(ehl,1,nhit);

    max = nhit*9/10;

    /* Check that we have enough calibrated tubes. */
    if (max < 5) return 0;

    sum = 0.0;
    for (i = 0; i < max; i++) {
        sum += ehl[i];
    }

    qhl_ratio = sum/max;

    return qhl_ratio < QRATIO_THRESHOLD;
}

/* Returns the time difference between the neck PMT hit and the average ECA
 * time of the PSUP PMTs required to tag an event as a neck event.
 *
 * These values come from the neck tube cut table in prod/filter_flth.dat. */
double get_neck_tube_cut_time_diff(event *ev)
{
    if (ev->dte < 19991216 || (ev->dte == 19991216 && ev->hmsc < 20370000))
        return 70.0;
    else if (ev->dte < 20040923 || (ev->dte < 20040923 && ev->hmsc < 100222))
        return 15.0;
    else
        return -85.0;
}

/* Returns 1 if the event is classified as a neck event. The definition of
 * a neck event is a combination of the neck flag in the DAMN word:
 *
 *     This cuts events containing neck tubes. It requires that either both
 *     tubes in the neck fire, or that one of those tubes fires and it has a
 *     high charge and is early. High charge is defined by a neck tube having a
 *     pedestal subtracted charge greater than 70 or less than -110. Early if
 *     defined by the neck tube having an ECA time 70ns or more before the
 *     average ECA time of the PSUP PMTS with z les than 0. After the cable
 *     changes to the neck tubes this time difference changes to 15ns. 
 *
 * and the added requirement that 50% of the normal PMTs have z <= -425.0 or
 * 50% of the ECA calibrated charge for normal PMTs have z <= -425.0.
 *
 * I added the requirement that 50% of the normal PMTs be near the bottom since
 * for the high energy events I'm interested in, it's entirely possible that a
 * high energy electron or muon travelling upwards could cause multiple neck
 * PMTs to fire.
 *
 */
int is_neck_event(event *ev)
{
    size_t i;
    int n;
    int high_charge;
    double avg_ept_psup;
    int nhit, nhit_bottom;
    double time_diff;
    double ehs, ehs_bottom;

    /* First, we compute the average ECA time for all the PSUP PMTs with z <= 0. */

    nhit = 0;
    avg_ept_psup = 0.0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue;

        /* Require good calibrations. */
        if (ev->pmt_hits[i].pf & (KPF_NO_CAL | KPF_BAD_CAL)) continue;

        if (ev->pmt_hits[i].hit && pmts[i].pos[2] < 0) {
            avg_ept_psup += ev->pmt_hits[i].ept;
            nhit += 1;
        }
    }

    if (nhit < 1) return 0;

    avg_ept_psup /= nhit;

    time_diff = get_neck_tube_cut_time_diff(ev);

    /* Now, we check if two or more neck tubes fired *or* a single neck tube
     * fired with a high charge and 70 ns before the average time of the normal
     * PMTs. */
    n = 0;
    high_charge = 0;
    nhit = 0;
    nhit_bottom = 0;
    ehs = 0;
    ehs_bottom = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (!ev->pmt_hits[i].hit) continue;

        switch (pmts[i].pmt_type) {
        case PMT_NECK:
            /* Require good calibrations. */
            if (ev->pmt_hits[i].ept <= -100) break;

            n += 1;

            /* FIXME: The SNOMAN companion says "Early if defined by the neck
             * tube having an ECA time 70ns or more before the average ECA time
             * of the PSUP PMTS with z les than 0." When does this change take
             * place? */
            if ((ev->pmt_hits[i].ehs > 70 || ev->pmt_hits[i].ehs < -110) &&
                (ev->pmt_hits[i].ept < avg_ept_psup - time_diff)) high_charge = 1;
            break;
        case PMT_NORMAL:
            if (pmts[i].pos[2] < NECK_BOTTOM_HEIGHT) {
                nhit_bottom += 1;
                ehs_bottom += ev->pmt_hits[i].ehs;
            }
            nhit += 1;
            ehs += ev->pmt_hits[i].ehs;
            break;
        }
    }

    if ((n >= 2 || high_charge) && ((nhit_bottom >= nhit*NECK_BOTTOM_FRACTION) || (ehs_bottom >= ehs*NECK_BOTTOM_FRACTION))) return 1;

    return 0;
}

/* Returns 1 if the event is classified as a breakdown event. The definition
 * for a breakdown event is:
 *
 *   - nhit > 1000
 *
 *   - a single crate has at least 256 hits and the median TAC for that crate
 *     is more than 500 from the median TAC for all the other crates with at
 *     least 20 hits
 *
 */
int is_breakdown(event *ev)
{
    int i;
    double tac[20][512];
    int nhit[20] = {0};
    double median_tac[20];
    size_t index[20];

    /* Breakdown event must have an nhit greater than 1000. */
    if (ev->nhit < 1000) return 0;

    for (i = 0; i < MAX_PMTS; i++) {
        if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_NORMAL) continue;

        tac[i/512][nhit[i/512]++] = ev->pmt_hits[i].tac;
    }

    for (i = 0; i < 19; i++) {
        if (nhit[i] <= MIN_NHIT_CRATE) {
            median_tac[i] = 0.0;
            continue;
        }
        gsl_sort(&tac[i][0],1,nhit[i]);
        median_tac[i] = gsl_stats_median_from_sorted_data(&tac[i][0],1,nhit[i]);
    }

    fargsort(median_tac,19,index);

    return (nhit[index[18]] >= MIN_NHIT_BREAKDOWN && median_tac[index[18]] > median_tac[index[17]] + 500);
}

/* Returns 1 if the average time for the hits in the same card as
 * `flasher_pmt_id` is at least 10 ns earlier than all the hits within 4
 * meters of the potential flasher.
 *
 * This cut is designed to help flag flashers in the is_flasher() function for
 * events in which the flashing channel may be missing from the event, so there
 * is no high charge channel to flag it. */
int is_slot_early(event *ev, int flasher_pmt_id)
{
    int i;
    double t_slot, t_nearby;
    int n_slot, n_nearby;
    double flasher_pos[3];
    double pmt_dir[3];
    double distance;

    COPY(flasher_pos,pmts[flasher_pmt_id].pos);

    n_slot = 0;
    n_nearby = 0;
    t_slot = 0.0;
    t_nearby = 0.0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_NORMAL) continue;

        /* Require good calibrations.
         *
         * Note: For some reason PMT hits in the flasher slot sometimes are
         * marked as having bad calibrations and get a time of -9999. I still
         * don't understand why this is, so for now we don't check the
         * calibration bits and instead just check to see that the ECA time is
         * reasonable.
         *
         * Here is what Javi said via email:
         *
         * If SNO calibrations are the same than SNO+ (and I'm pretty sure they
         * are, at last for this matter) there is a hit-level time validation,
         * i.e., if the TAC value is determined to be in the curl region of the
         * time slope (high TAC = very early hits), they get flagged and a
         * calibration is not provided (which is compatible with a flashing
         * channel). If this is the case, you can fall back to raw TACs and QHS
         * when calibration is not available.
         *
         * Another update:
         *
         * I looked into this further and there are cases when the TAC is too
         * high, but in these cases *both* the ECA time and the ECA+PCA time
         * are set to -9999.0. However, there are still cases where *only* the
         * ECA+PCA time is set to -9999.0 but the ECA time is OK. I sent
         * another email to Javi and he figured it out:
         *
         * I just recall that in SNO+ we have a similar hit-level test for PCA:
         * if a channel has a very low QHS, having a good time-walk calibration
         * is hard so we just set the time to invalid. There might have been
         * something similar in SNO, although I didn't find anything after a
         * quick look in the SNO docs.
         *
         * I looked into it and *all* of the hits I saw were indeed caused by a
         * very low QHS value.
         *
         * Here, we use the pt1 value which is the ECA+PCA time except without
         * walk correction which is what was causing the issue with the bad
         * regular time.
         *
         * Another update:
         *
         * When testing out the data cleaning cuts I found an event (run 20664
         * GTID 573780) where an obvious flasher was missed because all the
         * channels in the slot had pt1 = -9999. So now we just use ept which
         * should be good enough. */
        if (ev->pmt_hits[i].ept < -100) continue;

        if (i/32 == flasher_pmt_id/32) {
            /* This hit is in the same slot as the potential flasher. */
            t_slot += ev->pmt_hits[i].ept;
            n_slot += 1;
            continue;
        }

        /* Calculate the distance from the current channel to the high charge
         * channel. */
        SUB(pmt_dir,pmts[i].pos,flasher_pos);
        distance = NORM(pmt_dir);

        if (distance < 400.0) {
            t_nearby += ev->pmt_hits[i].ept;
            n_nearby += 1;
            continue;
        }
    }

    /* If there are no PMTs within 4 meters of the slot, it's almost certainly
     * a flasher. */
    if (n_slot > 4 && n_nearby == 0) return 1;

    if (n_slot == 0 || n_nearby == 0) return 0;

    t_slot /= n_slot;
    t_nearby /= n_nearby;

    return t_slot < t_nearby - 10.0;
}

/* Returns 1 if the event is a flasher and 0 if it isn't. The definition of
 * a flasher event comes from
 * http://detector.sno.laurentian.ca/~detector/private/snoop/doc/eventIDs.html,
 * with a few small changes.
 *
 * Here we define a flasher as an event with the following criteria:
 *
 * - 31 <= nhit <= 1000
 *
 * - Look for the paddle card with the most PMTs hit. This paddle card should
 *   have 4 or more hits.
 *
 * - If such a paddle card exists, look for one channel in the card that hosts
 *   this paddle card which has either QHS or QHL values (uncalib'd) > 4000,
 *   whereas the QHS and QHL values for the other channels in the card are <
 *   800. Note that QHS and QHL values below 300 are sent above 4096.
 *
 *   Note: This is changed from the original definition. We look for all the
 *   other channels to have a QHL and QHS below 800 instead of 700 since the
 *   latter was failing to tag many obvious flashers.
 *
 * - At least 70% of the (regular) pmts in the event which are not in the
 *   flasher slot should have fired more than 50 ns after the high charge
 *   channel.
 *
 * - At least 70% of the (regular) pmts in the event which are not in the
 *   flasher slot should be at a distance >= 12.0 m from the high charge pmt.
 *
 * The one criteria we *don't* use that is on the webpage is the AMB cut:
 *
 * - AMB Int >= 200
 *
 * since I can't seem to figure out where that is stored in the data structure.
 */
int is_flasher(event *ev)
{
    int hits_pc[1280] = {0};
    int qhs, qhl, qlx;
    int qhl_pc[8] = {0};
    int qhs_pc[8] = {0};
    int qlx_pc[8] = {0};
    double t_pc[8] = {0};
    int channel_pc[8] = {0};
    size_t crate, card, channel, id;
    int i, j;
    size_t index[1280], index_qhs[8], index_qhl[8], index_qlx[8];
    double t;
    int flasher_pc;
    double flasher_pos[3];
    double pmt_dir[3];
    double distance;
    int nhit, nhit_late, nhit_far;

    /* Flasher event must have an nhit between 31 and 1000. */
    if (ev->nhit < 31 || ev->nhit > 1000) return 0;

    for (i = 0; i < MAX_PMTS; i++) {
        if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_NORMAL) continue;

	crate = i/512;
	card = (i % 512)/32;
	channel = i % 32;

	id = crate*64 + card*4 + channel/8;

	hits_pc[id] += 1;
    }

    argsort(hits_pc, LEN(hits_pc), index);

    /* The paddle card with the most hits must have >= 4 hits. */
    if (hits_pc[index[1279]] < 4) return 0;

    /* Loop over the most hit paddle cards in order from most hits to least
     * hits. We do this because it's possible that more than one paddle card
     * has the same number of channels hit, and so we need to test all of them. */
    for (j = 1279; j >= 0; j--) {
	id = index[j];

        /* Make sure this paddle card has at least 4 hits. */
        if (hits_pc[id] < 4) return 0;

        nhit = 0;
	for (i = 0; i < MAX_PMTS; i++) {
            if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_NORMAL) continue;

            crate = i/512;
            card = (i % 512)/32;
            channel = i % 32;

	    id = crate*64 + card*4 + channel/8;

	    if (id != index[j]) continue;

            /* Get the uncalibrated QHS and QHL charges. */
	    qhs = ev->pmt_hits[i].qihs;
	    qhl = ev->pmt_hits[i].qihl;
	    qlx = ev->pmt_hits[i].qilx;

            if (qhs < 300) {
                /* QHS values below 300 are set to 4095. */
                qhs_pc[nhit] = 4095;
            } else {
                qhs_pc[nhit] = qhs;
            }

            if (qhl < 300) {
                /* QHL values below 300 are set to 4095. */
                qhl_pc[nhit] = 4095;
            } else {
                qhl_pc[nhit] = qhl;
            }

            if (qlx < 300) {
                /* QLX values below 300 are set to 4095. */
                qlx_pc[nhit] = 4095;
            } else {
                qlx_pc[nhit] = qlx;
            }

            if (ev->pmt_hits[i].ept < -100) {
                /* This can happen if the channel has bad calibration or is in
                 * the TAC curl region.
                 *
                 * It's not really obvious what to do in this situation. To be
                 * conservative, we will just set the time to a negative value.
                 * This means that if this channel ends up being the highest
                 * charge channel, it will basically be guaranteed to be
                 * earlier than the rest of the PMTs, and therefore this part
                 * of the next check will always pass. However it will still
                 * need to pass the cut that 70% of the other PMTs need to be
                 * at least 12 meters away. */
                t_pc[nhit] = -100;
            } else {
                t_pc[nhit] = ev->pmt_hits[i].ept;
            }

            channel_pc[nhit] = channel;

            nhit += 1;
	}

        if (nhit < 2) {
            fprintf(stderr, "only 2 hits on %zu/%zu/PC%zu but expected at least 4!\n",
                    id/64, (id % 64)/4, id % 4);
            continue;
        }

        /* Sort the QHS and QHL values by size. */
	argsort(qhs_pc, nhit, index_qhs);
	argsort(qhl_pc, nhit, index_qhl);
	argsort(qlx_pc, nhit, index_qlx);

        flasher_pc = index[j];
        crate = index[j]/64;
        card = (index[j] % 64)/4;
        channel = channel_pc[index_qhs[nhit-1]];

        /* Check if this paddle card has QHS, QHL, or QLX values > 4000, and the QHS,
         * QHL, and QLX values for the other channels in the card are less than 700.
         *
         * Alternatively, if there isn't an obvious charge outlier, it's still
         * possible that the flasher PMT hit is missing from the event or that
         * the QHS/QHL/QLX values were just outside of the threshold. In this
         * case, we look to see if the TAC for the proposed slot is early with
         * respect to neighboring PMTs. */
        if ((qhs_pc[index_qhs[nhit-1]] < 4000 || qhs_pc[index_qhs[nhit-2]] > 800) &&
	    (qhl_pc[index_qhl[nhit-1]] < 4000 || qhl_pc[index_qhl[nhit-2]] > 800) &&
	    (qlx_pc[index_qlx[nhit-1]] < 1000 || qlx_pc[index_qlx[nhit-2]] > 800) &&
            !is_slot_early(ev, crate*512 + card*32 + channel)) continue;

        if (qhs_pc[index_qhs[nhit-1]] > 4000) {
            channel = channel_pc[index_qhs[nhit-1]];
        } else if (qhl_pc[index_qhl[nhit-1]] > 4000) {
            channel = channel_pc[index_qhl[nhit-1]];
        } else if (qlx_pc[index_qlx[nhit-1]] > 1000) {
            channel = channel_pc[index_qlx[nhit-1]];
        }
        COPY(flasher_pos,pmts[crate*512 + card*32 + channel].pos);

        /* Calculate the median time for PMT hits in the paddle card.
         *
         * Note: Initially this algorithm just used the time of the highest
         * charge channel. However, in looking at events in run 10,000 I
         * noticed two flashers which were not flagged by the cut because the
         * flasher channel had a normal time. Therefore, we now use the median
         * time for all channels in the paddle card. */
        gsl_sort(t_pc,1,nhit);
        t = gsl_stats_median_from_sorted_data(t_pc,1,nhit);

        /* Check that 70% of the regular PMTs in the event which are not in the
         * flasher slot fired more than 50 ns after the high charge channel and are
         * more than 12 meters away from the high charge channel. */
        nhit = 0;
        /* Number of PMTs which fired 50 ns after the high charge channel. */
        nhit_late = 0;
        /* Number of PMTs which are 12 meters or more away from the high charge
         * channel. */
        nhit_far = 0;
        for (i = 0; i < MAX_PMTS; i++) {
            /* Skip PMTs which weren't hit. */
            if (!ev->pmt_hits[i].hit) continue;

            /* Only count regular PMTs without any flags. */
            if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;

            crate = i/512;
            card = (i % 512)/32;
            channel = i % 32;

            id = crate*64 + card*4 + channel/8;

            /* Skip PMTs in the same card as the high charge channel. */
            if (id/4 == flasher_pc/4) continue;

            /* Require good calibrations. */
            if (ev->pmt_hits[i].ept < -100) continue;

            nhit += 1;

            /* Calculate the distance from the current channel to the high charge
             * channel. */
            SUB(pmt_dir,pmts[i].pos,flasher_pos);
            distance = NORM(pmt_dir);

            /* If this channel is 12 meters or more away, increment nhit_far. */
            if (distance >= 1200.0) nhit_far += 1;

            /* If this channel fired more than 50 ns after the high charge channel,
             * increment nhit_late. */
            if (ev->pmt_hits[i].ept > t + 50.0) nhit_late += 1;
        }

        /* If at least 70% of the regular PMTs fired within 50 ns and were at
         * least 12 meters away from the high charge channel, then it's a
         * flasher! */
        if (nhit_late >= 7*nhit/10 && nhit_far >= 7*nhit/10) return 1;
    }

    return 0;
}