#include <Randomize.hh>

#include <RAT/PMTTransitTime.hh>
#include <RAT/DB.hh>
using namespace RAT;

#include <vector>
using namespace std;

void
PMTTransitTime::BeginOfRun()
{
  // fixme
  DBLinkPtr transitDB = DB::Get()->GetLink( "PMTTRANSIT", "r1408" );
  fX = transitDB->GetDArray( "x" );
  fY = transitDB->GetDArray( "y" );
  fCableDelay = transitDB->GetD( "cable_delay" );

  DBLinkPtr transitDB_hqe = DB::Get()->GetLink( "PMTTRANSIT", "r5912" );
  fX_hqe = transitDB_hqe->GetDArray( "x" );
  fY_hqe = transitDB_hqe->GetDArray( "y" );

}

double
PMTTransitTime::CalculateFrontEndTime( const double PETime, const unsigned int pmtID )
{

  isHQE(pmtID); // Check whether the PMT if HQE to determine what timing to use

  const double rand = G4UniformRand();

  size_t up = 1;
  for( size_t i = 1; i < time.size(); i++ )
    {
      if ( rand < cdf[i] )
        {
          up = i;
          i = time.size();
        }
    }
  const double transitTime = ( rand - cdf[(up-1)] ) * ( time[up] - time[(up-1)] ) / ( cdf[up] - cdf[(up-1)] ) + time[(up-1)];
  return PETime + transitTime + fCableDelay;
}

void
PMTTransitTime::isHQE(const unsigned pmtID){

  const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
  if(pmtInfo.GetType(pmtID) == DU::PMTInfo::HQE){
     time = fX_hqe;
     cdf  = fY_hqe;
  }
  else{
     time = fX;
     cdf  = fY;
  }
}


void
PMTTransitTime::GetTimeLimits(double& smallest, double& largest) {
    smallest = fX[0] + fCableDelay;
    largest =  fX[fX.size()-1] + fCableDelay;
}

double
PMTTransitTime::GetIntervalProbability( double ta, double tb, bool cableDelay){
    // cableDelay=true (the default) to indicate that the input times ta & tb
    // include cable delay
    if(tb<ta) {
        double temp = ta;
        ta = tb;
        tb = temp;
    }
    if(cableDelay) {
        ta = ta - fCableDelay;
        tb = tb - fCableDelay;
    }
    size_t nLast = fX.size() - 1;
    if( tb<fX[0] || ta>fX[nLast] ) return 0.;
    if( ta<fX[0] ) ta = fX[0];
    if( tb>fX[nLast] ) tb = fX[nLast];

    int i = 0;
    while (fX[i+1] <= ta) { i++; }
    // Cumulative probability value at ta:
    double cumPa = fY[i] + (fY[i+1] - fY[i])*(ta-fX[i])/(fX[i+1]-fX[i]);

    while (fX[i+1] < tb) { i++;}
    // Cumulative probability value at tb:
    double cumPb = fY[i] + (fY[i+1] - fY[i])*(tb-fX[i])/(fX[i+1]-fX[i]);

    return (cumPb - cumPa);
}