#include <RAT/TriggerEfficiencyProc.hh>
#include <RAT/Log.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DU/ChanHWStatus.hh>

#include <cmath>

using namespace RAT;

void low_pass_filter(float *yin, float *yout, int n, float dt, float rise_time)
{
    /* Simulate a simple RC low pass filter. See
     * https://wikipedia.org/wiki/Low-pass_filter#Simple_infinite_impulse_response_filter. */
    int i;
    float rc = rise_time/2.197;
    float alpha = dt/(rc + dt);

    yout[0] = alpha*yin[0];

    for (i = 1; i < n; i++) {
        yout[i] = alpha*yin[i] + (1-alpha)*yout[i-1];
    }
}

void reconstruct_trigger(DS::EV *ev, unsigned int n, float dt, float *n100, float *n20, float n100_width, float n20_width)
{
    /* Reconstructs an idealized trigger signal from the calibrated PMT hit
     * times of each channel with triggers enabled.
     *
     * Note: The reconstructed trigger signal is positive going and the scale is
     * in units of nhit. */
    unsigned int i, j;
    int crate, card, id;
    float time;
    float t;
    bool nhit_100_enabled, nhit_20_enabled;

    for (i = 0; i < n; i++) {
        n100[i] = 0;
        n20[i] = 0;
    }

    const RAT::DU::ChanHWStatus &channelHardwareStatus = DU::Utility::Get()->GetChanHWStatus();

    const RAT::DS::CalPMTs &calPMTs = ev->GetCalPMTs();

    for (j = 0; j < calPMTs.GetAllCount(); j++) {
        const RAT::DS::PMTCal &pmt = calPMTs.GetAllPMT(j);
        crate = pmt.GetCrate();
        card = pmt.GetCard();
        time = pmt.GetTime();

        if ((crate == 3 || crate == 13 || crate == 17 || crate == 18) && card == 15) {
            /* Skip the FEC/D and OWL slots since the triggers aren't added
             * to the normal N100/N20 trigger signals. */
            continue;
        }

        id = pmt.GetLCN();

        nhit_100_enabled = channelHardwareStatus.IsNHit100Enabled(id);
        nhit_20_enabled = channelHardwareStatus.IsNHit20Enabled(id);

        if (!nhit_100_enabled && !nhit_20_enabled) continue;

        for (i = 0; i < n; i++) {
            t = i*dt;
            if (nhit_100_enabled && (t > time) && (t < (time + n100_width))) n100[i] += 1.0;
            if (nhit_20_enabled && (t > time) && (t < (time + n20_width))) n20[i] += 1.0;
        }
    }
}

void TriggerEfficiencyProc::BeginOfRun(DS::Run&)
{
    DBLinkPtr fDB = DB::Get()->GetLink("TRIGGER_EFFICIENCY");

    fN100Width = fDB->GetD("n100_width");
    fN20Width = fDB->GetD("n20_width");
    fN100RiseTime = fDB->GetD("n100_rise_time");
    fN20RiseTime = fDB->GetD("n20_rise_time");
    fTimeWindow = fDB->GetD("time_window");
    fDeltaT = fDB->GetD("delta_t");
}

void TriggerEfficiencyProc::GetInTimeHits(DS::EV *ev, float *max_nhit_100, float *max_nhit_20)
{
    /* Calculate the maximum number of channels which contributed to the
     * trigger.
     *
     * This function reconstructs an idealized trigger signal and then finds the
     * peak of the trigger signal. Since it includes a low pass filter, the
     * maximum nhit returned is a floating point number. */
    unsigned int i;
    unsigned int trigger_signal_length = (unsigned int) ceil(fTimeWindow/fDeltaT);
    float *n100 = new float[trigger_signal_length];
    float *n20 = new float[trigger_signal_length];
    float *n100_low_pass = new float[trigger_signal_length];
    float *n20_low_pass = new float[trigger_signal_length];

    reconstruct_trigger(ev, trigger_signal_length, fDeltaT, n100, n20, fN100Width, fN20Width);

    low_pass_filter(n100, n100_low_pass, trigger_signal_length, fDeltaT, fN100RiseTime);
    low_pass_filter(n20, n20_low_pass, trigger_signal_length, fDeltaT, fN20RiseTime);

    *max_nhit_100 = 0.0;
    *max_nhit_20 = 0.0;
    for (i = 0; i < trigger_signal_length; i++) {
        if (n100_low_pass[i] > *max_nhit_100) *max_nhit_100 = n100_low_pass[i];
        if (n20_low_pass[i] > *max_nhit_20) *max_nhit_20 = n20_low_pass[i];
    }

    delete[] n100;
    delete[] n20;
    delete[] n100_low_pass;
    delete[] n20_low_pass;
}

Processor::Result TriggerEfficiencyProc::DSEvent(DS::Run&, DS::Entry& ds)
{
    unsigned int i;
    float max_nhit_100, max_nhit_20;

    for (i = 0; i < ds.GetEVCount(); i++) {
        DS::EV &ev = ds.GetEV(i);

        GetInTimeHits(&ev, &max_nhit_100, &max_nhit_20);

        ev.SetInTimeHits100(max_nhit_100);
        ev.SetInTimeHits20(max_nhit_20);
    }

    return Processor::OKTRUE;
}