#include <CLHEP/Units/PhysicalConstants.h>
using CLHEP::pi; using CLHEP::twopi; using CLHEP::halfpi;

#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/DU/GroupVelocity.hh>
using namespace RAT;
using namespace RAT::DU;

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

ClassImp( RAT::DU::GroupVelocity )

void
GroupVelocity::LoadGroupVelocity( DBLinkPtr dbTable,
                                  TGraph& property )
{
  const double hbarc = 197.32705e-12;
  const double c_light = 299.792458;

  enum EWavelength { eWavelength, eDyWavelength, eEnergy };
  int wavelengthOption = eEnergy;
  const string option = dbTable->GetS( "RINDEX_option" );
  if( option == string( "wavelength" ) )
    wavelengthOption = eWavelength;
  else if( option == string( "dy_dwavelength" ) )
    wavelengthOption = eDyWavelength;
  else if( option == string( "energy" ) )
    wavelengthOption = eEnergy;

  vector<double> gvel_wavelengths;
  vector<double> gvel_velocities;
  try
  {
    gvel_wavelengths = dbTable->GetDArray( "GROUPVEL_value1" );
    gvel_velocities = dbTable->GetDArray( "GROUPVEL_value2" );

    warn << "GroupVelocity::LoadGroupVelocity: Explicit group velocities found for entry: " + dbTable->GetName() + "[" + dbTable->GetIndex() + "] - these will be used \n";
    Log::Assert( gvel_wavelengths.size() == gvel_velocities.size(), "GroupVelocity::LoadGroupVelocity: Value arrays are miss-sized for GROUPVEL in " + dbTable->GetName() + "[" + dbTable->GetIndex() + "]" );

    for (size_t w = 0; w < gvel_wavelengths.size(); w++)
    {
      double energy_MeV;
      if ( gvel_wavelengths[w] == 0.0 )
      {
        warn << "GroupVelocity::LoadGroupVelocity: Warning, zero wavelength provided!\n";
        energy_MeV = 0.0;
      }
      energy_MeV = (twopi * hbarc) / (gvel_wavelengths[w]  * 1.0e-06);

      double groupVel_MmPerNs = gvel_velocities[w];

      property.SetPoint( w, energy_MeV, groupVel_MmPerNs );
    }
  }

  catch ( DBNotFoundError& e )
  {
    warn << "GroupVelocity::LoadGroupVelocity: No explicit group velocities available for entry: " + dbTable->GetName() + "[" + dbTable->GetIndex() + "] - calculating using old method \n";

    vector<double> val1 = dbTable->GetDArray( "RINDEX_value1" );
    vector<double> val2 = dbTable->GetDArray( "RINDEX_value2" );

    Log::Assert( val1.size() == val2.size() && !val1.empty(), "Util::LoadPropertyVector: Value arrays are miss-sized for RINDEX in " + dbTable->GetName() + "[" + dbTable->GetIndex() + "]" );
    int start = 0, dir = +1;
    if( wavelengthOption != eEnergy && val1.front() < val1.back() ) // Must invert the order as energy is 1/wavelength
      {
        start = val1.size() - 1;
        dir = -1;
      }
    else if( wavelengthOption == eEnergy && val1.front() > val1.back() ) // Must invert the energy order
      {
        start = val1.size() - 1;
        dir = -1;
      }
    vector<double> energies, refIndicies;
    for( int index = start; index >= 0 && index < static_cast<int>( val1.size() ); index += dir )
      {
        double energy = val1[index];
        double value = val2[index];
        if( wavelengthOption == eDyWavelength )
          {
            energy = twopi * hbarc / ( val1[index] * 1e-6 );
            value *= val1[index] / energy; // Times value by the
          }
        else if( wavelengthOption == eWavelength )
          energy = twopi * hbarc / ( val1[index] * 1e-6 );
        energies.push_back( energy );
        refIndicies.push_back( value );
      }
    int point = 0;
    double e0 = energies[0];
    double n0 = refIndicies[0];
    double e1 = energies[1];
    double n1 = refIndicies[1];

    double vg = c_light/(n0+(n1-n0)/std::log(e1/e0));
    // allow only for 'normal dispersion' -> dn/d(logE) > 0
    //
    if((vg<0) || (vg>c_light/n0))  { vg = c_light/n0; }
    property.SetPoint( point++, e0, vg );
   for( size_t index = 2; index < energies.size(); index++ )
      {
        vg = c_light/( 0.5*(n0+n1)+(n1-n0)/std::log(e1/e0));
        // allow only for 'normal dispersion' -> dn/d(logE) > 0
        //
        if((vg<0) || (vg>c_light/(0.5*(n0+n1))))  { vg = c_light/(0.5*(n0+n1)); }
        property.SetPoint( point++, e0, vg );

        e0 = e1;
        n0 = n1;
        e1 = energies[index];
        n1 = refIndicies[index];
      }
    vg = c_light/(n1+(n1-n0)/std::log(e1/e0));
    // allow only for 'normal dispersion' -> dn/d(logE) > 0
    //
    if((vg<0) || (vg>c_light/n1))  { vg = c_light/n1; }
    property.SetPoint( point++, e0, vg );
  }
}

void
GroupVelocity::BeginOfRun()
{
  DB* db = DB::Get();
  try { LoadGroupVelocity( db->GetLink( "OPTICS", db->GetLink( "GEO", "inner_av" )->GetS( "material" ) ), fInnerAVGroupVelocity ); }
  catch( DBNotFoundError& e ) {
    try { LoadGroupVelocity( db->GetLink( "OPTICS", db->GetLink( "GEO", "inner_av" )->GetS( "material_top" ) ), fInnerAVGroupVelocity ); }
    catch ( DBNotFoundError& e ) { warn << "GroupVelocity::BeginOfRun: inner_av material refractive index cannot be loaded.\n"; }
  }
  try { LoadGroupVelocity( db->GetLink( "OPTICS", db->GetLink( "GEO", "av" )->GetS( "material" ) ), fAVGroupVelocity ); }
  catch( DBNotFoundError& e ) { warn << "GroupVelocity::BeginOfRun: av material refractive index cannot be loaded.\n"; }
  try { LoadGroupVelocity( db->GetLink( "OPTICS", db->GetLink( "GEO", "cavity" )->GetS( "material" ) ), fWaterGroupVelocity ); }
  catch( DBNotFoundError& e ) { warn << "GroupVelocity::BeginOfRun: cavity material refractive index cannot be loaded.\n"; }
}

double
GroupVelocity::PMTBucketTime( double incidentAngle )
{

  // If the incident angle is within the functional domain: [0,47]-degrees
  // then calculate the bucket time according
  // to a functional (degree-4) approximation.
  // If the incident angle is outside of this domain
  // (i.e. greater than 47-degrees) then approximate at 47-degrees
  // to avoid discontinuities in the return value.
  if ( incidentAngle > 47.0 || incidentAngle < 0.0 ){ incidentAngle = 47.0; }
  return 0.4543
    + 0.003148 * incidentAngle
    - 0.0002565 * pow( incidentAngle, 2 )
    + 9.626e-6 * pow( incidentAngle, 3 )
    - 6.865e-8 * pow( incidentAngle, 4 );
  // Details on where the above coefficients come from
  // are available in SNO+-doc-3138.

}