#include <G4OpticalPhoton.hh>
#include <G4OpticalSurface.hh>
#include <G4TransportationManager.hh>
#include <Randomize.hh>
#include <G4Region.hh>
#include <G4UnitsTable.hh>
#include <G4GeometryTolerance.hh>
#include <local_g4compat.hh>

#include <RAT/PMTOpticalModel.hh>
#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/PhotonThinning.hh>
#include <RAT/Optics.hh>
#include <RAT/HitPMTCollection.hh>
#include <RAT/PMTVariation.hh>
#include <RAT/DS/MCPE.hh>
using namespace RAT;

using std::vector;
using std::complex;
using CLHEP::twopi;
using CLHEP::hbarc;

PMTOpticalModel::PMTOpticalModel( const G4String name,
                                  G4Region* const region,
                                  const std::string params,
                                  G4LogicalVolume* const pmtVolume,
                                  G4LogicalVolume* const envelopeVolume,
                                  G4OpticalSurface* const photocathodeSurface,
                                  G4OpticalSurface* const mirrorSurface )
  : G4VFastSimulationModel( name, region )
{
  fProcessName = "PMTOpticalModel";

  if( pmtVolume->GetNoDaughters() != 3 )
    Log::Die( "PMTOpticalModel::PMTOpticalModel Setup on a pmt without 3 inner volumes" );

  fGlassVolume = pmtVolume->GetSolid();
  fPCVacuumVolume = pmtVolume->GetDaughter(0)->GetLogicalVolume()->GetSolid();
  fPCAlVacuumVolume = pmtVolume->GetDaughter(1)->GetLogicalVolume()->GetSolid();
  fAlVacuumVolume = pmtVolume->GetDaughter(2)->GetLogicalVolume()->GetSolid();
  fDynodeVolume = pmtVolume->GetDaughter(2)->GetLogicalVolume()->GetDaughter(0)->GetLogicalVolume()->GetSolid(); // Assumes dynode is in the AlVacuum volume.
  fDynodeZTranslation = pmtVolume->GetDaughter(2)->GetLogicalVolume()->GetDaughter(0)->GetTranslation().z(); // Dynode is in the Al coord system not PMT

  fPMTVolNames.push_back( pmtVolume->GetName() );
  fPMTVolNames.push_back( pmtVolume->GetDaughter(0)->GetName() );
  fPMTVolNames.push_back( pmtVolume->GetDaughter(1)->GetName() );
  fPMTVolNames.push_back( pmtVolume->GetDaughter(2)->GetName() );
  fPMTVolNames.push_back( pmtVolume->GetDaughter(2)->GetLogicalVolume()->GetDaughter(0)->GetName() );
  fPMTVolNames.push_back( pmtVolume->GetName() + "_out" );

  // Now the Optical Parameters
  fWaterRIndex = envelopeVolume->GetMaterial()->GetMaterialPropertiesTable()->GetProperty("RINDEX");
  fGlassRIndex = pmtVolume->GetMaterial()->GetMaterialPropertiesTable()->GetProperty("RINDEX");
  fPCRIndex = photocathodeSurface->GetMaterialPropertiesTable()->GetProperty("RINDEX");
  fPCKIndex = photocathodeSurface->GetMaterialPropertiesTable()->GetProperty("KINDEX");
  fPCThickness = photocathodeSurface->GetMaterialPropertiesTable()->GetProperty("THICKNESS");
  fAlRIndex = mirrorSurface->GetMaterialPropertiesTable()->GetProperty("RINDEX");
  fAlKIndex = mirrorSurface->GetMaterialPropertiesTable()->GetProperty("KINDEX");
  fGlassGVelocity = pmtVolume->GetMaterial()->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");
  fPCVacuumGVelocity = pmtVolume->GetDaughter(0)->GetLogicalVolume()->GetMaterial()->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");
  fPCAlVacuumGVelocity = pmtVolume->GetDaughter(1)->GetLogicalVolume()->GetMaterial()->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");
  fAlVacuumGVelocity = pmtVolume->GetDaughter(2)->GetLogicalVolume()->GetMaterial()->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");

  // Load from DB, Optics[something]...
  DBLinkPtr pmtModelTable = DB::Get()->GetLink( "PMT_OPTICAL_PARAMETERS", params );
  fCollectionEfficiency = pmtModelTable->GetD( "collection_efficiency" );
  fDynodeReflectivity = pmtModelTable->GetD( "stack_reflectivity" );
  fPEEscapeProbability = Optics::LoadWavelengthPropertyVector( "pe_escape_probability", pmtModelTable );

  fPCThickness = Optics::LoadPropertyVector( "pc_thickness", pmtModelTable );
}

