#include <RAT/CAENCut.hh>
#include <RAT/DB.hh>
#include <RAT/BitManip.hh>
#include <RAT/Processor.hh>
#include <RAT/CAENBits.hh>
#include <vector>
using namespace std;

namespace RAT {

  void CAENCut::BeginOfRun(DS::Run&)
  {
    char buffer[128];
    const char* sigmoid = "(1/(1+TMath::Exp(-1*(x-[1])/[0])))";
    snprintf(buffer,128,"((1-%s)*[2]+%s*([3]+[4]*x))",sigmoid, sigmoid);

    DBLinkPtr dc_db = DB::Get()->GetLink("DATACLEANING",fName);
    DBLinkPtr daq_rl_db = DB::Get()->GetLink("DAQ_RUN_LEVEL");
    std::vector<int> attenuation_channels = daq_rl_db->GetIArray("attenuationChannels");
    bool attenuated = true;
    for (size_t i=0; i < attenuation_channels.size(); i++) {
        if(attenuation_channels[i] == 0) {
            attenuated = false;
            break;
        }
    }

    fMinNhit = dc_db->GetI("min_nhit");

    fStartSample = dc_db->GetI("start_sample");
    fEndSample = dc_db->GetI("end_sample");

    fPrePulseNumPedBins = dc_db->GetD("pre_pulse_ped_samples");
    fPostPulseNumPedBins = dc_db->GetD("post_pulse_ped_samples");

    double transition_nhit = dc_db->GetI("transition_nhit");
    double transition_width = dc_db->GetI("transition_width");;

    fPeakBounds[0] = TF1("peak_lower",buffer,0,10000);
    fPeakBounds[1] = TF1("peak_upper",buffer,0,10000);
    fIntegralBounds[0] = TF1("integral_lower",buffer,0,10000);
    fIntegralBounds[1] = TF1("integral_upper",buffer,0,10000);

    double peak_low_nhit_constant[2];
    double peak_linear_constant[2];
    double peak_linear_slope[2];

    double integral_low_nhit_constant[2];
    double integral_linear_constant[2];
    double integral_linear_slope[2];

    string prefix = "";

    if(attenuated) {
        prefix = "attenuated_";
    }
    peak_low_nhit_constant[0] = dc_db->GetD(prefix + "peak_low_nhit_lower_bound");
    peak_low_nhit_constant[1] = dc_db->GetD(prefix + "peak_low_nhit_upper_bound");
    peak_linear_constant[0]   = dc_db->GetD(prefix + "peak_lower_constant");
    peak_linear_constant[1]   = dc_db->GetD(prefix + "peak_upper_constant");
    peak_linear_slope[0]      = dc_db->GetD(prefix + "peak_lower_slope");
    peak_linear_slope[1]      = dc_db->GetD(prefix + "peak_upper_slope");


    integral_low_nhit_constant[0] = dc_db->GetD(prefix + "integral_low_nhit_lower_bound");
    integral_low_nhit_constant[1] = dc_db->GetD(prefix + "integral_low_nhit_upper_bound");
    integral_linear_constant[0]   = dc_db->GetD(prefix + "integral_lower_constant");
    integral_linear_constant[1]   = dc_db->GetD(prefix + "integral_upper_constant");
    integral_linear_slope[0]      = dc_db->GetD(prefix + "integral_lower_slope");
    integral_linear_slope[1]      = dc_db->GetD(prefix + "integral_upper_slope");

    for(int i=0;i<2;i++) {
        fPeakBounds[i].SetParameter(0,transition_width);
        fPeakBounds[i].SetParameter(1,transition_nhit);
        fPeakBounds[i].SetParameter(2,peak_low_nhit_constant[i]);
        fPeakBounds[i].SetParameter(3,peak_linear_constant[i]);
        fPeakBounds[i].SetParameter(4,peak_linear_slope[i]);
        fIntegralBounds[i].SetParameter(0,transition_width);
        fIntegralBounds[i].SetParameter(1,transition_nhit);
        fIntegralBounds[i].SetParameter(2,integral_low_nhit_constant[i]);
        fIntegralBounds[i].SetParameter(3,integral_linear_constant[i]);
        fIntegralBounds[i].SetParameter(4,integral_linear_slope[i]);
    }


    fCutMode = dc_db->GetI("cut_mode");
    fUseLowerCut = (bool)dc_db->GetI("use_lower_bound");

    fLDAQ = DB::Get()->GetLink("DAQ_RUN_LEVEL");
    fEsumTrig = ESHiLo;//This cut (currently) only works with ESUMH-Lo
    fHasWarned = false;
  }

