#include <RAT/QPDT.hh>
#include <RAT/DS/FitResult.hh>
#include <RAT/DS/FitVertex.hh>
#include <RAT/DB.hh>
#include <RAT/PMTCharge.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/DU/LightPathCalculator.hh>
#include <RAT/DU/EffectiveVelocity.hh>
using namespace RAT;
using namespace RAT::Classifiers;
using namespace RAT::DS;

#include <math.h>
using namespace std;

void
QPDT::BeginOfRun( DS::Run& )
{
  DBLinkPtr dbQPDTLink = DB::Get()->GetLink("CLASSIFIER_QPDT");
  fLowerTimeLimit = dbQPDTLink->GetD("t1");
  fUpperTimeLimit = dbQPDTLink->GetD("t2");
  fQCutoffMax = dbQPDTLink->GetD("qcutoff_max");
  fQCutoffMin = dbQPDTLink->GetD("qcutoff_min");
  fLightPath = DU::Utility::Get()->GetLightPathCalculator();
}

void
QPDT::DefaultSeed()
{
  fEventPos = TVector3();
  fEventTime = 0.0;
  fSeedVertex.SetPosition( fEventPos, false, true );
  fSeedVertex.SetTime( fEventTime, false, true );
}

void
QPDT::SetSeed( const DS::FitResult& seed )
{
  try
    {
      fSeedVertex = seed.GetVertex(0);
    }
  catch( FitResult::NoVertexError& error )
    {
      // Nothing to seed from, strange request
      return;
    }
  try { fEventPos = fSeedVertex.GetPosition(); }
  catch( FitVertex::NoValueError& error ) { /* No data to seed from */ }
  try { fEventTime = fSeedVertex.GetTime(); }
  catch( FitVertex::NoValueError& error ) { /* No data to seed from */ }
}

DS::ClassifierResult
QPDT::GetClassification()
{
  fClassifierResult.Reset();
  if( fPMTData.empty() )
    return fClassifierResult;

  SelectPMTData( fSeedVertex );

  int numEarly = 0; // Number of hits in early time window
  double qmax = 0; // Maximum charge in early time window
  double lowestProb=1.0; // Lowest probability hit in early time window, based on integrating qhs distribtuions

  const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
  // fSelectedPMTData uses the PMTCal selector
  for( vector<FitterPMT>::const_iterator iPMT = fSelectedPMTData.begin(); iPMT != fSelectedPMTData.end(); ++iPMT )
    {
      const TVector3 pmtPos = pmtInfo.GetPosition( iPMT->GetID() );

      fLightPath.CalcByPosition(fEventPos, pmtPos);
      double distInInnerAV = fLightPath.GetDistInInnerAV();
      double distInAV = fLightPath.GetDistInAV();
      double distInWater = fLightPath.GetDistInWater();

      const double transitTime = DU::Utility::Get()->GetEffectiveVelocity().CalcByDistance( distInInnerAV, distInAV, distInWater );
      const double corTime = iPMT->GetTime() - transitTime - fEventTime;

      if( corTime < fUpperTimeLimit && corTime > fLowerTimeLimit )
        {
          numEarly++;
          int q = iPMT->GetQHS();
          if (q < fQCutoffMax && q >= fQCutoffMin)
            {
              double integral1 = PMTCharge::Get()->PIntegrate(q, iPMT->GetID());
              double integral2 = PMTCharge::Get()->PIntegrate(fQCutoffMax, iPMT->GetID());
              double prob = 1. - (integral1 / integral2);
              if (prob < lowestProb)
                {lowestProb = prob;}
              if(q > qmax)
                {qmax = q;}
            }
          else
            {lowestProb=0;}
        }
    }
  double probCorrected = 1.0;
  if (numEarly > 0)
    {probCorrected = 1.0-pow((1.0-lowestProb),numEarly);}

  fClassifierResult.SetClassification( "QPDT", probCorrected);
  fClassifierResult.SetClassification("NHits_early",numEarly);
  fClassifierResult.SetClassification("QMax", qmax);
  fClassifierResult.SetValid( true );
  return fClassifierResult;
}