PMTOpticalModel::~PMTOpticalModel()
{

}

G4bool
PMTOpticalModel::IsApplicable( const G4ParticleDefinition& particleType )
{
  return ( &particleType == G4OpticalPhoton::OpticalPhotonDefinition() );
}

G4bool
PMTOpticalModel::ModelTrigger( const G4FastTrack& )
{
  // Tracks all photons in the PMT
  return true;
}

void
PMTOpticalModel::DoIt( const G4FastTrack& fastTrack,
                       G4FastStep& fastStep )
{

  // Assume Geant4 has properly passed a track that respects the region set
  // when placing PMTs

  G4ThreeVector position = fastTrack.GetPrimaryTrackLocalPosition();
  G4ThreeVector direction = fastTrack.GetPrimaryTrackLocalDirection();
  G4ThreeVector polarisation = fastTrack.GetPrimaryTrackLocalPolarization();
  G4double time = fastTrack.GetPrimaryTrack()->GetGlobalTime();
  G4double weight = fastTrack.GetPrimaryTrack()->GetWeight();
  G4int pmtid = GetPMTID(fastTrack);

  TrackPhoton(fastTrack, pmtid, position, direction, polarisation, time, weight);

  fastStep.SetPrimaryTrackFinalPosition( position );
  fastStep.SetPrimaryTrackFinalTime( time );
  fastStep.SetPrimaryTrackFinalMomentum( direction );
  fastStep.SetPrimaryTrackFinalPolarization( polarisation );
  if( weight <= 0.0 )
    fastStep.ProposeTrackStatus( fStopAndKill );
}

