#include <RAT/QvNhit.hh>
#include <RAT/DB.hh>

using namespace std;

namespace RAT {

  void QvNhit::BeginOfRun(DS::Run&)
  {
    fLClean = DB::Get()->GetLink("DATACLEANING",fName);
    fLPmtCalib = DB::Get()->GetLink("PMTCALIB");
    fQHL_HHP = fLPmtCalib->GetDArray("QHL_hhp");
    fRatio = fLClean->GetD("ratio");
    fECAhitmask = fLClean->GetI("ECA_mask");
    fPCAhitmask = fLClean->GetI("PCA_mask");
    fHHPlow = fLClean->GetD("HHP_low");
    fHHPhigh = fLClean->GetD("HHP_high");
  }

  Processor::Result QvNhit::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 QvNhit::Event(DS::Entry&, DS::EV& ev)
  {

    int id,nhit_cal,ncut,qlarge_index;
    float qhl,qhl_large,qhl_tot,qhl_ratio;
    std::vector<float> qhlarray;
    std::vector<int> qlargearray;
    std::vector<int> idarray;

    fPassFlag = true;

    // get the number of normal pmts hit (that calibrations have been attempted on)
    size_t calcount = ev.GetCalPMTs().GetCount();
    // initialization
    qhlarray.clear();
    qlargearray.clear();
    idarray.clear();
    nhit_cal = 0;

    for (size_t ipmt=0;ipmt<calcount;ipmt++){
      const DS::PMTCal& pmt = ev.GetCalPMTs().GetPMT(ipmt);
      id = pmt.GetID();
      // we use qhl because pickup will be ac noise and will integrate to 0
      qhl = pmt.GetQHL();
      if (abs(qhl) < 8000){ // if we have a reasonable QHL
        // we build up an array of all the physical charges
        if (qhl > -100){ // if charge is physical
          int ECA_status = (int)pmt.GetStatus().GetBits(0,32);
          int PCA_status = (int)pmt.GetStatus().GetBits(32,32);
          if (ECA_status & fECAhitmask)
            continue;
          if (PCA_status & fPCAhitmask)
            continue;
          //Check for a sane HHP value
          if ((fHHPlow < fQHL_HHP[id]) && (fQHL_HHP[id] < fHHPhigh)){
            nhit_cal++;
            qlargearray.push_back(0);
            idarray.push_back(id);
            qhlarray.push_back(qhl);
          }
        }
      }
    }

    // now we find the 10% of hits with the highest charge
    ncut = static_cast<int>(floor( nhit_cal*0.1 ));
    if ((nhit_cal-ncut)<5){
      // there are not enough calibrated good tubes to see if this event should be cut
      // so it passes
      UpdateMask(ev);
      return fPassFlag ? OKFALSE : OKTRUE;
    }

    for (int i=0;i<ncut;i++){
      qhl_large = -9999;
      qlarge_index = 0;
      for (int j=0;j<nhit_cal;j++){
        if (qhlarray[j] > qhl_large && qlargearray[j] == 0){
          qhl_large = qhlarray[j];
          qlarge_index = j;
        }
      }
      qlargearray[qlarge_index] = 1;
    }

    // now that we have flagged the highest charge 10%, we sum up the total charge
    // from the remaining lower 90%
    qhl_tot = 0;
    for (int i=0;i<nhit_cal;i++){
      if (qlargearray[i] == 0){
        // we want the charge in adc counts so we scale by hhp
        qhl_tot += (float) qhlarray[i] / fQHL_HHP[idarray[i]];
      }
    }

    // now calculate charge/nhit ratio, decide if it is pickup
    qhl_ratio = qhl_tot / (float) (nhit_cal-ncut);

    // now update the mask
    if (qhl_ratio < fRatio)
      fPassFlag = false;

    UpdateMask(ev);

    return fPassFlag ? OKFALSE : OKTRUE;
  }

} // namespace RAT