#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
using namespace CLHEP;

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

#include <string>
using namespace std;

ConcentratorSNOMANOpticalModel::ConcentratorSNOMANOpticalModel( const G4String name,
                                                                G4Region* const region,
                                                                const std::string params,
                                                                G4LogicalVolume* const volume )
  : ConcentratorOpticalModel( name, region, params, volume )
{
  fProcessName = "ConcentratorSNOMANOpticalModel";

  // First load the overall model, Optics[something]...
  DBLinkPtr concModelTable = DB::Get()->GetLink( "CONC_OPTICAL_PARAMETERS", params );
  string promptParams;
  try
    {
      promptParams = concModelTable->GetS( "prompt_parameters" );
    }
  catch( RAT::DBNotFoundError &e )
    {
      RAT::Log::Die( "ConcentratorSNOMANOpticalModel::ConcentratorSNOMANOpticalModel: Cannot Find CONC_OPTICAL_PARAMETERS Parameters" );
    }

  // Now load the correct prompt parameters
  DBLinkPtr concPromptTable = DB::Get()->GetLink( "CONC_PROMPT_PARAMETERS", promptParams.c_str() );
  try
    {
      fPrompt = concPromptTable->GetDArray( "prompt_parameters" );
      fPromptNormalisation = concPromptTable->GetD( "prompt_normalisation" );
    }
  catch ( RAT::DBNotFoundError &e)
    {
      RAT::Log::Die( "ConcentratorSNOMANOpticalModel::ConcentratorSNOMANOpticalModel: Cannot Find CONC_PROMPT_PARAMETERS Parameters" );
    }
}

void
ConcentratorSNOMANOpticalModel::DoIt( const G4FastTrack& fastTrack,
                                      G4FastStep& fastStep )
{
  G4ThreeVector dir =  fastTrack.GetPrimaryTrackLocalDirection();
  G4int trackID = fastTrack.GetPrimaryTrack()->GetTrackID();
  G4int ipmt = GetPMTID( fastTrack );

  const G4ThreeVector bucketNorm( 0.0, 0.0, 1.0 );
  double cosBucket = dir * bucketNorm;
  OpticalModelBase::SetBucketCos( trackID, ipmt, cosBucket );

  ConcentratorOpticalModel::DoIt( fastTrack, fastStep );
}

void
ConcentratorSNOMANOpticalModel::GetReflectionProb( const G4ThreeVector& position,
                                                   const G4double theta, // In Radians
                                                   const G4double energy, // In MeV
                                                   const EPolarisation polarisation, // As enum defined in OpticalModelBase.hh
                                                   G4double& specularR,
                                                   G4double& diffuseR )
{
  const G4double wavelength = twopi * hbarc / energy;
  // Modelled on pmt_ref_alum.for
  ConcentratorOpticalModel::GetReflectionProb( position, theta, energy, polarisation, specularR, diffuseR );
  const G4int iBucketTheta = static_cast<G4int>( acos( OpticalModelBase::GetBucketCos() ) / degree + 0.5 );

  G4int iWavelength = static_cast<G4int>( ( wavelength * mm ) / ( 10.0 * nm ) + 0.5 ); // Want in units of 10nm
  if( iWavelength >= 23 && iWavelength <= 70 )
    {
      iWavelength -= 23;
      specularR *= fPromptNormalisation / ( 1.0 + fPrompt[ iBucketTheta + iWavelength * 91] * ( 1.0 - fPromptNormalisation ) );
    }
  specularR = specularR - diffuseR;
}