  Processor::Result CAENCut::DSEvent(DS::Run&, DS::Entry& ds)
  {
    bool pass = true; // If any triggered event fails, they all fail
    for( size_t iEV = 0; iEV < ds.GetEVCount(); iEV++ )
      if( Event(ds, ds.GetEV(iEV)) != OKTRUE )
        pass = false;
    return pass ? OKTRUE : OKFALSE;
  }

  Processor::Result CAENCut::Event(DS::Entry&, DS::EV& ev)
  {
    double integral(0), min(0), pedValue(0);
    const int nhit = ev.GetNhits();
    const int firstInitialPedBin = fStartSample;
    const int lastInitialPedBin = fStartSample + fPrePulseNumPedBins;
    const int firstEndingPedBin = fEndSample - fPostPulseNumPedBins;
    const int lastEndingPedBin = fEndSample;
    const int numSamples = fEndSample-(fStartSample+1);
    const int numPedBins = fPrePulseNumPedBins + fPostPulseNumPedBins;
    const int numPulseBins = numSamples - numPedBins;
    bool hasEsum = true;
    fPassFlag = true;
    if( numPedBins <= 0 && !fHasWarned)
    {
            warn<<"Running CAEN Cut without pedestal subtraction"<<newline;
            fHasWarned = true;
    }


    try { hasEsum = ev.GetDigitiser().ExistsWaveform(fEsumTrig); }
    catch( RAT::DS::DataNotFound& e ) { hasEsum = false; }
    // check if nhit is high enough to attempt cut
    // and that we are actually digitizing esumhi
    if (nhit < fMinNhit || fEsumTrig < 0 || !hasEsum){
      UpdateMask(ev);
      return fPassFlag ? OKFALSE : OKTRUE;
    }

    // calculate integral and peak size
    DS::Digitiser& dig = ev.GetDigitiser();

    try {
        integral = dig.Integral(fEsumTrig, lastInitialPedBin, firstEndingPedBin);
        min = dig.Min(fEsumTrig,lastInitialPedBin,firstEndingPedBin);
        if(numPedBins > 0) {
            double start_ped = dig.Average(fEsumTrig, firstInitialPedBin,
                                                      lastInitialPedBin);
            start_ped *= fPrePulseNumPedBins/(double)(numPedBins);
            double end_ped = dig.Average(fEsumTrig, firstEndingPedBin,
                                                    lastEndingPedBin);
            end_ped *= fPostPulseNumPedBins/(double)(numPedBins);
            pedValue = start_ped+end_ped;
        }
    } catch( std::out_of_range& oor) {
        return FAIL;
    }


    //subtract the pedestal from the peak
    double peak = pedValue - min;

    // subtract pedestal from integral
    integral = integral-pedValue*numPulseBins;

    // check if each passes
    int peakFail = 0;
    int integralFail = 0;
    if( ((fPeakBounds[0].Eval(nhit) > peak) && fUseLowerCut) ||
         (fPeakBounds[1].Eval(nhit) < peak))
    {
        peakFail = 1;
    }
    if( ((fIntegralBounds[0].Eval(nhit) > integral) && fUseLowerCut) ||
         (fIntegralBounds[1].Eval(nhit) < integral) )
    {
        integralFail = 1;
    }

    // determine if CAENCut is passed or not
    if (fCutMode == 1){
      if (peakFail || integralFail)
        fPassFlag = false;
    }else{
      if (peakFail && integralFail)
        fPassFlag = false;
    }

    // now update the mask
    UpdateMask(ev);
    return fPassFlag ? OKFALSE : OKTRUE;
  }

} // namespace RAT