#include <TMath.h>
#include <TVector3.h>
using namespace ROOT;

#include <vector>
using namespace std;

#include <RAT/NearAVAngularClassifier.hh>
#include <RAT/DB.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
using namespace RAT;
using namespace RAT::Classifiers;

#include <CLHEP/Units/PhysicalConstants.h>  // required for the value of pi
using namespace CLHEP;


void NearAVAngularClassifier::BeginOfRun(DS::Run&)
{
  DB *db = DB::Get();

  // Parameters are the same for all materials, so just use the "labppo_scintillator" entry
  DBLinkPtr dbNAVLink = db->GetLink("CLASSIFIER_NEAR_AV_ANGULAR", "labppo_scintillator");

  fNhitsCut = dbNAVLink->GetI("nhits_cut");
  fDipRegionStart = dbNAVLink->GetD("dip_region_start");
  fDipRegionEnd = dbNAVLink->GetD("dip_region_end");
}

void NearAVAngularClassifier::EndOfRun(DS::Run&) {}

DS::ClassifierResult NearAVAngularClassifier::GetClassification()
{
  fClassifierResult.Reset();

  if (fPMTData.empty())
  {
    throw ClassifierFailError("NearAVAngularClassifier::GetClassification() - no Nhits in event - cannot calculate classifier");
  }

  // Classifier cannot reliably calculate ratio if not enough PMTs were hit
  if (fPMTData.size() < fNhitsCut)
  {
    return fClassifierResult;
  }

  // Calculate the reconstructed event axis w.r.t. to AV centre
  TVector3 eventAxis (0.0, 0.0, 0.0);
  const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();

  // First, get the 10th percerntile hit PMT
  vector<double> pmtTimes(fPMTData.size(), 0);
  for (vector<FitterPMT>::const_iterator iPMT = fPMTData.begin(); iPMT != fPMTData.end(); ++iPMT){
    size_t index = iPMT - fPMTData.begin();
    pmtTimes[index] = iPMT->GetTime();
  }

  int pmtStart = static_cast<int>(0.1 * fPMTData.size());
  nth_element(pmtTimes.begin(), pmtTimes.begin() + pmtStart, pmtTimes.end());
  double windowStart = pmtTimes[pmtStart] - 10; // ns
  double windowEnd = pmtTimes[pmtStart] + 10; // ns

  for (vector<FitterPMT>::const_iterator iPMT = fPMTData.begin(); iPMT != fPMTData.end(); ++iPMT)
    {
      double pmtTime = iPMT->GetTime();
      if(pmtTime >= windowStart && pmtTime <= windowEnd) {
        // Only include PMTs within the time window
        TVector3 pmtPos = pmtInfo.GetPosition(iPMT->GetID()).Unit();
        eventAxis += pmtPos;
      }
    }

  eventAxis = eventAxis.Unit(); // normalise

  // Calculate ratio: #PMTs in this event with PMT-AVcentre-event angle in dip region / total #PMTs in this event
  double counter = 0.0, total = 0.0;

  for (vector<FitterPMT>::const_iterator iPMT = fPMTData.begin(); iPMT != fPMTData.end(); ++iPMT)
  {
    const TVector3 pmtVector = pmtInfo.GetPosition(iPMT->GetID());
    double angle = acos(eventAxis.Dot(pmtVector) / (eventAxis.Mag() * pmtVector.Mag())) * (180 / pi);

    if ((angle >= fDipRegionStart) && (angle <= fDipRegionEnd))
    {
      counter += 1.0;
    }
    total += 1.0;
  }
  double ratio = (counter / total);

  fClassifierResult.SetClassification("ratio", ratio);
  fClassifierResult.SetValid(true);
  return fClassifierResult;
}