/* 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 "event.h"
#include "sort.h"
#include <stddef.h> /* for size_t */
#include "misc.h"
#include "pmt.h"
#include "vector.h"

/* 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;
}