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"

/* 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.
 *
 * Note: This function may return -1 if there is an error traversing the PMT
 * banks. */
int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev)
{
    size_t 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 (!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 1 if the event is classified as a neck event. The definition of
 * a neck event comes from the SNOMAN companion:
 *
 *     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. 
 *
 */
int is_neck_event(event *ev)
{
    size_t i;
    int n;
    int high_charge;
    double avg_ept_psup;
    int nhit;

    /* First, we compute the average ECA time for all normal PMTs. */

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

        if (ev->pmt_hits[i].hit) {
            avg_ept_psup += ev->pmt_hits[i].ept;
            nhit += 1;
        }
    }

    if (nhit < 1) return 0;

    avg_ept_psup /= nhit;

    /* 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;
    for (i = 0; i < MAX_PMTS; i++) {
        if (pmts[i].pmt_type == PMT_NECK && ev->pmt_hits[i].hit) {
            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 - 70.0)) high_charge = 1;
        }
    }

    if (n >= 2 || high_charge) return 1;

    return 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;
    int qhl_pc[8] = {0};
    int qhs_pc[8] = {0};
    double t_pc[8] = {0};
    size_t crate, card, channel, id;
    int i, j;
    size_t index[1280], index_qhs[8], index_qhl[8];
    double t;
    int flasher_pc;
    double flasher_pos[3];
    double pmt_dir[3];
    int flasher = 0;
    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];

        /* If this paddle card doesn't have as many hits as the most hit paddle
         * card and we haven't found a flasher yet, then it's not a flasher. */
        if (hits_pc[id] < hits_pc[index[1279]]) return 0;

	for (i = 0; i < MAX_PMTS; i++) {
            if (!ev->pmt_hits[i].hit) 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;

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

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

            t_pc[channel % 8] = ev->pmt_hits[i].t;
	}

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

        /* Check if this paddle card has QHS or QHL values > 4000, and the QHS
         * and QHL values for the other channels in the card are less than 700. */
        if (qhs_pc[index_qhs[7]] > 4000 || qhl_pc[index_qhl[7]] > 4000) {
	    if (qhs_pc[index_qhs[6]] < 800 && qhl_pc[index_qhl[6]] < 800) {
                /* Flasher! */
                flasher = 1;
                break;
	    }
	}
    }

    /* If we didn't find a flasher, return 0. */
    if (!flasher) return 0;

    flasher_pc = index[j];
    crate = index[j]/64;
    card = (index[j] % 64)/4;
    if (qhs_pc[index_qhs[7]] > 4000) {
        channel = (index[j] % 4)*8 + index_qhs[7];
        t = t_pc[index_qhs[7]];
    } else {
        channel = (index[j] % 4)*8 + index_qhl[7];
        t = t_pc[index_qhl[7]];
    }
    COPY(flasher_pos,pmts[crate*512 + card*32 + channel].pos);

    /* 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 paddle card as the high charge channel. */
        if (id == flasher_pc) 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].t > t + 50.0) nhit_late += 1;
    }

    /* If less than 70% of the regular PMTs fired within 50 ns of the high
     * charge channel, then it's not a flasher. */
    if (nhit_late < nhit*0.7) return 0;

    /* If less than 70% of the regular PMTs were closer than 12 meters away
     * from the high charge channel, then it's not a flasher. */
    if (nhit_far < nhit*0.7) return 0;

    return 1;
}