#include <RAT/FTSCut.hh>
#include <RAT/DB.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>

using namespace std;

namespace RAT {

  void FTSCut::BeginOfRun(DS::Run&)
  {
    fLClean = DB::Get()->GetLink("DATACLEANING",fName);
    fClusterMode = fLClean->GetI("clustermode");
    fDeltaTThresh = fLClean->GetI("deltatthresh");
    fDistThresh = fLClean->GetI("distthresh");
    fCountThresh = fLClean->GetI("countthresh");
    fMedianCut = fLClean->GetD("mediancut");
    fMaxPairs = 10000;
  }

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

    int cal_count,pair_count,cluster;
    int iccc,jccc;
    float tubedist, deltat, dtmedian;
    vector<DS::PMTCal *> pmt_list;
    vector<int> ccc_hits(19*16*32,0);
    vector<float> deltat_list;

    fPassFlag = true;

    // get a list of the pmts with good
    // time calibrations
    cal_count = ev.GetCalPMTs().GetCount();
    for (int ipmt=0;ipmt<cal_count;ipmt++){
      DS::PMTCal* pmt = &ev.GetCalPMTs().GetPMT(ipmt);
      // if time is calibrated
      // add to the list
      if (abs(pmt->GetTime()) < 8000){
        pmt_list.push_back(pmt);
        ccc_hits[pmt->GetID()] = 1;
      }
    }

    pair_count = 0;
    cluster = 0;
    // loop over all pairs of hit pmts
    for (size_t ipmt=0;ipmt<pmt_list.size();ipmt++){
      if (pair_count >= fMaxPairs)
        break;
      for (size_t jpmt=ipmt;jpmt<pmt_list.size();jpmt++){
        if (pair_count >= fMaxPairs)
          break;
        if (ipmt == jpmt)
          continue;

        iccc = pmt_list[ipmt]->GetID();
        jccc = pmt_list[jpmt]->GetID();
        TVector3 ipos = DU::Utility::Get()->GetPMTInfo().GetPosition(iccc);
        TVector3 jpos = DU::Utility::Get()->GetPMTInfo().GetPosition(jccc);

        // compute the linear distance between the
        // pmts and the absolute value of their time
        // difference
        tubedist = sqrt(pow((ipos.X() - jpos.X()),2) + pow((ipos.Y() - jpos.Y()),2) + pow((ipos.Z() - jpos.Z()),2));
        deltat = abs(pmt_list[ipmt]->GetTime() - pmt_list[jpmt]->GetTime());

        // count this pmt pair as a valid pair and store the time diff if
        // deltat < threshold in ns
        // pmts closer than threshold in cm
        // pmt channels not in a cluster
        if ((deltat < fDeltaTThresh) and (tubedist < fDistThresh) and (!cluster)){


          // if cluster pair removal is enabled, check to see if
          // every tube between these two has been hit
          if (fClusterMode){
            int high_ccc,low_ccc;
            if (iccc > jccc){
              high_ccc = iccc;
              low_ccc = jccc;
            }else{
              high_ccc = jccc;
              low_ccc = iccc;
            }

            cluster = 1;
            for (int iclust=low_ccc+1;iclust<high_ccc;iclust++){
              if (ccc_hits[iclust] == 0){
                cluster = 0;
                break;
              }
            }
          } // end if cluster mode

          if (!cluster){
            pair_count++;
            deltat_list.push_back(deltat);
          }
        }
      } // end loop over j pmts
    } // end loop over i pmts


    // now we have constructed our list of deltats.
    // if we have enough statistics, compute the median.
    // fail the event if median is out of range
    if (pair_count >= fCountThresh){
      int middle = deltat_list.size()/2;
      nth_element(deltat_list.begin(),deltat_list.begin()+middle,deltat_list.end());
      if (deltat_list.size()%2){
        nth_element(deltat_list.begin(),deltat_list.begin()+middle-1,deltat_list.end());
        dtmedian = (deltat_list[middle]+deltat_list[middle-1])/2.0;
      }else{
        dtmedian = deltat_list[middle];
      }
      if (dtmedian > fMedianCut){
        fPassFlag = false;
      }
    }

    UpdateMask(ev);

    return fPassFlag ? OKFALSE : OKTRUE;
  }
} // namespace RAT