#include <G4Material.hh>
#include <G4MaterialPropertyVector.hh>
#include <G4PhysicsOrderedFreeVector.hh>

#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/SimpleTiming.hh>
#include <RAT/FitterPMT.hh>
#include <RAT/DS/FitVertex.hh>
#include <CLHEP/Units/SystemOfUnits.h>
#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 CLHEP;
using namespace RAT::DS;

#include <string>
using namespace std;

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

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

  G4MaterialPropertyVector* theWaveForm = G4Material::GetMaterial( material )->GetMaterialPropertiesTable()->GetProperty( "SCINTWAVEFORM" );
  const double maxTime = -3.0*(theWaveForm->GetMinLowEdgeEnergy()); // Smallest decay constant, window is three times this
  const double binWidth = 1.0; // in ns for simplicity
  const int pBins = static_cast<int>( 50.0 / binWidth ); // 50ns of bins before 0
  const int nBins = static_cast<int>( maxTime / binWidth );
  double* tval = new double[nBins + pBins]; // Time values
  double* pval = new double[nBins + pBins]; // Probability values
  for( int iLoop = 0; iLoop < nBins + pBins; iLoop++ )
    {
      tval[iLoop] = ( iLoop - pBins ) * binWidth;
      pval[iLoop] = 0.0;
    }

  for (unsigned int iref = 1;
       iref < theWaveForm->GetVectorLength();iref++)
  {
      const double ampl = (*theWaveForm)[iref];     // Decay component fraction/amplitude
      const double decy = theWaveForm->Energy(iref); // Decay constant
      for( int iLoop = pBins; iLoop < nBins + pBins; iLoop++ ) // Only start from 0ns (in pBins)
        {
          double temp = ampl * exp( tval[iLoop] / decy ) * ( 1.0 - exp( binWidth / decy ) ); // Fraction at this time for this component
          temp /=  theWaveForm->GetVectorLength() * ampl * ( 1.0 - exp( maxTime / decy ) ); // Normalisation
          pval[iLoop] += temp;
        }
      }

  for( int iLoop = 0; iLoop < pBins; iLoop++ )
    pval[iLoop] =  pval[pBins] * exp( -tval[iLoop] * tval[iLoop] / ( 2.0 * 2.5 ) ); // sqrt(2.5) bin width spread

  // Now normalise
  double normalisation = 0.0;
  for( int iLoop = 0; iLoop < pBins; iLoop++ )
    normalisation += pval[iLoop];
  for( int iLoop = 0; iLoop < pBins; iLoop++ )
    pval[iLoop] /= normalisation;

  fProbaility = new G4PhysicsOrderedFreeVector( tval, pval, nBins + pBins );

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

  // in Geant4.0.0, G4PhysicsOrderedFreeVector makes its own copy
  // of any array passed to its constructor, so ...
  delete[] tval;
  delete[] pval;

  return;
}

double
SimpleTiming::GetProbability( const FitterPMT& pmt,
                              const FitVertex& vertex )
{
  TVector3 startPos;
  double eventTime;
  try
    {
      startPos = vertex.GetPosition();
      eventTime = vertex.GetTime();
    }
  catch( FitVertex::NoValueError& error )
    {
      warn << "SimpleTiming::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;
  return fProbaility->Value( correctedTime );
}