#include <G4PhysicsOrderedFreeVector.hh>

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

#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/ET1D.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/EffectiveVelocity.hh>
using namespace RAT;
using namespace RAT::PDFs;
using namespace RAT::DS;

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

void
ET1D::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 labppo_scintillator as default
  string defaultIndex = "labppo_scintillator";
  RAT::DBLinkPtr promptTimingTable = RAT::DB::Get()->GetLink( "ET1D",defaultIndex );

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

  vector<double> times;
  vector<double> probability;

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

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

  fLightPath = DU::Utility::Get()->GetLightPathCalculator();
}
void
ET1D::EndOfRun( DS::Run& )
{
  delete fProbability;
}

double
ET1D::GetProbability( const FitterPMT& pmt,
                      const FitVertex& vertex )
{
  TVector3 startPos;
  double eventTime;
  try
    {
      startPos = vertex.GetPosition();
      eventTime = vertex.GetTime();
    }
  catch( FitVertex::NoValueError& error )
    {
      warn << "ET1D::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()->GetEffectiveVelocity().CalcByDistance( distInInnerAV, distInAV, distInWater );
  const double correctedTime = pmt.GetTime() - transitTime - eventTime;

  return fProbability->Value( correctedTime );
}