#include <RAT/QvTCut.hh>
#include <RAT/DB.hh>
#include <RAT/BitManip.hh>
#include <RAT/ECACal.hh>

using namespace std;

namespace RAT {

  void QvTCut::BeginOfRun(DS::Run& run)
  {
    fLClean = DB::Get()->GetLink("DATACLEANING",fName);
    fTimeMode = fLClean->GetI("timemode");
    fCutMode = fLClean->GetI("cutmode");
    fEarlyCTimeBound = fLClean->GetD("earlycaltimebound");
    fLateCTimeBound = fLClean->GetD("latecaltimebound");
    fEarlyUTimeBound = fLClean->GetD("earlyuncaltimebound");
    fLateUTimeBound = fLClean->GetD("lateuncaltimebound");
    fQHBound = fLClean->GetI("qhbound");
    fQLBound = fLClean->GetI("qlbound");
    fQHSECAmask = fLClean->GetI("QHSECAmask");
    fQHLECAmask = fLClean->GetI("QHLECAmask");
    fQLXECAmask = fLClean->GetI("QLXECAmask");
    fPCAmask = fLClean->GetI("PCAmask");
    fBadTACHigh = fLClean->GetI("badTACHigh");
    fBadTACLow = fLClean->GetI("badTACLow");
    fLowTCAL = fLClean->GetI("lowTCAL");
    fMinPMTs = fLClean->GetI("minPMTs");
    fNegativeRail = fLClean->GetI("negativeRail");
    fTCurlRaw = fLClean->GetI("tCurlRaw");
    fECACal.BeginOfRun(run);
  }