void
PMTOpticalModel::TrackPhoton( const G4FastTrack& fastTrack,
                              const G4int ipmt,
                              G4ThreeVector& position,
                              G4ThreeVector& direction,
                              G4ThreeVector& polarisation,
                              G4double& time,
                              G4double& weight )
{
  // Logic Summary
  // 1) Find where the track is inside the PMT
  // 2) Track to the next volume
  // 3) On crossing a border call the TrackAcrossBorder routine
  // 4) Track to the next volume, etc...
  // 5) If tracks out of the PMT glass then return to Geant4

  const G4int trackID = fastTrack.GetPrimaryTrack()->GetTrackID();
  const G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();

  // group velocities
  const G4double gvGlass = fGlassGVelocity->Value( energy );
  const G4double gvPCVacuum = fPCVacuumGVelocity->Value( energy );
  const G4double gvPCAlVacuum = fPCAlVacuumGVelocity->Value( energy );
  const G4double gvAlVacuum = fAlVacuumGVelocity->Value( energy );

  EPMTVolumes currentVolume, nextVolume;

  if (fGlassVolume->Inside(position) == kInside) {
    // These volumes are all inside glass, check them first before assuming glass
    if (fPCVacuumVolume->Inside(position) == kInside) {
      currentVolume = ePCVacuum;
    } else if (fPCAlVacuumVolume->Inside(position) == kInside) {
      currentVolume = ePCAlVacuum;
    // Dynode is insie AlVacuum, so must check it first.
    // Also dynode is offset from PMT local coordinates.
    } else if (fDynodeVolume->Inside(position - G4ThreeVector( 0.0, 0.0, fDynodeZTranslation )) == kInside) {
      currentVolume = eDynode;
    } else if (fAlVacuumVolume->Inside(position) == kInside) {
      currentVolume = eAlVacuum;
    } else {
      currentVolume = eGlass;
    }
  } else {
    currentVolume = eOut;
  }

  RecordTrajectory(fastTrack, position, direction, 0.0, time, energy, fPMTVolNames[currentVolume], fPMTVolNames[currentVolume]);

  int iIterations;
  for ( iIterations = 0; iIterations < fMaxIterations; iIterations++ ) {

    // Identify and move to next boundary based on curent location
    switch (currentVolume) {
      case ePCVacuum: {
        double dist = fPCAlVacuumVolume->DistanceToIn( position, direction ); // Track to surface
        if (dist == kInfinity) { // Must be going back into glass
          nextVolume = eGlass;
          dist = fPCVacuumVolume->DistanceToOut( position, direction ); // Track to surface
        } else {
          nextVolume = ePCAlVacuum;
        }
        //Advance the photon to the next boundary
        position += direction * dist;
        time += dist / gvPCVacuum;
        break;
      }
      case ePCAlVacuum: {
        vector<G4double> distances(2); // Store distances to the volumes we can potentially hit
        const EPMTVolumes volume[2] = { ePCVacuum, eAlVacuum }; //PCAlVacuum touches these
        distances[0] = fPCVacuumVolume->DistanceToIn( position, direction ); // Track to surface
        distances[1] = fAlVacuumVolume->DistanceToIn( position, direction ); // Track to surface
        uint32_t closest = min_element( distances.begin(), distances.end() ) - distances.begin(); //index of the smallest distance
        double dist = distances[closest];
        if (dist == kInfinity) { // Must be going back into glass
          nextVolume = eGlass;
          dist = fPCAlVacuumVolume->DistanceToOut( position, direction ); // Track to surface
        } else {
          nextVolume = volume[closest];
        }
        //Advance the photon to the next boundary
        position += direction * dist;
        time += dist / gvPCAlVacuum;
        break;
      }
      case eAlVacuum: {
        vector<G4double> distances(2); // Store distances to the volumes we can potentially hit
        const EPMTVolumes volume[2] = { ePCAlVacuum, eDynode }; //AlVacuum touches these
        distances[0] = fPCAlVacuumVolume->DistanceToIn( position, direction ); // Track to surface
        G4ThreeVector dynodePosition = position - G4ThreeVector( 0.0, 0.0, fDynodeZTranslation );
        distances[1] = fDynodeVolume->DistanceToIn( dynodePosition, direction ); // Track to surface
        uint32_t closest = min_element( distances.begin(), distances.end() ) - distances.begin(); //index of the smallest distance
        double dist = distances[closest];
        if (dist == kInfinity) { // Must be going back into glass
          nextVolume = eGlass;
          dist = fAlVacuumVolume->DistanceToOut( position, direction ); // Track to surface
        } else {
          nextVolume = volume[closest];
        }
        //Advance the photon to the next boundary
        position += direction * dist;
        time += dist / gvAlVacuum;
        break;
      }
      case eGlass: {
        vector<G4double> distances(4); // Store distances to the volumes we can potentially hit
        const EPMTVolumes volume[4] = { eOut, ePCVacuum, ePCAlVacuum, eAlVacuum }; //Glass touches these
        distances[0] = fGlassVolume->DistanceToOut( position, direction ); // Track to surface - eOut is special
        distances[1] = fPCVacuumVolume->DistanceToIn( position, direction ); // Track to surface
        distances[2] = fPCAlVacuumVolume->DistanceToIn( position, direction ); // Track to surface
        distances[3] = fAlVacuumVolume->DistanceToIn( position, direction ); // Track to surface
        uint32_t closest = min_element( distances.begin(), distances.end() ) - distances.begin(); //index of the smallest distance
        nextVolume = volume[closest]; // What we're about to hit
        //Advance the photon to the next boundary
        position += direction * distances[closest];
        time += distances[closest] / gvGlass;
        break;
      }
      case eDynode: {
        //Dynode is opaque, kill track
        RecordTrajectory(fastTrack, position, direction, 0.0, time, energy, fPMTVolNames[currentVolume], fPMTVolNames[currentVolume]);
        weight = 0.0;
        return;
      }
      case eOut: {
        double dist = fGlassVolume->DistanceToIn( position, direction );
        if (dist == kInfinity) { //never enters, or is exiting
          RecordTrajectory(fastTrack, position, direction, 0.0, time, energy, fPMTVolNames[currentVolume], fPMTVolNames[currentVolume]);
          return;
        } else {
          nextVolume = eGlass;
          //Advance the photon to the next boundary
          position += direction * dist;
          time += dist / gvGlass;
        }
        break;
      }
      default: {
        Log::Die("Impossible volume in PMTOpticalModel::TrackPhoton");
      }
    }

    RecordTrajectory(fastTrack, position, direction, 0.0, time, energy, fPMTVolNames[currentVolume], fPMTVolNames[currentVolume]);

    //Figure out what happens at the identified boundary
    bool reflected = TrackAcrossBorder( fastTrack, position, energy, time, currentVolume, nextVolume, ipmt, trackID, direction, polarisation, weight );

    if( weight == 0.0 ) { // weight set to 0 when absorbed
      RecordTrajectory(fastTrack, position, direction, 0.0, time, energy, fPMTVolNames[currentVolume], fPMTVolNames[nextVolume]);
      return;
    }

    if (reflected) {
        RecordTrajectory(fastTrack, position, direction, 0.0, time, energy, fPMTVolNames[currentVolume], fPMTVolNames[currentVolume]);
        //never left current volume
    } else {
        RecordTrajectory(fastTrack, position, direction, 0.0, time, energy, fPMTVolNames[currentVolume], fPMTVolNames[nextVolume]);
        currentVolume = nextVolume;
    }

  }

  if( iIterations == fMaxIterations ) {
    debug << "PMTOpticalModel::TrackPhoton, photon is stuck. Killing." << newline;
    weight = 0.0;
  }

}

