#include <G4PhysicsOrderedFreeVector.hh>

#include <vector>
#include <string>
using namespace std;

#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/GV1D.hh>
#include <RAT/FitterPMT.hh>
#include <RAT/DS/FitVertex.hh>
#include <RAT/DS/Run.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/DU/LightPathCalculator.hh>
#include <RAT/DU/GroupVelocity.hh>
using namespace RAT;
using namespace RAT::PDFs;
using namespace RAT::DS;

void
GV1D::Initialise( const std::string& param )
{
  fIndex = param;
}

void
GV1D::BeginOfRun( DS::Run& )
{
  DB *db = DB::Get();
  string index;
  if( fIndex.empty() == true )
    {
      index = db->GetLink( "GEO", "inner_av" )->GetS( "material" );
    }
  else
    index = fIndex;
  // Set PDF to internalwater_snoplus as default
  const string defaultIndex = "internalwater_snoplus";
  RAT::DBLinkPtr promptTimingTable = RAT::DB::Get()->GetLink( "GV1D", defaultIndex );

  // then check for material specific PDF
  DBLinkGroup grp=db->GetLinkGroup( "GV1D" );
  DBLinkGroup::iterator it;
  for( it = grp.begin(); it != grp.end(); ++it )
    {
      if( index == it->first )
        {
          promptTimingTable = RAT::DB::Get()->GetLink("GV1D",index);
          break;
        }
    }

  vector<double> times;
  vector<double> probability;
  try
    {
	  times = promptTimingTable->GetDArray( "time" );
	  probability = promptTimingTable->GetDArray( "probability" );
	}
  catch ( RAT::DBNotFoundError &e) { RAT::Log::Die( "GV1D::Initialise: Cannot Find Prompt PDF Parameters" ); }

  fProbability = new G4PhysicsOrderedFreeVector( &(times[0]), &(probability[0]), times.size() );

  fLightPath = DU::Utility::Get()->GetLightPathCalculator();
}

void
GV1D::EndOfRun( DS::Run& )
{
  delete fProbability;
}

double
GV1D::GetProbability( const FitterPMT& pmt,
                      const FitVertex& vertex )
{
  TVector3 startPos;
  double eventTime;
  try
    {
      startPos = vertex.GetPosition();
      eventTime = vertex.GetTime();
    }
  catch( FitVertex::NoValueError& error )
    {
      warn << "GV1D::GetProbability: This PDF does nothing without a proposed position & time" << newline;
      return 0.0;
    }

  const double mag2 = startPos.Mag2();
  if( std::isinf( mag2 ) || std::isnan( mag2 ) )
    return 0.0;

  const TVector3 pmtPos = DU::Utility::Get()->GetPMTInfo().GetPosition( pmt.GetID() );
  double distInInnerAV = 0.0;
  double distInAV = 0.0;
  double distInWater = 0.0;

  fLightPath.CalcByPosition( startPos, pmtPos );
  if( fLightPath.GetPathValid() ) {
    distInInnerAV = fLightPath.GetDistInInnerAV();
    distInAV = fLightPath.GetDistInAV();
    distInWater = fLightPath.GetDistInWater();
  } else  // approximation: total distance in only inner AV material
    distInInnerAV = ( startPos - pmtPos ).Mag();

  const double transitTime = DU::Utility::Get()->GetGroupVelocity().CalcByDistance( distInInnerAV, distInAV, distInWater );
  const double correctedTime = pmt.GetTime() - transitTime - eventTime;

  if( correctedTime > -100.0 && correctedTime < 300.0 )
	return fProbability->Value( correctedTime );
  // "Punish" if way out of range
  else if( correctedTime < -100.0 )
	return exp( correctedTime + 100.0 ) * fProbability->Value( correctedTime );
  else
	return exp( 300.0 - correctedTime ) * fProbability->Value( correctedTime );
}