#include <RAT/ITCTimeSpreadCut.hh>
#include <RAT/DB.hh>
#include <RAT/BitManip.hh>
#include <RAT/Processor.hh>
#include <algorithm>

using namespace std;

namespace RAT {

  void ITCTimeSpreadCut::BeginOfRun(DS::Run&)
  {
    fLClean = DB::Get()->GetLink("DATACLEANING",fName);
    fWindowSize = fLClean->GetD("window_size");
    fWindowStep = fLClean->GetD("window_step");
    fRatioCut = fLClean->GetD("ratio_cut");
  }

  Processor::Result ITCTimeSpreadCut::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 ITCTimeSpreadCut::Event(DS::Entry&, DS::EV& ev)
  {
    vector<double> pmt_times;
    fPassFlag = true;
    double min_time,max_time;
    int max_count = -1;

    for (size_t ipmt=0;ipmt<ev.GetCalPMTs().GetCount();ipmt++){
      const DS::PMTCal& pmt = ev.GetCalPMTs().GetPMT(ipmt);
      // if pmt time is actually calibrated
      if (abs(pmt.GetTime()) < 8000)
        pmt_times.push_back(pmt.GetTime());
    }

    // also skip if no tubes with actual time calibrations
    if (pmt_times.size() == 0){
      UpdateMask(ev);
      return fPassFlag ? OKFALSE : OKTRUE;
    }

    sort(pmt_times.begin(),pmt_times.end());
    int npmt = pmt_times.size()-1;

    min_time = pmt_times[0];
    max_time = pmt_times[npmt];

    // if all the hits are within one window
    // we are done
    if ((max_time - min_time) <= fWindowSize){
      max_count = npmt+1;
    }else{

      double itime = min_time;
      int i_win_low,i_win_high;
      while (itime <= (max_time - fWindowSize + fWindowStep)){
        // itime is the time in ns for the start of our new
        // window. We have to find which pmt this corresponds to.
        // We will do a binary search for speed
        int ilow = 0;
        int ihigh = npmt;
        while ((ihigh - ilow) > 1){
          int imid = (ilow + ihigh) / 2;
          if (itime == pmt_times[imid]){
            ilow = imid-1;
            ihigh = ilow;
          }else if (itime > pmt_times[imid]){
            ilow = imid;
          }else{
            ihigh = imid;
          }
        }
        i_win_low = ilow;

        // again binary search for top position of window
        ilow = 0;
        ihigh = npmt;
        while((ihigh-ilow) > 1){
          int imid = (ilow + ihigh) / 2;
          if (((itime + fWindowSize) >= pmt_times[imid]) &&
              ((itime + fWindowSize) <= (pmt_times[imid] + fWindowStep))){
            ilow = imid-1;
            ihigh = ilow;
          }else if ((itime + fWindowSize) < pmt_times[imid]){
            ihigh = imid;
          }else{
            ilow = imid;
          }
        }
        i_win_high = ilow;

        if ((i_win_high - i_win_low) > max_count)
          max_count = (i_win_high-i_win_low);

        itime += fWindowStep;
      } // end while itime <= max_time-win_size+win_step
    } // end if/else hits are all in one window

    double max_ratio = (double) max_count / (double) pmt_times.size();
    if (max_ratio < fRatioCut)
      fPassFlag = false;

    UpdateMask(ev);
    return fPassFlag ? OKFALSE : OKTRUE;
  }

} // namespace RAT