bool
PMTOpticalModel::TrackAcrossBorder( const G4FastTrack& fastTrack,
                                    const G4ThreeVector& position,
                                    const G4double energy,
                                    const G4double time,
                                    const EPMTVolumes previousVolume,
                                    const EPMTVolumes currentVolume,
                                    const G4int ipmt,
                                    const G4int trackID,
                                    G4ThreeVector& direction,
                                    G4ThreeVector& polarisation,
                                    G4double& weight )
{
  // Logic Summary
  // 1) Establish which border was crossed, then track accordingly

  G4double n1;
  complex<double> n2, n3;
  G4double T = 0.0, R = 0.0;
  G4ThreeVector normal;

  if( currentVolume == ePCVacuum && previousVolume == eGlass ) // From glass to PC
    {
      n1 = fGlassRIndex->Value( energy );
      n2 = complex<double>( fPCRIndex->Value( energy ), fPCKIndex->Value( energy ) );
      n3 = complex<double>( 1.0, 0.0 );
      normal = fPCVacuumVolume->SurfaceNormal( position ).unit();
      const double pcThickness = fPCThickness->Value( position.z() );
      CalculateTRThinFilm( direction, polarisation, normal, n1, n2, n3, energy, pcThickness, T, R );
    }
  else if( currentVolume == eGlass && previousVolume == ePCVacuum ) // From PC to glass
    {
      n1 = 1.0;
      n2 = complex<double>( fPCRIndex->Value( energy ), fPCKIndex->Value( energy ) );
      n3 = complex<double>( fGlassRIndex->Value( energy ), 0.0 );
      normal = -fPCVacuumVolume->SurfaceNormal( position ).unit(); // Minus sign such that norm dot direction is < 0
      const double pcThickness = fPCThickness->Value( position.z() );
      CalculateTRThinFilm( direction, polarisation, normal, n1, n2, n3, energy, pcThickness, T, R );
    }
  else if( currentVolume == ePCAlVacuum && previousVolume == eGlass ) // From glass to PCAl
    {
      n1 = fGlassRIndex->Value( energy );
      n2 = complex<double>( fAlRIndex->Value( energy ), fAlKIndex->Value( energy ) );
      normal = fPCAlVacuumVolume->SurfaceNormal( position ).unit();
      CalculateTRFresnel( direction, polarisation, normal, n1, n2, T, R );
      T = 0; // No Transmission for this interface
    }
  else if( currentVolume == eGlass && previousVolume == ePCAlVacuum ) // From PCAl to glass
    {
      n1 = 1.0;
      n2 = complex<double>( fPCRIndex->Value( energy ), fPCKIndex->Value( energy ) );
      n3 = complex<double>( fAlRIndex->Value( energy ), fAlKIndex->Value( energy ) );
      normal = -fPCAlVacuumVolume->SurfaceNormal( position ).unit(); // Minus sign such that normal dot direction is < 0
      const double pcThickness = fPCThickness->Value( position.z() );
      CalculateTRThinFilm( direction, polarisation, normal, n1, n2, n3, energy, pcThickness, T, R );
    }
  else if( currentVolume == eAlVacuum && previousVolume == eGlass ) // From glass to Al
    {
      n1 = fGlassRIndex->Value( energy );
      n2 = complex<double>( fAlRIndex->Value( energy ), fAlKIndex->Value( energy ) );
      normal = fAlVacuumVolume->SurfaceNormal( position ).unit();
      CalculateTRFresnel( direction, polarisation, normal, n1, n2, T, R );
      T = 0; // No Transmission for this interface
      }
  else if( currentVolume == eGlass && previousVolume == eAlVacuum ) // From Al to glass
    {
      n1 = 1.0;
      n2 = complex<double>( fAlRIndex->Value( energy ), fAlKIndex->Value( energy ) );
      normal = -fAlVacuumVolume->SurfaceNormal( position ).unit(); // Minus sign such that normal dot direction is < 0
      CalculateTRFresnel( direction, polarisation, normal, n1, n2, T, R );
      T = 0; // No Transmission for this interface
    }
  else if( currentVolume == eDynode && previousVolume == eAlVacuum ) // From Al to dynode
    {
      n1 = 1.0;
      n2 = complex<double>( fAlRIndex->Value( energy ), fAlKIndex->Value( energy ) );
      normal = fDynodeVolume->SurfaceNormal( position ).unit();
      R = fDynodeReflectivity;
      T = 0; // No Transmission for this interface
    }
  else if( currentVolume == eOut && previousVolume == eGlass ) // Leaving
    {
      n1 = fGlassRIndex->Value( energy );
      n3 = complex<double>( fWaterRIndex->Value( energy ), 0.0 );
      normal = -fGlassVolume->SurfaceNormal( position ).unit(); // Minus sign such that normal dot direction is < 0
      CalculateTRFresnel( direction, polarisation, normal, n1, n3, T, R );
    }
  else
    return false; // Nothing to do, crossing between vacuums, thus return

  G4double randR = G4UniformRand();

  if( randR < R ) // Reflect
    {
      NormalReflect( normal, direction, polarisation );
      return true;
    }
  else if( randR < ( R + T ) ) // Transmit
    {
      if( currentVolume == ePCVacuum || previousVolume == ePCVacuum || currentVolume == eOut ) // Others do not transmit
        Refract( normal, n1, real( n3 ), direction, polarisation );
      else
        weight = 0.0; // Just absorbed if not transmissive.
    }
  else // Absorbed
    {
      if( currentVolume == ePCVacuum
          || ( currentVolume == eGlass && previousVolume == ePCAlVacuum )
          || ( currentVolume == eGlass && previousVolume == ePCVacuum ) )
        {
          // Could cause a hit
          G4double signalProb = GetSignalProbability(ipmt, position, energy, weight );

          const double uniRand = G4UniformRand();
          if( uniRand < signalProb )
            {
              DS::MCPE photoelectron;
              photoelectron.SetPosition( TVector3( position.x(), position.y(), position.z() ) );
              photoelectron.SetCreationTime( time );
              photoelectron.SetCharge( 1 );
              photoelectron.SetPhotonTrackID( trackID );
              ExtractPhotonHistory(fastTrack,photoelectron);
              HitPMTCollection::Get()->DetectPhotoelectron( ipmt, photoelectron );
              RecordPhotoelectron( ipmt );
            }
        }
      weight = 0.0;
    }
  // Must do nothing here
  return false;
}

