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

#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/RadialPDF.hh>
#include <RAT/FitterPMT.hh>
#include <RAT/DS/FitVertex.hh>
#include <RAT/ET1D.hh>
using namespace RAT;
using namespace RAT::PDFs;
using namespace RAT::DS;

RadialPDF::~RadialPDF()
{
  for( unsigned int uPDF = 0; uPDF < fPDFs.size(); uPDF++ )
    delete fPDFs[uPDF];
  fPDFs.clear();
}

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

void
RadialPDF::BeginOfRun( DS::Run& run )
{
  DB *db = DB::Get();
  string index;
  if( fIndex.empty() == true )
    {
      index = db->GetLink("GEO", "inner_av")->GetS("material");
    }
  else
    index = fIndex;

  RAT::DBLinkPtr radialSummaryTable = RAT::DB::Get()->GetLink( "RADIAL_PDFS", index );

  try
    {
      fRadialCentre = radialSummaryTable->GetDArray( "radial_centre" );
    }
  catch ( RAT::DBNotFoundError &e)
    {
      RAT::Log::Die( "RadialPDF::Initialise: Cannot Find PDF Parameters" );
    }

  for( unsigned int uPDF = 0; uPDF < fRadialCentre.size(); uPDF++ )
    {
      ET1D* newPDF = new ET1D();
      newPDF->BeginOfRun( run );
      stringstream pdfName;
      pdfName << "P" << uPDF;
      newPDF->Initialise( pdfName.str() );
      fPDFs.push_back( newPDF );
    }
}

void
RadialPDF::EndOfRun( DS::Run& run )
{
  for( size_t iPDF = 0; iPDF < fPDFs.size(); iPDF++ )
    {
      fPDFs[iPDF]->EndOfRun( run );
      delete fPDFs[iPDF];
    }
}

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

  unsigned int pdf = 0;
  for( unsigned int uPDF = 1; uPDF < fRadialCentre.size(); uPDF++ )
    if( startPos.Mag() > fRadialCentre[uPDF-1] && startPos.Mag() < fRadialCentre[uPDF] )
      pdf = uPDF;

  double probability = 1e-10;
  if( pdf == 0 || pdf == fRadialCentre.size() )
    probability = fPDFs[pdf]->GetProbability( pmt, vertex );
  else
    {
      const double prob1 = fPDFs[pdf]->GetProbability( pmt, vertex );
      const double prob0 = fPDFs[pdf - 1]->GetProbability( pmt, vertex );
      probability = ( startPos.Mag() - fRadialCentre[pdf - 1] ) * ( prob1 - prob0 ) / ( fRadialCentre[pdf] - fRadialCentre[pdf - 1] ) + prob0;
    }
  return probability;
}