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

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

PMTSNOMANOpticalModel::PMTSNOMANOpticalModel( const G4String name,
                                              G4Region* const region,
                                              const std::string params,
                                              G4LogicalVolume* const pmtVolume,
                                              G4LogicalVolume* const envelopeVolume,
                                              G4OpticalSurface* const photocathodeSurface,
                                              G4OpticalSurface* const mirrorSurface )
  : PMTOpticalModel( name, region, params, pmtVolume, envelopeVolume, photocathodeSurface, mirrorSurface )
{
  fProcessName = "PMTSNOMANOpticalModel";
  // First load the overall model, Optics[something]...
  DBLinkPtr pmtModelTable = DB::Get()->GetLink( "PMT_OPTICAL_PARAMETERS", params );
  try
    {
      fA1 = pmtModelTable->GetD("a1");
      fA2 = pmtModelTable->GetD("a2");
      fA3 = pmtModelTable->GetD("a3");
      fB1 = pmtModelTable->GetD("b1");
      fB2 = pmtModelTable->GetD("b2");
    }
  catch( RAT::DBNotFoundError &e )
    {
      RAT::Log::Die( "PMTSNOMANOpticalModel::PMTSNOMANOpticalModel: Cannot Find PMT_OPTICAL_PARAMETERS Parameters" );
    }

  // Now rewrite the aluminium params
  fAlRIndex = new G4MaterialPropertyVector();
  fAlKIndex = new G4MaterialPropertyVector();
  fAlRIndex->InsertValues( twopi * hbarc/( 810 * nanometer ), 0.590321 );
  fAlRIndex->InsertValues( twopi * hbarc/( 60 * nanometer ), 0.590321 );
  fAlKIndex->InsertValues( twopi * hbarc/( 810 * nanometer ),  5.3618 );
  fAlKIndex->InsertValues( twopi * hbarc/( 60 * nanometer ),  5.3618 );

  fGlassRIndex = new G4MaterialPropertyVector();
  fGlassRIndex->InsertValues( twopi * hbarc/( 810 * nanometer ), 1.49 );
  fGlassRIndex->InsertValues( twopi * hbarc/( 60 * nanometer ), 1.49 );
}

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

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

  PMTOpticalModel::DoIt( fastTrack, fastStep );
}

G4double
PMTSNOMANOpticalModel::GetSignalProbability( G4int ipmt,
                                             const G4ThreeVector pos,
                                             const G4double energy,
                                             const G4double weight )
{
  // Based on pmt_resp_func.for
  const G4double kPMTTopOffset = 71.4; //Position of the inside glass on the z axis at the front of the pmt
  G4double signalProb = PMTOpticalModel::GetSignalProbability( ipmt, pos, energy, weight ) / 0.55; // NOTE 0.55 Here to match SNOMAN
  G4double distFromTop4 = ( pos.z() - kPMTTopOffset ) / 10.0; // 10.0 As MCE works in cm units
  distFromTop4 = distFromTop4 * distFromTop4 * distFromTop4 * distFromTop4;
  const G4double mce = (fA1 - fA2 * exp( fB1 * distFromTop4 ) ) * exp( fB2 * distFromTop4 ) + fA3;
  signalProb *= mce;

  return signalProb;
}