G4double
PMTOpticalModel::GetSignalProbability( const G4int ipmt,
                                       const G4ThreeVector,
                                       const G4double energy,
                                       const G4double weight )
{
  G4double signalProb = fPEEscapeProbability->Value( energy ) * fCollectionEfficiency;
  signalProb *= RAT::PhotonThinning::GetFactor() * weight;
  signalProb *= PMTVariation::Get()->GetVariation( ipmt );
  return signalProb;
}

void PMTOpticalModel::CalculateTRFresnel( const G4ThreeVector direction,
                                          const G4ThreeVector polarisation,
                                          const G4ThreeVector normal,
                                          const G4double n1,
                                          const std::complex<double> n2,
                                          G4double& T,
                                          G4double& R )
{
  double cos1 = direction.dot( normal );
  if( cos1 < 0.0 )
    cos1 = -cos1;
  const double sin1_2 = 1.0 - cos1 * cos1;
  const complex<double> n12 = n1 / n2;
  const complex<double> cos2 = sqrt( complex<double>(1.0,0.0) - ( n1 / n2 ) * ( n1 / n2 ) * sin1_2 );

  const complex<double> Ds = ( n2 * cos2 + n1 * cos1 );
  const complex<double> Dp = ( n2 * cos1 + n1 * cos2 );


  const double Rsc = abs( ( n1 * cos1 - n2 * cos2 ) / Ds );
  const double Rpc = abs( ( n1 * cos2 - n2 * cos1 ) / Dp );
  const double Rs = Rsc * Rsc;
  const double Rp = Rpc * Rpc;

  const complex<double> Nt = complex<double>( 4.0, 0.0 ) * n12 * cos1 * cos2;
  const double Ts = abs( n2 * n2 * Nt / ( Ds * Ds ) );
  const double Tp = abs( n2 * n2 * Nt / ( Dp * Dp ) );

  const double sFraction = GetSPolarisationFraction( normal, direction, polarisation );
  if( sFraction < 0.0 || sFraction > 1.0 )
    warn << "PMTOpticalModel::CalculateTRFresnel incorrect polarisation calculated as :" << sFraction << endl;
  T = sFraction * Ts + ( 1.0 - sFraction ) * Tp;
  R = sFraction * Rs + ( 1.0 - sFraction ) * Rp;
  if( R > 1.1 || R < 0.0 ) // T can be > 1.0
    warn << "PMTOpticalModel::CalculateTRFresnel R is too high " << Rs << " " << Rp << " " << R << endl;
}

