#include <RAT/FlasherGeoCut.hh>
#include <RAT/DB.hh>
#include <RAT/BitManip.hh>
#include <RAT/Processor.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>

using namespace std;

namespace RAT {

  void FlasherGeoCut::BeginOfRun(DS::Run&)
  {
    fLClean = DB::Get()->GetLink("DATACLEANING",fName);
    fMaxNhits = static_cast<unsigned int>( fLClean->GetI("max_nhits") );
    fMinNhits = static_cast<unsigned int>( fLClean->GetI("min_nhits") );
    fNhitsPsup = static_cast<unsigned int>( fLClean->GetI("nhits_psup") );
    fNhitsCrate = static_cast<unsigned int>( fLClean->GetI("nhits_crate") );
    fCCCLower = fLClean->GetI("ccc_lower");
    fCCCUpper = fLClean->GetI("ccc_upper");

    fRadiusCluster = fLClean->GetD("radius_cluster");
    fRatioCut = fLClean->GetD("ratio_cut");
    fTacCut = fLClean->GetD("tac_cut");
    fDistanceCut = fLClean->GetD("distance_cut");
  }

  Processor::Result FlasherGeoCut::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 FlasherGeoCut::Event(DS::Entry&, DS::EV& ev)
  {
    fPassFlag = true;

    // only attempt the cut if nhit is in the right range
    if (ev.GetNhits() > fMaxNhits || ev.GetNhits() < fMinNhits){
      UpdateMask(ev);
      return fPassFlag ? OKFALSE : OKTRUE;
    }

    double distance;
    int n_cluster;
    double t_cluster;
    int i_other;
    double t_other,dis_other;
    int bad_found;
    double ratio;
    double distance_max = -9999;
    float cqhs,cqhl,cqlx,uqhs,uqhl,uqlx;

    const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
    size_t pmt_count = pmtInfo.GetCount();
    vector<int> pmtids;
    vector<float> pmtx(pmt_count);
    vector<float> pmty(pmt_count);
    vector<float> pmtz(pmt_count);
    vector<float> pmtt(pmt_count);
    vector<int> badhits(pmt_count);
    vector<int> crates(pmt_count);
    vector<int> cards(pmt_count);

    int calibrated = 0;
    if (ev.GetCalPMTs().GetCount() == ev.GetUncalPMTs().GetCount()) // FIXME almost guaranteed to give the wrong answer
      calibrated = 1;

    int uncal_count = ev.GetUncalPMTs().GetCount();

    for (int i=0;i<uncal_count;i++){
      const DS::PMTUncal& upmt = ev.GetUncalPMTs().GetPMT(i);
      pmtids.push_back(upmt.GetID());
      pmtx[pmtids[i]] = pmtInfo.GetPosition(pmtids[i]).X();
      pmty[pmtids[i]] = pmtInfo.GetPosition(pmtids[i]).Y();
      pmtz[pmtids[i]] = pmtInfo.GetPosition(pmtids[i]).Z();
      pmtt[pmtids[i]] = upmt.GetTime();

      uqhs = upmt.GetQHS();
      uqhl = upmt.GetQHL();
      uqlx = upmt.GetQLX();
      if (calibrated){
        const DS::PMTCal& cpmt = ev.GetCalPMTs().GetPMT(i);
        cqhs = cpmt.GetQHS();
        cqhl = cpmt.GetQHL();
        cqlx = cpmt.GetQLX();
      }else{
        cqhs = -9999;
        cqhl = -9999;
        cqlx = -9999;
      }
      if ((abs(cqhs) > 8000 && uqhs < 250)
          || (abs(cqhs) < 8000 && cqhs < -100)
          || (uqhs > 4000)
          || (abs(cqhl) > 8000 && uqhl < 250)
          || (abs(cqhl) < 8000 && cqhl < -100)
          || (uqhl > 4000)
          || (abs(cqlx) > 8000 && uqlx < 300)
          || (abs(cqlx) < 8000 && cqlx < -100)
          || (uqlx > 4000)){
        badhits[pmtids[i]] = 1;
      }else{
        badhits[pmtids[i]] = 0;
      }
      crates[pmtids[i]] = BitManip::GetCrate(pmtids[i]);
      cards[pmtids[i]] = BitManip::GetCard(pmtids[i]);
    }

    // first we check for any clusters in real pmt
    // space. We check for bad charges and calculate
    // the distances
    for (int i=0;i<uncal_count;i++){
      n_cluster = 0;
      bad_found = 0;
      t_cluster = 0;
      i_other = 0;
      dis_other = 0;
      t_other = 0;
      for (int j=0;j<uncal_count;j++){
        distance = sqrt(pow(pmtx[pmtids[j]]-pmtx[pmtids[i]],2)
            +pow(pmty[pmtids[j]]-pmty[pmtids[i]],2)
            +pow(pmtz[pmtids[j]]-pmtz[pmtids[i]],2));
        if (distance < fRadiusCluster){
          n_cluster++;
          t_cluster += pmtt[pmtids[j]];
          if (badhits[pmtids[j]])
            bad_found = 1;
        }else{
          i_other++;
          dis_other += distance;
          t_other += pmtt[pmtids[j]];
        }
      }

      ratio = (double) (i_other) / (double) (i_other + n_cluster);
      if ((n_cluster >= fNhitsPsup) && (ratio > fRatioCut)){
        dis_other /= (double) i_other;
        t_other /= (double) i_other;
        t_cluster /= (double) n_cluster;
        if ((bad_found || t_cluster > (fTacCut + t_other))
            && (dis_other > distance_max)){
          distance_max = dis_other;
        }
      }
    } // end loop over i

    if (distance_max > fDistanceCut){
      fPassFlag = false;
      UpdateMask(ev);
      return fPassFlag ? OKFALSE : OKTRUE;
    }

    sort(pmtids.begin(), pmtids.end());

    // now we do the same thing looking
    // for clusters in electronics space
    distance_max = -9999;
    int cluster_found;
    float distance_other_max;
    int continueflag;
    for (int i=0;i<(uncal_count-fNhitsCrate+1);i++){
      continueflag = 0;
      cluster_found = 1;
      bad_found = 0;
      t_cluster = 0;
      for (int j=0;j<(fNhitsCrate-1);j++){
        if ((pmtids[i+j] != pmtids[i]+j) ||
            (crates[pmtids[i]] != crates[pmtids[i+j]]) ||
            (cards[pmtids[i]] != cards[pmtids[i+j]]))
          cluster_found = 0;
        t_cluster += pmtt[pmtids[i+j]];
        if (badhits[pmtids[i+j]])
          bad_found = 1;
      } // end loop over j pmts

      t_cluster = t_cluster / (double) fNhitsCrate;
      if (cluster_found){
        i_other = 0;
        t_other = 0;
        distance_other_max = 0;
        for (int j=0;j<fNhitsCrate;j++){
          dis_other = 0;
          for (int k=0;k<uncal_count;k++){
            if ((pmtids[k] < pmtids[i] - fCCCLower) ||
                (pmtids[k] > pmtids[i] + fCCCUpper + fNhitsCrate)){
              i_other++;
              t_other += pmtt[pmtids[k]];
              dis_other += sqrt(pow(pmtx[pmtids[k]]-pmtx[pmtids[i+j]],2)
            +pow(pmty[pmtids[k]]-pmty[pmtids[i+j]],2)
            +pow(pmtz[pmtids[k]]-pmtz[pmtids[i+j]],2));
            }
          } // loop over k

          ratio = (double) (i_other) / (double) (i_other + fNhitsCrate);
          if (ratio > fRatioCut)
            dis_other /= (double) i_other;
          else
            continueflag = 1;
          if (continueflag)
            break;

          if (dis_other > distance_other_max)
            distance_other_max = dis_other;
        } // loop over j

        if (continueflag)
          continue;

        t_other = t_other / (double) i_other;
        if (bad_found || (t_cluster > (fTacCut + t_other)))
          if (distance_other_max > distance_max)
            distance_max = distance_other_max;

      } // end if cluster found

    } // end loop over i pmts

    if (distance_max > fDistanceCut){
      fPassFlag = false;
      UpdateMask(ev);
      return fPassFlag ? OKFALSE : OKTRUE;
    }

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

} // namespace RAT