#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 lightwater_sno as default
  const string defaultIndex = "lightwater_sno";
  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;
    }

  fLightPath.CalcByPosition( startPos, DU::Utility::Get()->GetPMTInfo().GetPosition( pmt.GetID() ) );
  double distInInnerAV = fLightPath.GetDistInInnerAV();
  double distInAV = fLightPath.GetDistInAV();
  double distInWater = fLightPath.GetDistInWater();

  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 );
}