void
PMTOpticalModel::CalculateTRThinFilm( const G4ThreeVector direction,
                                      const G4ThreeVector polarisation,
                                      const G4ThreeVector normal,
                                      const G4double lN1,
                                      const std::complex<double> n2,
                                      const std::complex<double> n3,
                                      const G4double energy,
                                      const G4double pcThickness,
                                      G4double& T,
                                      G4double& R )
{
  const double wavelength = twopi * hbarc / energy;

  const complex<double> n1( lN1, 0.0 );

  double cosT1 = direction.dot( normal );
  if( cosT1 < 0.0 )
    cosT1 = -cosT1;
  const double theta = acos( cosT1 );

  complex<double> cos1 = cos( theta );
  complex<double> sin1 = sin( theta );

  double eta = twopi * pcThickness / wavelength;
  complex<double> ratio13sin = (n1/n3) * (n1/n3) * sin1 * sin1;
  complex<double> cos3 = sqrt( complex<double>( 1.0, 0.0 ) - ratio13sin );
  complex<double> ratio12sin = (n1/n2) * (n1/n2) * sin1 * sin1;
  complex<double> cos2 = sqrt( complex<double>( 1.0, 0.0 ) - ratio12sin );
  double u = real( n2 * cos2 );
  double v = imag( n2 * cos2 );

  complex<double> r12;
  complex<double> r23;
  complex<double> t12;
  complex<double> t23;
  complex<double> g;

  double Ts, Rs;
  {
    // s Polarisation
    complex<double> n1c1 = n1 * cos1;
    complex<double> n2c2 = n2 * cos2;
    complex<double> n3c3 = n3 * cos3;
    r12 = ( n1c1 - n2c2 ) / ( n1c1 + n2c2 );
    r23 = ( n2c2 - n3c3 ) / ( n3c3 + n2c2 );
    t12 = complex<double>( 2.0, 0.0 ) * n1c1 / ( n1c1 + n2c2 );
    t23 = complex<double>( 2.0, 0.0 ) * n2c2 / ( n2c2 + n3c3 );
    g = n3c3 / n1c1;

    double abs_r12 = abs( r12 );
    double abs_r23 = abs( r23 );
    double abs_t12 = abs( t12 );
    double abs_t23 = abs( t23 );
    double arg_r12 = arg( r12 );
    double arg_r23 = arg( r23 );
    double exp1 = exp( 2.0 * v * eta);

    ToTR( u, eta, g, abs_r12, abs_r23, abs_t12, abs_t23, arg_r12, arg_r23, exp1, Ts, Rs );
  }
  double Tp, Rp;
  {
    // p Polarisation
    complex<double> n2c1 = n2 * cos1;
    complex<double> n3c2 = n3 * cos2;
    complex<double> n2c3 = n2 * cos3;
    complex<double> n1c2 = n1 * cos2;
    r12 = ( n2c1 - n1c2 ) / ( n2c1 + n1c2 );
    r23 = ( n3c2 - n2c3 ) / ( n3c2 + n2c3 );
    t12 = complex<double>( 2.0, 0.0 ) * n1 * cos1 / ( n2c1 + n1c2 );
    t23 = complex<double>( 2.0, 0.0 ) * n2 * cos2 / ( n3c2 + n2c3 );
    g = n3 * cos3 / ( n1 * cos1 );

    double abs_r12 = abs( r12 );
    double abs_r23 = abs( r23 );
    double abs_t12 = abs( t12 );
    double abs_t23 = abs( t23 );
    double arg_r12 = arg( r12 );
    double arg_r23 = arg( r23 );
    double exp1 = exp( 2.0 * v * eta);

    ToTR( u, eta, g, abs_r12, abs_r23, abs_t12, abs_t23, arg_r12, arg_r23, exp1, Tp, Rp );
  }
  const double sFraction = GetSPolarisationFraction( normal, direction, polarisation );

  T = sFraction * Ts + ( 1.0 - sFraction ) * Tp;
  R = sFraction * Rs + ( 1.0 - sFraction ) * Rp;
  if( T > 1.0 || R > 1.0 || T < 0.0 || R < 0.0 )
    warn << "PMTOpticalModel::CalculateTRThinFilm T and R are wrong, R:" << R << " T:" << T << endl;
}

void
PMTOpticalModel::ToTR( const double u,
                       const double eta,
                       const std::complex<double> g,
                       const double abs_r12,
                       const double abs_r23,
                       const double abs_t12,
                       const double abs_t23,
                       const double arg_r12,
                       const double arg_r23,
                       const double exp1,
                       double& T,
                       double& R )
{
  const double exp2 = 1/exp1;
  const double denom = exp1 + abs_r12 * abs_r12 * abs_r23 * abs_r23 * exp2 + 2.0 * abs_r12 * abs_r23 * cos( arg_r23 + arg_r12 + 2.0 * u * eta );
  R = abs_r12 * abs_r12 * exp1 + abs_r23 * abs_r23 *exp2 + 2.0 * abs_r12 * abs_r23 * cos( arg_r23 - arg_r12 + 2.0 * u * eta );
  R = R / denom;
  T = real( g ) * abs_t12 * abs_t12 * abs_t23 * abs_t23;
  T = T / denom;
}