#include <RAT/Beta14.hh>
#include <RAT/DB.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/FitVertex.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
using namespace RAT;
using namespace RAT::Classifiers;
using namespace RAT::DS;

using namespace std;

void
Beta14::DefaultSeed()
{
  fEventPos = TVector3();
}

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

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

  const size_t numPMTHits = fPMTData.size();
  double p1 = 0.0;
  double p4 = 0.0;
  const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
  for( size_t ipmt = 0; ipmt < numPMTHits - 1; ipmt++ )
    {
      const TVector3 pmtDir1 = ( pmtInfo.GetPosition( fPMTData[ipmt].GetID() ) - fEventPos ).Unit();
      for( size_t ipmt2 = ipmt + 1; ipmt2 < numPMTHits; ipmt2++ )
        {
          const TVector3 pmtDir2 = ( pmtInfo.GetPosition( fPMTData[ipmt2].GetID() ) - fEventPos ).Unit();
          const double cosTheta12 = pmtDir1.Dot( pmtDir2 );
          p1 += cosTheta12; // Legendre Polynomial order 1
          const double cosThetaSquared = cosTheta12 * cosTheta12;
          p4 += ( 35.0 * cosThetaSquared * cosThetaSquared - 30.0 * cosThetaSquared + 3.0 ) / 8.0; // Legendre Polynomial order 4
        }
    }
  p1 = 2.0 * p1 / static_cast<double>( numPMTHits * ( numPMTHits - 1 ) );
  p4 = 2.0 * p4 / static_cast<double>( numPMTHits * ( numPMTHits - 1 ) );
  const double beta14 = p1 + 4 * p4;

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