  Processor::Result QvTCut::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 QvTCut::Event(DS::Entry&, DS::EV& ev)
  {
    vector<unsigned short> utime;
    int utacCount = 0;
    vector<double> etime;
    int etacCount = 0;
    vector<double> ctime;
    int tacCount = 0;

    float qhsTot = 0;
    int qhsCount = 0;
    float qhsMax = -9999;
    double qhsMax_tac = 0;
    double qhsMax_etac = 0;
    unsigned short qhsMax_utac = 0;

    float qhlTot = 0;
    int qhlCount = 0;
    float qhlMax = -9999;
    double qhlMax_tac = 0;
    double qhlMax_etac = 0;
    unsigned short qhlMax_utac = 0;

    float qlxTot = 0;
    int qlxCount = 0;
    float qlxMax = -9999;
    double qlxMax_tac = 0;
    double qlxMax_etac = 0;
    unsigned short qlxMax_utac = 0;

    float qhsAvg, qhlAvg, qlxAvg;
    double lateBound;
    double earlyBound;
    int qhsFail = 0;
    int qhlFail = 0;
    int qlxFail = 0;
    double medianTac;



    int badECATAC = -9999; // this is how they are flagged in ECA processor right now
    int noECATAC = -9998;
    int rail = 4100; // charge set to above max ADC; makes discernible from nonrail
    fPassFlag = true;

    // get the number of normal pmts hit
    const size_t calcount = ev.GetCalPMTs().GetCount();

    for (size_t ipmt=0;ipmt<calcount;ipmt++){
      const DS::PMTCal& cpmt = ev.GetCalPMTs().GetPMT(ipmt);
      const DS::PMTUncal& upmt = ev.GetUncalPMTs().GetPMT(ipmt);

      double tac;
      double etac;
      unsigned short utac;

      float qlx = cpmt.GetQLX();
      float qhl = cpmt.GetQHL();
      float qhs = cpmt.GetQHS();

      // fill one array with uncalibrated data in case we're missing a calibration
      utac = upmt.GetTime();
      if ((utac < fBadTACHigh) && (utac > fBadTACLow)){
        utime.push_back(utac);
        utacCount++;
      }

      if (fTimeMode > 1){
        // apply only ECA calibrations in case PCA calibration is bad
        fECACal.CalibrateHit(upmt.GetQHS(),upmt.GetQHL(),upmt.GetQLX(),utac,upmt.GetCrateCardChannelCell());
        // if tac is good
        etac = fECACal.GetTime();
        if (etac > fLowTCAL){
          etime.push_back(etac);
          etacCount++;
        }
      }

      // now fill an array with full calibration (at least whatever is in pmtcal)
      tac = cpmt.GetTime();
      if (tac > fLowTCAL){
        ctime.push_back(tac);
        tacCount++;
      }

      //using good Q calibrated tubes, find the highest Q tube
      int ECA_status = (int)cpmt.GetStatus().GetBits(0,32);
      int PCA_status = (int)cpmt.GetStatus().GetBits(32,32);

      if (((ECA_status & fQHSECAmask)==0) && ((PCA_status & fPCAmask)==0)){
        qhsTot += qhs;
        qhsCount++;
        if (qhs < fNegativeRail)
          qhs = rail;
        if ((qhs >= qhsMax) && (utac < fBadTACHigh) && (utac > fBadTACLow)){
          qhsMax = qhs;
          qhsMax_tac = tac;
          qhsMax_utac = utac;
          qhsMax_etac = etac;
        }
      }

      if (((ECA_status & fQHLECAmask)==0) && ((PCA_status & fPCAmask)==0)){
        qhlTot += qhl;
        qhlCount++;
        if (qhl < fNegativeRail)
          qhl = rail;
        if ((qhl >= qhlMax) && (utac < fBadTACHigh) && (utac > fBadTACLow)){
          qhlMax = qhl;
          qhlMax_tac = tac;
          qhlMax_utac = utac;
          qhlMax_etac = etac;
        }
      }

      if (((ECA_status & fQLXECAmask)==0) && ((PCA_status & fPCAmask)==0)){
        qlxTot += qlx;
        qlxCount++;
        if (qlx < fNegativeRail)
          qlx = rail;
        if ((qlx >= qlxMax) && (utac < fBadTACHigh) && (utac > fBadTACLow)){
          qlxMax = qlx;
          qlxMax_tac = tac;
          qlxMax_utac = utac;
          qlxMax_etac = etac;
        }
      }
    } // end loop over pmts

    // For each cut mode, check for the required number of hits
    // and if the max charge tube TAC and eTAC are valid
    bool maxtacs_bad = false;
    bool maxetacs_bad = false;
    bool minQHitsNeeded = false;
    bool minTHitsNeeded = false;

    if ((tacCount > fMinPMTs) || (etacCount > fMinPMTs) ||
         (utacCount > fMinPMTs))
      minTHitsNeeded = true;

    if ((fCutMode == 1) || (fCutMode == 2)){
      if ((qhlCount > fMinPMTs) && (qlxCount > fMinPMTs))
        minQHitsNeeded = true;
      if (qlxMax_tac == badECATAC || qlxMax_tac == noECATAC ||
           qhlMax_tac == badECATAC || qhlMax_tac == noECATAC)
        maxtacs_bad = true;
      if (qlxMax_etac == badECATAC || qlxMax_etac == noECATAC ||
           qhlMax_etac == badECATAC || qhlMax_etac == noECATAC)
        maxetacs_bad = true;
    }else if ((fCutMode == 3) || (fCutMode == 4)){
      if ((qhsCount > fMinPMTs) && (qlxCount > fMinPMTs))
        minQHitsNeeded = true;
      if (qlxMax_tac == badECATAC || qlxMax_tac == noECATAC ||
           qhlMax_tac == badECATAC || qhlMax_tac == noECATAC)
        maxtacs_bad = true;
      if (qlxMax_etac == badECATAC || qlxMax_etac == noECATAC ||
           qhlMax_etac == badECATAC || qhlMax_etac == noECATAC)
        maxetacs_bad = true;
    }else {
      if ((qhlCount > fMinPMTs) && (qlxCount > fMinPMTs) &&
           (qhsCount > fMinPMTs))
        minQHitsNeeded = true;
      if (qlxMax_tac == badECATAC || qlxMax_tac == noECATAC ||
           qhlMax_tac == badECATAC || qhlMax_tac == noECATAC ||
           qhsMax_tac == badECATAC || qhsMax_tac == noECATAC)
        maxtacs_bad = true;
      if (qlxMax_etac == badECATAC || qlxMax_etac == noECATAC ||
           qhlMax_etac == badECATAC || qhlMax_etac == noECATAC ||
           qhsMax_etac == badECATAC || qhsMax_etac == noECATAC)
        maxetacs_bad = true;
    }

    // now only do calculation if there is enough data
    if (minQHitsNeeded && minTHitsNeeded){
      qhsAvg = qhsTot / (float) (qhsCount - 1);
      qhlAvg = qhlTot / (float) (qhlCount - 1);
      qlxAvg = qlxTot / (float) (qlxCount - 1);

      // calculate median time. If PCA calibrated time is bad, try ECA only
      // if that is also bad, default to uncalibrated
      if ((fTimeMode > 1) && maxtacs_bad){
        // check if ECA values are also bad
        if ((fTimeMode > 2) && maxetacs_bad){
          //Use uncalibrated times
          if (utacCount%2){
            nth_element(utime.begin(),utime.begin()+utime.size()/2,utime.end());
            medianTac = utime[utime.size()/2];
          }else{
            nth_element(utime.begin(),utime.begin()+utime.size()/2,utime.end());
            nth_element(utime.begin(),utime.begin()+utime.size()/2-1,utime.end());
            medianTac = (utime[utime.size()/2] + utime[utime.size()/2-1]);
          }
          qhsMax_tac = qhsMax_utac;
          qhlMax_tac = qhlMax_utac;
          qlxMax_tac = qlxMax_utac;
          // sign of cut is reversed for uncalibrated
          lateBound = -fEarlyUTimeBound;
          earlyBound = -fLateUTimeBound;
        }else{
          // use ECA calibrated times
          if (etacCount%2){
            nth_element(etime.begin(),etime.begin()+etime.size()/2,etime.end());
            medianTac = etime[etime.size()/2];
          }else{
            nth_element(etime.begin(),etime.begin()+etime.size()/2,etime.end());
            nth_element(etime.begin(),etime.begin()+etime.size()/2-1,etime.end());
            medianTac = (etime[etime.size()/2] + etime[etime.size()/2-1]);
          }
          qhsMax_tac = qhsMax_etac;
          qhlMax_tac = qhlMax_etac;
          qlxMax_tac = qlxMax_etac;
          earlyBound = fEarlyCTimeBound;
          lateBound = fLateCTimeBound;
        }
      }else{
        // use fully calibrated tac times
        if (tacCount%2){
          nth_element(ctime.begin(),ctime.begin()+ctime.size()/2,ctime.end());
          medianTac = ctime[ctime.size()/2];
        }else if(ctime.size()!=0){
          nth_element(ctime.begin(),ctime.begin()+ctime.size()/2,ctime.end());
          nth_element(ctime.begin(),ctime.begin()+ctime.size()/2-1,ctime.end());
          medianTac = (ctime[ctime.size()/2] + ctime[ctime.size()/2-1]);
        } else {
	//no ctimes
	}
        earlyBound = fEarlyCTimeBound;
        lateBound = fLateCTimeBound;
      } // done deciding which tac times to use

      // if maximum high gain charge tube has a bad calibrated tac time, and the
      // raw tac time is in the tac curl, it is considered to be in the early
      // time box and its time is set to the earliest edge of that box (the
      // late time boundary) if the user has opted for time mode > 0
      if (fTimeMode > 0 && qhsMax_tac == badECATAC && qhsMax_utac > fTCurlRaw)
        qhsMax_tac = medianTac - fLateCTimeBound;
      if (fTimeMode > 0 && qhlMax_tac == badECATAC && qhlMax_utac > fTCurlRaw)
        qhlMax_tac = medianTac - fLateCTimeBound;

      // now apply cuts to qhs, qhl, and qlx
      if ((qhsMax-qhsAvg >= fQHBound) &&
          ((medianTac-qhsMax_tac) >= earlyBound) && ((medianTac-qhsMax_tac) <= lateBound)){
        qhsFail = 1;
      }

      if ((qhlMax-qhlAvg >= fQHBound) &&
          ((medianTac-qhlMax_tac) >= earlyBound) && ((medianTac-qhlMax_tac) <= lateBound)){
        qhlFail = 1;
      }

      if ((qlxMax-qlxAvg >= fQLBound) &&
          ((medianTac-qlxMax_tac) >= earlyBound) && ((medianTac-qlxMax_tac) <= lateBound)) {
        qlxFail = 1;
      }

      if (fCutMode == 1){ // use or of qhl and qlx cuts
        if (qhlFail || qlxFail)
          fPassFlag = false;
      }else if (fCutMode == 2){ // use and of qhl and qlx cuts
        if (qhlFail && qlxFail)
          fPassFlag = false;
      }else if (fCutMode == 3){ // use or of qhs and qlx cuts
        if (qhsFail || qlxFail)
          fPassFlag = false;
      }else if (fCutMode == 4){ // use and of qhs and qlx cuts
        if (qhsFail && qlxFail)
          fPassFlag = false;
      }else if (fCutMode == 5){ // use qhl&qlx || qhs&qlx
        if ((qhsFail && qlxFail) || (qhlFail && qlxFail))
          fPassFlag = false;
      }
    } // done checking if there is enough data/running cut

    UpdateMask(ev);

    return fPassFlag ? OKFALSE : OKTRUE;
  }

} // namespace RAT