G4StationFastCerenkov.cc

Go to the documentation of this file.

//
// ********************************************************************
// * DISCLAIMER                                                       *
// *                                                                  *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in:                      *
// *   http://cern.ch/geant4/license                                  *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.                                                             *
// *                                                                  *
// * This  code  implementation is the  intellectual property  of the *
// * GEANT4 collaboration.                                            *
// * By copying,  distributing  or modifying the Program (or any work *
// * based  on  the Program)  you indicate  your  acceptance of  this *
// * statement, and all its terms.                                    *
// ********************************************************************
//
// MODIFIED by P Skelton November 2003 - only use this as part of the
// G4FasterSimulatorModule in the Auger observatory simulation.
//
// Due to size, the code is split into several files. The subsidiary ones
// are #included into this file Subsidiary file: Propagate.icc

#include "G4StationFastCerenkov.h"
#include "G4StationConstruction.h"
#include "G4Utils.h"

#include <fwk/RunController.h>
#include <fwk/CentralConfig.h>
#include <utl/Reader.h>
#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/AugerException.h>
#include <utl/Math.h>

#include <G4PhysicalConstants.hh>
#include <G4SystemOfUnits.hh>

#include <algorithm>

using namespace fwk;
using namespace std;
using utl::Sqr;


#define IN_TANK 0
#define DOWN_IN_TANK 1
#define IN_DOME_1 11
#define IN_DOME_2 12
#define IN_DOME_3 13
#define IN_DOME_4 14
#define IN_INTERFACE_1 21
#define IN_INTERFACE_2 22
#define IN_INTERFACE_3 23
#define IN_INTERFACE_4 24
#define INWARD 31
#define OUTWARD 32

#define TARGET_PMT1 0
#define TARGET_PMT2 1
#define TARGET_PMT3 2
#define TARGET_SPMT 3
#define TARGET_WALL 4
#define TARGET_FLOOR 5
#define BOUNCED 10


namespace G4StationSimulatorOG {

  double G4StationFastCerenkov::fTankHeight;
  double G4StationFastCerenkov::fTankHalfHeight;
  double G4StationFastCerenkov::fTankThickness;
  double G4StationFastCerenkov::fTankRadius;
  double G4StationFastCerenkov::fTankRadiusSq;
  double G4StationFastCerenkov::fDomeRadius;
  double G4StationFastCerenkov::fDomeRadiusz;
  double G4StationFastCerenkov::fInterfaceRadius;
  double G4StationFastCerenkov::fInterfaceRadiusz;
  double G4StationFastCerenkov::fPMTRadius;
  double G4StationFastCerenkov::fPMTRadiusz;
  double G4StationFastCerenkov::fDomeRadiusSq;
  double G4StationFastCerenkov::fDomeRadiuszSq;
  double G4StationFastCerenkov::fInterfaceRadiusSq;
  double G4StationFastCerenkov::fInterfaceRadiuszSq;
  double G4StationFastCerenkov::fPMTRadiusSq;
  double G4StationFastCerenkov::fPMTRadiuszSq;
  double G4StationFastCerenkov::fHeightz;
  
  bool G4StationFastCerenkov::fHasSPMT;
  double G4StationFastCerenkov::fDomeRadius_SPMT;
  double G4StationFastCerenkov::fDomeThick_SPMT;
  double G4StationFastCerenkov::fInterfaceRadius_SPMT;
  double G4StationFastCerenkov::fInterfaceThick_SPMT;
  double G4StationFastCerenkov::fPMTRadius_SPMT;
  double G4StationFastCerenkov::fPMTActiveRadius_SPMT;
  double G4StationFastCerenkov::fPMTThick_SPMT;
  double G4StationFastCerenkov::fDomeRadius_SPMT_Sq;
  double G4StationFastCerenkov::fInterfaceRadius_SPMT_Sq;
  double G4StationFastCerenkov::fPMTRadius_SPMT_Sq;
  double G4StationFastCerenkov::fPMTActiveRadius_SPMT_Sq;

  double G4StationFastCerenkov::fSigmaAlpha;

  utl::TabulatedFunction G4StationFastCerenkov::fQE;
  utl::TabulatedFunction G4StationFastCerenkov::fQE_SPMT;

  utl::TabulatedFunction G4StationFastCerenkov::fWaterAbsLength;

  utl::TabulatedFunction G4StationFastCerenkov::fInterfaceRIndex;
  utl::TabulatedFunction G4StationFastCerenkov::fInterfaceAbsLength;

  utl::TabulatedFunction G4StationFastCerenkov::fDomeRIndex;
  utl::TabulatedFunction G4StationFastCerenkov::fDomeAbsLength;

  utl::TabulatedFunction G4StationFastCerenkov::fLinerReflectivity;
  utl::TabulatedFunction G4StationFastCerenkov::fLinerSpecularLobe;
  utl::TabulatedFunction G4StationFastCerenkov::fLinerSpecularSpike;
  utl::TabulatedFunction G4StationFastCerenkov::fLinerBackscatter;
  G4ThreeVector G4StationFastCerenkov::fPMTPos[5];

  double G4StationFastCerenkov::fRoofPos[5];


  G4StationFastCerenkov::G4StationFastCerenkov(const G4String& processName) :
    G4VContinuousProcess(processName),
    fG4StationSimulator(
      dynamic_cast<G4StationSimulator&>(RunController::GetInstance().GetModule("G4StationSimulatorOG"))
    )
  {
    if (verboseLevel > 0)
      G4cerr << GetProcessName() << " is created " << G4endl;
    BuildThePhysicsTable();

    // Rayleigh Scattering data taken from G4OpRayleigh.cc
    // isothermal compressibility of water
    const G4double betat = 7.658e-23 * m3 / MeV;
    // K Boltzman
    const G4double kboltz = 8.61739e-11 * MeV / kelvin;
    // Temperature of water is 10 degrees celsius
    const G4double temp = 283.15 * kelvin;
    const G4double hc4 = pow(h_Planck * c_light, 4); // constants defined
    // in G4PhysicalConstants
    fRayleighConst = 27 * hc4 / (8 * pow(utl::kPi, 3) * betat * temp * kboltz);
  }


  G4StationFastCerenkov::~G4StationFastCerenkov()
  {
    if (fPhysicsTable) {
      fPhysicsTable->clearAndDestroy();
      delete fPhysicsTable;
      fPhysicsTable = nullptr;
    }
  }


  void
  G4StationFastCerenkov::GetDataFromConstruction()
  {
    const sdet::Station& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();

    fTankRadius = G4StationConstruction::fgTankRadius;
    fTankHalfHeight = G4StationConstruction::fgTankHalfHeight;
    fTankHeight = fTankHalfHeight * 2;
    fTankThickness = G4StationConstruction::fgTankThickness;

    fPMTRadius = G4StationConstruction::fgPmtRmax;
    fPMTRadiusz = G4StationConstruction::fgPmtRzmax;

    fInterfaceRadius = G4StationConstruction::fgInterfaceRmax;
    fInterfaceRadiusz = G4StationConstruction::fgInterfaceRzmax;
    fInterfaceRIndex = G4StationConstruction::fgInterfaceRINDEX;
    fInterfaceAbsLength = G4StationConstruction::fgInterfaceABSORPTION;

    fDomeRadius = G4StationConstruction::fgDomeRmax;
    fDomeRadiusz = G4StationConstruction::fgDomeRzmax;
    fDomeRIndex = G4StationConstruction::fgPmtdomeRINDEX;
    fDomeAbsLength = G4StationConstruction::fgPmtdomeABSORPTION;

    fHeightz = G4StationConstruction::fgHeightz;

    fWaterAbsLength = G4StationConstruction::fgWaterABSORPTION;

    fLinerReflectivity = G4StationConstruction::fgLinerREFLECTIVITY;
    fLinerSpecularLobe = G4StationConstruction::fgLinerSPECULARLOBECONSTANT;
    fLinerSpecularSpike = G4StationConstruction::fgLinerSPECULARSPIKECONSTANT;
    fLinerBackscatter = G4StationConstruction::fgLinerBACKSCATTERCONSTANT;
    fSigmaAlpha = G4StationConstruction::fgSIGMA_ALPHA;

    // Fill Q.E. x COLL. EFF. of LPMTs 
    const auto& pmt = dStation.GetPMT(1);
    fQE = pmt.GetQuantumEfficiency();
    const double collectionEfficiency = pmt.GetCollectionEfficiency();
    for (auto& qe : fQE.YRange())
      qe *= collectionEfficiency;
    for (auto& en : fQE.XRange())
      en *= CLHEP::eV / utl::eV;

    const auto& cs = dStation.GetLocalCoordinateSystem();
    for (int i = 1; i <= 3; ++i) {
      fPMTPos[i] = ToG4Vector(dStation.GetPMT(i).GetPosition(), cs, CLHEP::meter/utl::meter);
      fRoofPos[i] = fTankHeight + fTankThickness - fPMTPos[i].z();
    }

    // checked above - PMTs are in expected positions
    fTankRadiusSq = fTankRadius * fTankRadius;
    fDomeRadiusSq = fDomeRadius * fDomeRadius;
    fInterfaceRadiusSq = fInterfaceRadius * fInterfaceRadius;
    fPMTRadiusSq = fPMTRadius * fPMTRadius;
    fDomeRadiuszSq = fDomeRadiusz * fDomeRadiusz;
    fInterfaceRadiuszSq = fInterfaceRadiusz * fInterfaceRadiusz;
    fPMTRadiuszSq = fPMTRadiusz * fPMTRadiusz;

    // Fill SmallPMT quantities
    fHasSPMT = dStation.HasSmallPMT();

    if (fHasSPMT) {
      fPMTActiveRadius_SPMT = G4StationConstruction::fPmtRmin_SPMT;
      fPMTRadius_SPMT = G4StationConstruction::fPmtRmax_SPMT;
      fPMTThick_SPMT = G4StationConstruction::fGlassThickness_SPMT;
      fInterfaceRadius_SPMT = G4StationConstruction::fInterfaceRadius_SPMT;
      fInterfaceThick_SPMT = G4StationConstruction::fInterfaceThickness_SPMT;
      fDomeRadius_SPMT = G4StationConstruction::fDomeRadius_SPMT;
      fDomeThick_SPMT = G4StationConstruction::fDomeThickness_SPMT;

      // Fill Q.E. x COLL. EFF. of SPMT
      const auto& small_pmt = dStation.GetPMT(4);
      fQE_SPMT = small_pmt.GetQuantumEfficiency();
      const double collectionEfficiency_SPMT = small_pmt.GetCollectionEfficiency();      
      for (auto& qe : fQE_SPMT.YRange())
        qe *= collectionEfficiency_SPMT;      
      for (auto& en : fQE_SPMT.XRange())
        en *= CLHEP::eV / utl::eV;

      fPMTPos[4] = ToG4Vector(dStation.GetPMT(4).GetPosition(), cs, CLHEP::meter/utl::meter);
      fRoofPos[4] = fTankHeight + fTankThickness - fPMTPos[4].z();
      
      fPMTRadius_SPMT_Sq = Sqr(fPMTRadius_SPMT);
      fPMTActiveRadius_SPMT_Sq = Sqr(fPMTActiveRadius_SPMT);
      fInterfaceRadius_SPMT_Sq = Sqr(fInterfaceRadius_SPMT);
      fDomeRadius_SPMT_Sq = Sqr(fDomeRadius_SPMT);
    }
  }


  void
  G4StationFastCerenkov::DumpPhysicsTable()
    const
  {
    for (int i = 0, n = fPhysicsTable->entries(); i < n; ++i)
      static_cast<G4PhysicsOrderedFreeVector*>((*fPhysicsTable)[i])->DumpValues();
  }


  G4VParticleChange*
  G4StationFastCerenkov::AlongStepDoIt(const G4Track& track, const G4Step& step)
  {
    aParticleChange.Initialize(track);
    const G4DynamicParticle* particle = track.GetDynamicParticle();
    const G4Material* const material = track.GetMaterial();
    G4StepPoint* const preStepPoint  = step.GetPreStepPoint();
    G4StepPoint* const postStepPoint = step.GetPostStepPoint();

    const G4ThreeVector x0 = preStepPoint->GetPosition();
    const G4ThreeVector p0 = step.GetDeltaPosition().unit();
    const G4double t0 = preStepPoint->GetGlobalTime();

    if (t0 > 1e6)
      return G4VContinuousProcess::AlongStepDoIt(track, step);

    G4ThreeVector primaryDirection = preStepPoint->GetMomentumDirection();

    const G4double cosAlpha = primaryDirection.z();
    const G4double sinAlpha = sqrt(1 - Sqr(cosAlpha));
    const G4double cosBeta = primaryDirection.x() / sinAlpha;
    const G4double sinBeta = primaryDirection.y() / sinAlpha;

    G4MaterialPropertiesTable* const materialPropertiesTable = material->GetMaterialPropertiesTable();

    if (!materialPropertiesTable)
      return G4VContinuousProcess::AlongStepDoIt(track, step);

    G4MaterialPropertyVector* const rIndex = materialPropertiesTable->GetProperty("RINDEX");

    if (!rIndex)
      return G4VContinuousProcess::AlongStepDoIt(track, step);

    G4double meanNumPhotons = GetAverageNumberOfPhotons(particle, material, rIndex);
    if (meanNumPhotons <= 0) {
      aParticleChange.SetNumberOfSecondaries(0);
      return G4VContinuousProcess::AlongStepDoIt(track, step);
    }

    meanNumPhotons *= step.GetStepLength();

    const G4int numPhotons = static_cast<G4int>(G4Poisson(meanNumPhotons));

    // we don't want to produce secondaries within Geant4, since they will
    // be handled within this dedicated Cerenkov routine for Auger.
    aParticleChange.SetNumberOfSecondaries(0);
    if (numPhotons <= 0)
      return G4VContinuousProcess::AlongStepDoIt(track, step);

    const G4double pMin = rIndex->GetMinLowEdgeEnergy();
    const G4double pMax = rIndex->GetMaxLowEdgeEnergy();
    const G4double nMax = rIndex->GetMaxValue();

    const G4double dp = pMax - pMin;

    const G4double betaInverse = particle->GetTotalEnergy() / particle->GetTotalMomentum();

    const G4double maxCos = betaInverse / nMax;
    const G4double maxSin2 = 1 - maxCos * maxCos;

    double minQEMomentum = fQE.XFront();
    double maxQEMomentum = fQE.XBack();
    if (fHasSPMT) {
      minQEMomentum = std::max(fQE.XFront(), fQE_SPMT.XFront());
      maxQEMomentum = std::min(fQE.XBack(), fQE_SPMT.XBack());
    }
    // leave in stuff for dealing with change of refractive index with
    // wavelength - even if the current setup ignores dispersion in water...

    for (G4int i = 0; i < numPhotons; ++i) {
      G4double rand = 0;
      G4double sin2Theta = 0;
      G4double cosTheta = 0;
      do {
        rand = G4UniformRand();
        fSampledMomentum = pMin + rand * dp;
        fSampledRI = rIndex->Value(fSampledMomentum);
        cosTheta = betaInverse / fSampledRI;
        sin2Theta = 1 - cosTheta * cosTheta;
        rand = G4UniformRand();
      } while (rand * maxSin2 > sin2Theta);

      // Kill according to quantum efficiency of PMT and collection efficiency

      if (fSampledMomentum < minQEMomentum || fSampledMomentum > maxQEMomentum)
        continue;

      fPhotonMaxQE = fQE.Y(fSampledMomentum);
      if (fHasSPMT)
        fPhotonMaxQE = std::max(fPhotonMaxQE, fQE_SPMT.Y(fSampledMomentum));

      if (G4UniformRand() > fPhotonMaxQE)
        continue;  // Photon is killed

      // calculate absorption and Rayleigh scattering constants for this photon
      int iFirstValidBin_Wat = 0;
      int iLastValidBin_Wat = fWaterAbsLength.GetNPoints() - 1;

      if (fSampledMomentum < fWaterAbsLength[iFirstValidBin_Wat].X() || fSampledMomentum > fWaterAbsLength[iLastValidBin_Wat].X())
        continue;

      const double absorptionMFP = fWaterAbsLength.Y(fSampledMomentum);  // CHECK

      const G4double sampledRI_Sq = fSampledRI * fSampledRI;
      const G4double rIConst = (sampledRI_Sq + 2) * (sampledRI_Sq - 1);

      // does this need to be member function?
      const double rayleighMFP = fRayleighConst / (pow(rIConst, 2) * pow(fSampledMomentum, 4));

      const G4double invSumMFP = 1 / (absorptionMFP + rayleighMFP);

      fInvMFP = 1 / absorptionMFP + 1 / rayleighMFP;
      fRayleighFrac = absorptionMFP * invSumMFP;
      fAbsorbFrac = rayleighMFP * invSumMFP;

      // calculate photon vectors (direction, polarisation)
      const G4double sinTheta = sqrt(sin2Theta);
      rand = G4UniformRand();

      const G4double phi = 2 * utl::kPi * rand;
      const G4double sinThetaSinPhi = sin(phi) * sinTheta;
      const G4double sinThetaCosPhi = cos(phi) * sinTheta;

      const G4double px = cosAlpha * cosBeta * sinThetaCosPhi - sinBeta * sinThetaSinPhi + cosBeta * sinAlpha * cosTheta;
      const G4double py = cosAlpha * sinBeta * sinThetaCosPhi + cosBeta * sinThetaSinPhi + sinBeta * sinAlpha * cosTheta;
      const G4double pz = cosAlpha * cosTheta - sinAlpha * sinThetaCosPhi;

      fPhotonMomentum.set(px, py, pz); // This is a momentum direction
      fPhotonPolarization = primaryDirection - cosTheta * fPhotonMomentum;
      fPhotonPolarization /= fPhotonPolarization.mag();

      rand = G4UniformRand();
      const G4double delta = rand * step.GetStepLength();
      const G4double deltaTime = delta / (0.5 * (preStepPoint->GetVelocity() + postStepPoint->GetVelocity()));
      fPhotonTime = t0 + deltaTime;
      fPhotonPosition = x0 + rand * step.GetDeltaPosition();

      const G4double radiusTest = fPhotonPosition.x() * fPhotonPosition.x() + fPhotonPosition.y() * fPhotonPosition.y() - fTankRadiusSq;

      // if Cerenkov photon is inside tank, then track it
      if (radiusTest < 0 &&
          fPhotonPosition.z() > fTankThickness &&
          fPhotonPosition.z() < fTankHeight + fTankThickness)
        TrackPhotonInTank();
    }

    if (verboseLevel > 0) {
      G4cerr << "\n Exiting from G4Cerenkov::DoIt -- NumberOfSecondaries = "
             << aParticleChange.GetNumberOfSecondaries() << G4endl;
    }

    return G4VContinuousProcess::AlongStepDoIt(track, step);
  }


  void
  G4StationFastCerenkov::BuildThePhysicsTable()
  {
    if (fPhysicsTable)
      return;

    const G4MaterialTable* const materialTable = G4Material::GetMaterialTable();
    const G4int numOfMaterials = G4Material::GetNumberOfMaterials();
    fPhysicsTable = new G4PhysicsTable(numOfMaterials);

    for (G4int i = 0; i < numOfMaterials; ++i) {

      G4PhysicsOrderedFreeVector* const physicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();
      G4Material* const material = (*materialTable)[i];
      G4MaterialPropertiesTable* const materialPropertiesTable = material->GetMaterialPropertiesTable();

      if (materialPropertiesTable) {

        G4MaterialPropertyVector* const rIndexVector = materialPropertiesTable->GetProperty("RINDEX");

        if (rIndexVector) {

          G4double currentRI = (*rIndexVector)[0];

          if (currentRI > 1) {

            G4double currentPM = rIndexVector->Energy(0);

            G4double currentCAI = 0;
            physicsOrderedFreeVector->InsertValues(currentPM , currentCAI);
            G4double prevRI = currentRI;
            G4double prevPM = currentPM;
            G4double prevCAI = currentCAI;

            for (size_t ii = 1; ii < rIndexVector->GetVectorLength(); ++ii) {

              currentRI = (*rIndexVector)[ii];
              currentPM = rIndexVector->Energy(ii);
              currentCAI = 0.5 * (1 / (prevRI * prevRI) + 1 / (currentRI * currentRI));
              currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
              physicsOrderedFreeVector->InsertValues(currentPM, currentCAI);
              prevPM = currentPM;
              prevCAI = currentCAI;
              prevRI = currentRI;

            }

          }

        }

      }

      fPhysicsTable->insertAt(i, physicsOrderedFreeVector);

    }
  }


  G4double
  G4StationFastCerenkov::GetContinuousStepLimit(const G4Track& track, G4double, G4double, G4double&)
  {
    if (fMaxPhotons <= 0)
      return DBL_MAX;

    const G4Material* const material = track.GetMaterial();
    G4MaterialPropertiesTable* const materialPropertiesTable = material->GetMaterialPropertiesTable();

    if (!materialPropertiesTable)
      return DBL_MAX;

    G4MaterialPropertyVector* const rIndex = materialPropertiesTable->GetProperty("RINDEX");

    if (!rIndex)
      return DBL_MAX;

    const G4DynamicParticle* const particle = track.GetDynamicParticle();
    const G4double meanNumPhotons = GetAverageNumberOfPhotons(particle, material, rIndex);
    if (meanNumPhotons <= 0)
      return DBL_MAX;

    const G4double stepLimit = fMaxPhotons / meanNumPhotons;
    return stepLimit;
  }


  G4double
  G4StationFastCerenkov::GetAverageNumberOfPhotons(const G4DynamicParticle* const particle,
                                                const G4Material* const material,
                                                G4MaterialPropertyVector* const rIndex)
    const
  {
    const G4double rFact = 369.81 / (eV * cm);

    if (particle->GetTotalMomentum() <= 0)
      return 0;

    const G4double betaInverse = particle->GetTotalEnergy() / particle->GetTotalMomentum();

    const G4int materialIndex = material->GetIndex();

    G4PhysicsOrderedFreeVector* const cerenkovAngleIntegrals =
      (G4PhysicsOrderedFreeVector*)((*fPhysicsTable)(materialIndex));

    if (!cerenkovAngleIntegrals->IsFilledVectorExist())
      return 0;

    G4double pMin = rIndex->GetMinLowEdgeEnergy();
    const G4double pMax = rIndex->GetMaxLowEdgeEnergy();
    const G4double nMin = rIndex->GetMinValue();
    const G4double nMax = rIndex->GetMaxValue();

    const G4double cAImax = cerenkovAngleIntegrals->GetMaxValue();
    G4double dp = 0;
    G4double ge = 0;

    if (nMax < betaInverse) {

      dp = 0;
      ge = 0;

    } else if (nMin > betaInverse) {

      dp = pMax - pMin;
      ge = cAImax;

    } else {

      pMin = rIndex->GetEnergy(betaInverse);
      dp = pMax - pMin;
      G4bool isOutRange = false;
      const G4double cAImin = cerenkovAngleIntegrals->GetValue(pMin, isOutRange);
      ge = cAImax - cAImin;

      if (verboseLevel > 0) {
        G4cerr << "CAImin = " << cAImin << G4endl;
        G4cerr << "ge = " << ge << G4endl;
      }

    }

    const G4double charge = particle->GetDefinition()->GetPDGCharge();
    const G4double numPhotons = rFact * pow(charge / eplus, 2) * (dp - ge * betaInverse*betaInverse);

    return numPhotons;
  }


  //#include "Propagate.icc" // Contains the code for propogating photons -
  // part of the same class and namespace, but put in a different file
  // for clarity

  void
  G4StationFastCerenkov::DoRayleighScatter()
  {
    const G4double rand = 2 * G4UniformRand() - 1;
    const G4double cosTheta = (rand < 0) ?  -pow(-rand, 1 / 3.) : pow(rand, 1 / 3.);
    const G4double sinTheta = sqrt(1 - cosTheta * cosTheta);
    const G4double phi = 2 * utl::kPi * G4UniformRand();
    const G4double cosPhi = cos(phi);
    const G4double sinPhi = sin(phi);
    const G4double a = fPhotonPolarization.x();
    const G4double b = fPhotonPolarization.y();
    const G4double c = fPhotonPolarization.z();
    const G4double k = 1 / sqrt(1 - c * c);
    const G4double R11 = k * b;
    const G4double R12 = k * a * c;
    const G4double R21 = -k * a;
    const G4double R22 = k * b * c;
    const G4double R32 = -1 / k;

    G4double x = 0;
    G4double y = 0;
    G4double z = 0;

    if (G4UniformRand() < 0.5) {

      x = sinTheta * cosPhi;
      y = sinTheta * sinPhi;
      z = cosTheta;
      fPhotonPolarization.set((R11 * x + R12 * y + a * z),
                              (R21 * x + R22 * y + b * z),
                              (R32 * y + c * z));
      x = cosTheta * cosPhi;
      y = cosTheta * sinPhi;
      z = -sinTheta;
      fPhotonMomentum.set((R11 * x + R12 * y + a * z),
                          (R21 * x + R22 * y + b * z),
                          (R32 * y + c * z));

    } else {

      x = sinTheta * cosPhi;
      y = sinTheta * sinPhi;
      z = cosTheta;
      fPhotonPolarization.set((R11 * x + R12 * y + a * z),
                              (R21 * x + R22 * y + b * z),
                              (R32 * y + c * z));
      x = -cosTheta * cosPhi;
      y = -cosTheta * sinPhi;
      z = sinTheta;

      fPhotonMomentum.set((R11 * x + R12 * y + a * z),
                          (R21 * x + R22 * y + b * z),
                          (R32 * y + c * z));

    }
  }


  // included into the G4StationFastCerenkov class for handling cerenkov photons
  // in the Geant4 Auger tank simulation

  void
  G4StationFastCerenkov::TrackPhotonInTank()
  {
    // This routine makes use of knowledge about the design of the Auger
    // tanks that it shouldn't have (numbers not obtained through official
    // routes, but hard-coded). That's life, if you want to make a program
    // run faster.

    // The actual tracking is done in PropagateInTank and other routines -
    // this one just determines which volume we start in, and gets us to
    // the correct starting point in the nested methods. Once we return
    // to this routine we fall straight back down to the level below,
    // and move on to the next photon. Any detection of the photon has to
    // done in the routines called from this one.

    if (fPhotonMomentum.z() < 0 &&
        fPhotonPosition.z() < fTankHeight + fTankThickness - fDomeRadiusz) {
      PropagateInTank(DOWN_IN_TANK);
      return;
    }

    G4double closest = 0;
    G4ThreeVector displacement;

    // loop over the PMTs and return whenever a PMT is hit
    for (unsigned int i = 1; i <= 3; ++i) {
      displacement = fPhotonPosition - fPMTPos[i];
      closest = (Sqr(displacement.x()) + Sqr(displacement.y())) / fDomeRadiusSq + Sqr(displacement.z()) / fDomeRadiuszSq;
      if (closest < 1) {
        closest = (Sqr(displacement.x()) + Sqr(displacement.y())) / fInterfaceRadiusSq + Sqr(displacement.z()) / fInterfaceRadiuszSq;
        if (closest < 1) {
          closest = (Sqr(displacement.x()) + Sqr(displacement.y())) / fPMTRadiusSq + Sqr(displacement.z()) / fPMTRadiuszSq;
          if (closest < 1 && -displacement.z() > fHeightz) {  // z < fHeightz : insensitive area
            fG4StationSimulator.AddPhoton(i, fPhotonTime);
            return;
          }
          PropagateInTank(IN_INTERFACE_1 + i - 1);
          return;
        }
        PropagateInTank(IN_DOME_1 + i - 1);
        return;
      }
    }

    PropagateInTank(IN_TANK);
  }


  void
  G4StationFastCerenkov::PropagateInTank(G4int flag)
  {
    G4bool dead;
    if (flag > DOWN_IN_TANK) {
      if (flag == IN_DOME_1 || flag == IN_INTERFACE_1) {
        dead = TransitionToDome<1>(flag);
        if (dead)
          return;
      } else if (flag == IN_DOME_2 || flag == IN_INTERFACE_2) {
        dead = TransitionToDome<2>(flag);
        if (dead)
          return;
      } else if (flag == IN_DOME_3 || flag == IN_INTERFACE_3) {
        dead = TransitionToDome<3>(flag);
        if (dead)
          return;
      } else { // IN_DOME_4 or IN_INTERFACE_4
        dead = TransitionToDome<4>(flag);
        if (dead)
          return;
      }
      if (fPhotonMomentum.z() < -0.0463)
        flag = DOWN_IN_TANK;
      // A photon emerging anywhere from a dome with z component of
      // direction below this cannot possibly hit another dome
    }

    G4double distanceToFloor, distanceToWall, distanceTravelled, distance;
    G4double distanceToRoof, distanceToPMT1, distanceToPMT2, distanceToPMT3, distanceToSPMT;
    G4double alpha, beta, gamma;
    G4double testVal;
    G4bool repeat;
    G4double closest1, closest2, closest3, closest4;
    G4ThreeVector closest, displacement1, displacement2, displacement3, displacement4;
    G4double dotProduct, rand, distanceToPMT;
    G4int hitPMT, targetPMT = -1, target;

    do {
      repeat = false;
      if (flag == DOWN_IN_TANK) {
        // photon cannot hit domes, so only check for base and wall

        if (!fPhotonMomentum.z()) {
          WARNING("The z component of the momentul is 0, skipping photon.");
          return;
        }

        distanceToFloor = -(fPhotonPosition.z() - fTankThickness) / fPhotonMomentum.z();

        gamma = 1 / (Sqr(fPhotonMomentum.x()) + Sqr(fPhotonMomentum.y()));
        alpha = (fPhotonMomentum.x()*fPhotonPosition.x() + fPhotonMomentum.y()*fPhotonPosition.y()) * gamma;
        beta = (Sqr(fPhotonPosition.x()) + Sqr(fPhotonPosition.y()) - fTankRadiusSq) * gamma;

        const double distance = Sqr(alpha) - beta;
        if (distance < 0) {
          WARNING("Distance to wall is not finite, skipping photon");
          return;
        }

        distanceToWall = sqrt(distance) - alpha;

        // check for Rayleigh scattering or absorption en route
        rand = G4UniformRand();
        testVal = 1 - exp(-min(distanceToWall, distanceToFloor) * fInvMFP);

        if (rand < testVal) {
          // some interaction has occurred before reaching destination
          if (rand < testVal * fRayleighFrac) {
            // we have a Rayleigh scatter
            distanceTravelled = -log(1 - rand/fRayleighFrac) / fInvMFP;
            fPhotonPosition += distanceTravelled * fPhotonMomentum;
            fPhotonTime += distanceTravelled * fSampledRI / c_light;
            DoRayleighScatter();  // to pick new polarisation and direction
            if (fPhotonMomentum.z() < 0 &&
                fPhotonPosition.z() < fTankHeight + fTankThickness - fPMTRadiusz)
              flag = DOWN_IN_TANK;
            else
              flag = IN_TANK;
            repeat = true;
          } else
            return;  // we have absorption
        }
        if (!repeat) {
          // no absorption or scattering - we reach the destination surface
          if (distanceToWall < distanceToFloor) {
            // we hit the wall
            fPhotonPosition += distanceToWall * fPhotonMomentum;
            fPhotonTime += distanceToWall * fSampledRI / c_light;
            dead = ScatterOffWall();
          } else {
            // we hit the floor
            fPhotonPosition += distanceToFloor * fPhotonMomentum;
            fPhotonTime += distanceToFloor * fSampledRI / c_light;
            dead = ScatterOffFloor();
          }
          if (dead)
            return;
          repeat = true;
        }
      } else {
     
        hitPMT = 0;
        distanceToRoof =
          (fPhotonMomentum.z() > 0) ?
          (fTankHeight + fTankThickness - fPhotonPosition.z()) / fPhotonMomentum.z() :
          -(fPhotonPosition.z() - fTankThickness) / fPhotonMomentum.z();  // floor

        gamma = 1 / (Sqr(fPhotonMomentum.x()) + Sqr(fPhotonMomentum.y()));
        alpha = (fPhotonMomentum.x()*fPhotonPosition.x() + fPhotonMomentum.y()*fPhotonPosition.y()) * gamma;
        beta = (Sqr(fPhotonPosition.x()) + Sqr(fPhotonPosition.y()) - fTankRadiusSq) * gamma;
        distanceToWall = sqrt(Sqr(alpha) - beta) - alpha;

        displacement1 = fPhotonPosition - fPMTPos[1];
        dotProduct = displacement1.dot(fPhotonMomentum);
        closest = displacement1 - dotProduct * fPhotonMomentum;
        closest1 = (Sqr(closest.x()) + Sqr(closest.y())) / fDomeRadiusSq + Sqr(closest.z()) / fDomeRadiuszSq;
        if (closest1 < 1)
          hitPMT += 1;

        displacement2 = fPhotonPosition - fPMTPos[2];
        dotProduct = displacement2.dot(fPhotonMomentum);
        closest = displacement2 - dotProduct * fPhotonMomentum;
        closest2 = (Sqr(closest.x()) + Sqr(closest.y())) / fDomeRadiusSq + Sqr(closest.z()) / fDomeRadiuszSq;
        if (closest2 < 1)
          hitPMT += 2;

        displacement3 = fPhotonPosition - fPMTPos[3];
        dotProduct = displacement3.dot(fPhotonMomentum);
        closest = displacement3 - dotProduct * fPhotonMomentum;
        closest3 = (Sqr(closest.x()) + Sqr(closest.y())) / fDomeRadiusSq + Sqr(closest.z()) / fDomeRadiuszSq;
        if (closest3 < 1)
          hitPMT += 4;

        if (fHasSPMT && fPhotonMomentum.z() > 0) {
          displacement4 = fPhotonPosition + distanceToRoof * fPhotonMomentum - fPMTPos[4];
          closest4 = (Sqr(displacement4.x()) + Sqr(displacement4.y())) / fDomeRadius_SPMT_Sq;
          if (closest4 < 1)
            hitPMT += 8;
        }

        switch (hitPMT) {
        case 0:
          distanceToPMT = 1e99;
          break;
        case 1:
          distanceToPMT = GetEllipsoidIntersect(displacement1, fDomeRadius, fDomeRadiusz);
          targetPMT = TARGET_PMT1;
          break;
        case 2:
          distanceToPMT = GetEllipsoidIntersect(displacement2, fDomeRadius, fDomeRadiusz);
          targetPMT = TARGET_PMT2;
          break;
        case 3:
          distanceToPMT1 = GetEllipsoidIntersect(displacement1, fDomeRadius, fDomeRadiusz);
          distanceToPMT2 = GetEllipsoidIntersect(displacement2, fDomeRadius, fDomeRadiusz);
          if (distanceToPMT1 < distanceToPMT2) {
            targetPMT = TARGET_PMT1;
            distanceToPMT = distanceToPMT1;
          } else {
            targetPMT = TARGET_PMT2;
            distanceToPMT = distanceToPMT2;
          }
          break;
        case 4:
          distanceToPMT = GetEllipsoidIntersect(displacement3, fDomeRadius, fDomeRadiusz);
          targetPMT = TARGET_PMT3;
          break;
        case 5:
          distanceToPMT1 = GetEllipsoidIntersect(displacement1, fDomeRadius, fDomeRadiusz);
          distanceToPMT3 = GetEllipsoidIntersect(displacement3, fDomeRadius, fDomeRadiusz);
          if (distanceToPMT1 < distanceToPMT3) {
            targetPMT = TARGET_PMT1;
            distanceToPMT = distanceToPMT1;
          } else {
            targetPMT = TARGET_PMT3;
            distanceToPMT = distanceToPMT3;
          }
          break;
        case 6:
          distanceToPMT2 = GetEllipsoidIntersect(displacement2, fDomeRadius, fDomeRadiusz);
          distanceToPMT3 = GetEllipsoidIntersect(displacement3, fDomeRadius, fDomeRadiusz);
          if (distanceToPMT2 < distanceToPMT3) {
            targetPMT = TARGET_PMT2;
            distanceToPMT = distanceToPMT2;
          } else {
            targetPMT = TARGET_PMT3;
            distanceToPMT = distanceToPMT3;
          }
          break;
        case 8:
          distanceToPMT = distanceToRoof;
          targetPMT = TARGET_SPMT;
          break;
        case 9:
          distanceToPMT1 = GetEllipsoidIntersect(displacement1, fDomeRadius, fDomeRadiusz);
          distanceToSPMT = distanceToRoof;
          if (distanceToPMT1 < distanceToSPMT) {
            targetPMT = TARGET_PMT1;
            distanceToPMT = distanceToPMT1;
          } else {
            targetPMT = TARGET_SPMT;
            distanceToPMT = distanceToSPMT;
          }
          break;
        case 10:
          distanceToPMT2 = GetEllipsoidIntersect(displacement2, fDomeRadius, fDomeRadiusz);
          distanceToSPMT = distanceToRoof;
          if (distanceToPMT2 < distanceToSPMT) {
            targetPMT = TARGET_PMT2;
            distanceToPMT = distanceToPMT2;
          } else {
            targetPMT = TARGET_SPMT;
            distanceToPMT = distanceToSPMT;
          }
          break;
        case 12:
          distanceToPMT3 = GetEllipsoidIntersect(displacement3, fDomeRadius, fDomeRadiusz);
          distanceToSPMT = distanceToRoof;
          if (distanceToPMT3 < distanceToSPMT) {
            targetPMT = TARGET_PMT3;
            distanceToPMT = distanceToPMT3;
          } else {
            targetPMT = TARGET_SPMT;
            distanceToPMT = distanceToSPMT;
          }
          break;
        default:
          // hit 3 or more pmts - error
          WARNING("Error determining next photon intersection in PropagateInTank()");
          return;
        }

        if (distanceToPMT < distanceToWall) {
          if (distanceToPMT < distanceToRoof || 
              targetPMT == TARGET_SPMT) {
            target = targetPMT;
            distance = distanceToPMT;
          } else {
            target = TARGET_FLOOR;
            distance = distanceToRoof;
          }
        } else {
          if (distanceToRoof < distanceToWall) {
            target = TARGET_FLOOR;
            distance = distanceToRoof;
          } else {
            target = TARGET_WALL;
            distance = distanceToWall;
          }
        }
        // check for Rayleigh scattering or absorption en route
        rand = G4UniformRand();
        testVal = 1 - exp(-distance * fInvMFP);
        if (rand < testVal) {
          // some interaction has occurred before reaching destination
          if (rand >= testVal * fRayleighFrac)
            return;  // absorption
          // if we get here we have a Rayleigh scatter
          distanceTravelled = -log(1 - rand/fRayleighFrac) / fInvMFP;
          fPhotonPosition += distanceTravelled * fPhotonMomentum;
          fPhotonTime += distanceTravelled * fSampledRI / c_light;
          DoRayleighScatter(); // to pick new polarisation and direction
          if (fPhotonMomentum.z() < 0 &&
              fPhotonPosition.z() < fTankHeight + fTankThickness - fPMTRadiusz)
            flag = DOWN_IN_TANK;
          else
            flag = IN_TANK;
          repeat = true;
        } else {
          fPhotonPosition += distance * fPhotonMomentum;
          fPhotonTime += distance * fSampledRI / c_light;
          switch (target) {
          case TARGET_WALL:
            dead = ScatterOffWall();
            if (dead)
              return;
            break;
          case TARGET_FLOOR:
            if (fPhotonMomentum.z() < 0)
              dead = ScatterOffFloor();
            else
              dead = ScatterOffRoof();
            if (dead)
              return;
            break;
          case TARGET_PMT1:
            dead = TransitionToDome<1>(INWARD);
            if (dead)
              return;
            break;
          case TARGET_PMT2:
            dead = TransitionToDome<2>(INWARD);
            if (dead)
              return;
            break;
          case TARGET_PMT3:
            dead = TransitionToDome<3>(INWARD);
            if (dead)
              return;
            break;
          case TARGET_SPMT:
            dead = TransitionToDome<4>(INWARD);
            if (dead)
              return;
            break;
          default:
            WARNING("Error finding target in PropagateInTank");
            return;
          }
          repeat = true;
        }
      }

      // quick check if we can ignore PMTs on next track
      if (fPhotonPosition.z() < fTankHeight + fTankThickness - fDomeRadiusz &&
          fPhotonMomentum.z() < 0)
        flag = DOWN_IN_TANK;
      else
        flag = IN_TANK;

    } while (repeat);

    // never get here - we either get detected or absorbed inside the loop
    const std::string err = "Tank photon tracking ran into an impossible condition";
    ERROR(err);
    throw utl::OutOfBoundException(err);
  }


  // PMT ROUTINES

  template<int pmtId>
  G4bool
  G4StationFastCerenkov::TransitionToDome(G4int flag)
  {
    fPhotonPosition -= fPMTPos[pmtId];
    G4ThreeVector alpha, beta, delta, normal;
    G4bool reflected, dead;
    G4double n, n2, transmission, vDotN;
    G4double paraComp, perpComp, paraMag, perpMag, x;
    G4ThreeVector paraPolarization, perpPolarization;

    if (flag == IN_INTERFACE_1 ||
        flag == IN_DOME_1 ||
        flag == IN_INTERFACE_2 ||
        flag == IN_DOME_2 ||
        flag == IN_INTERFACE_3 ||
        flag == IN_DOME_3 ||
        flag == IN_INTERFACE_4 ||
        flag == IN_DOME_4 ) {
      dead = PropagateInDome<pmtId>(flag); // flag needs to be altered in PiD
      if (dead)
        return true;
    }

    do {
      if (flag == INWARD) {
        n = fSampledRI;
        n2 = fDomeRIndex.Y(fSampledMomentum);
      } else {
        n = fDomeRIndex.Y(fSampledMomentum);
        n2 = fSampledRI;
      }
      normal.set(
        2 * fPhotonPosition.x() / fDomeRadiusSq,
        2 * fPhotonPosition.y() / fDomeRadiusSq,
        2 * fPhotonPosition.z() / fDomeRadiuszSq
      );
      normal *= 1/normal.mag();
      if (pmtId == 4)
        normal.set(0, 0, -1);
      vDotN = -normal.dot(fPhotonMomentum);  // is actually cos(theta)
      if (vDotN > 0.999999) {
        transmission = 4 * n * n2 / pow(n + n2, 2);
        reflected = (G4UniformRand() > transmission);
        if (reflected) {
          fPhotonMomentum = -fPhotonMomentum;
          fPhotonPolarization = -fPhotonPolarization;
          if (flag == INWARD) {
            fPhotonPosition += fPMTPos[pmtId];
            return false;
          }
          flag = INWARD;
          dead = PropagateInDome<pmtId>(flag);
          if (dead)
            return true;
        } else {
          if (flag == OUTWARD) {
            fPhotonPosition += fPMTPos[pmtId];
            return false;
          }
          dead = PropagateInDome<pmtId>(flag);
          if (dead)
            return true;
        }
      } else {
        alpha = normal.cross(fPhotonMomentum);
        alpha *= 1/alpha.mag();
        x = Sqr(n2) - Sqr(n) * (1 - Sqr(vDotN));
        if (x <= 0)
          x = 0;  // total reflection
        else
          x = (vDotN > 0) ? -sqrt(x) : sqrt(x);
        paraMag = fPhotonPolarization.dot(alpha);
        beta = fPhotonMomentum.cross(alpha); // should be unit vector
        perpMag = fPhotonPolarization.dot(beta);
        perpComp = perpMag / (n * vDotN - x);
        paraComp = paraMag / (n2 * vDotN - (n * x / n2));
        transmission = -4 * n * vDotN * x * (Sqr(perpComp) + Sqr(paraComp));
        reflected = (G4UniformRand() >= transmission);
        if (x == 0)
          reflected = true;
        if (flag == INWARD) {
          if (reflected) {
            // reflected back into tank
            fPhotonMomentum += 2 * vDotN * normal;
            fPhotonMomentum *= 1/fPhotonMomentum.mag();
            beta = fPhotonMomentum.cross(alpha);
            fPhotonPolarization = paraComp * alpha + perpComp * beta;
            fPhotonPolarization *= 1/fPhotonPolarization.mag();
            fPhotonPosition *= 1.00001; // Fudge to avoid rounding errors
            fPhotonPosition += fPMTPos[pmtId];
            return false; // so we continue tracking in Tank
          } else {
            // transmitted in to dome
            fPhotonMomentum *= n;
            fPhotonMomentum += (n * vDotN + x) * normal;
            fPhotonMomentum *= 1/fPhotonMomentum.mag();
            beta = fPhotonMomentum.cross(alpha);
            fPhotonPolarization = paraComp * alpha + perpComp * beta;
            fPhotonPolarization *= 1/fPhotonPolarization.mag();
            dead = PropagateInDome<pmtId>(flag); // flag needs to be altered in PiD
            if (dead)
              return true;
          }
        } else {
          // reflected off inside of dome - remains in dome
          if (reflected) {
            fPhotonMomentum += 2 * vDotN * normal;
            fPhotonMomentum *= 1/fPhotonMomentum.mag();
            beta = fPhotonMomentum.cross(alpha);
            fPhotonPolarization = paraComp * alpha + perpComp * beta;
            fPhotonPolarization *= 1/fPhotonPolarization.mag();
            flag = INWARD;
            dead = PropagateInDome<pmtId>(flag);
            if (dead)
              return true;
          } else {
            // transmitted from dome back into tank
            fPhotonMomentum *= n;
            fPhotonMomentum += (n * vDotN + x) * normal;
            fPhotonMomentum *= 1/fPhotonMomentum.mag();
            beta = fPhotonMomentum.cross(alpha);
            fPhotonPolarization = paraComp * alpha + perpComp * beta;
            fPhotonPolarization *= 1/fPhotonPolarization.mag();
            fPhotonPosition *= 1.00001; // Fudge to avoid rounding errors
            fPhotonPosition += fPMTPos[pmtId];
            return false;
          }
        }
      }
    } while (true);
  }


  template<int pmtId>
  G4bool
  G4StationFastCerenkov::PropagateInDome(G4int& flag)
  {
    G4bool dead;
    if (flag == IN_INTERFACE_1 ||
        flag == IN_INTERFACE_2 ||
        flag == IN_INTERFACE_3 ||
        flag == IN_INTERFACE_4) {
      dead = TransitionToInterface<pmtId>(flag);
      if (dead)
        return true;
    }

    while (pmtId == 4) {
      // We are inside the SPMT dome, that is above the water level
      // So the z-component of the photon position relative to the tank roof is > 0
      G4double distanceToWater = (fRoofPos[pmtId] - fPhotonPosition.z()) / fPhotonMomentum.z();
      G4double distanceToInterface = (fRoofPos[pmtId] + fDomeThick_SPMT - fPhotonPosition.z()) / fPhotonMomentum.z();
 
      G4ThreeVector displacement;
      G4double distance;

      if (flag == INWARD) {

        if (distanceToInterface < 0) {
          const std::string err = "Photon tracking inside SPMT dome run into an impossible condition : INWARD & distToInterface<0";
          ERROR(err);
          throw utl::OutOfBoundException(err);
        }

        displacement = fPhotonPosition + distanceToInterface * fPhotonMomentum;
        distance = (Sqr(displacement.x()) + Sqr(displacement.y())) / fInterfaceRadius_SPMT_Sq;
        if (distance < 1) {
          dead = (G4UniformRand() > exp(-distanceToInterface / fDomeAbsLength.Y(fSampledMomentum)));
          if (dead)
            return true;
          fPhotonPosition += distanceToInterface * fPhotonMomentum;
          fPhotonTime += distanceToInterface * fDomeRIndex.Y(fSampledMomentum) / c_light;
          dead = TransitionToInterface<pmtId>(flag);
          if (dead)
            return true;
        } else 
          return true;  // photon killed

      } else {

        if (distanceToWater < 0) {
          const std::string err = "Photon tracking inside SPMT dome run into an impossible condition : OUTWARD & distToWater<0";
          ERROR(err);
          throw utl::OutOfBoundException(err);
        }

        displacement = fPhotonPosition + distanceToWater * fPhotonMomentum;
        distance = (Sqr(displacement.x()) + Sqr(displacement.y())) / fDomeRadius_SPMT_Sq;
        if (distance < 1) {
          dead = (G4UniformRand() > exp(-distanceToWater / fDomeAbsLength.Y(fSampledMomentum)));
          if (dead)
            return true;
          fPhotonPosition += distanceToWater * fPhotonMomentum;
          fPhotonTime += distanceToWater * fDomeRIndex.Y(fSampledMomentum) / c_light;
          return false;
        } else
          return true;  // photon killed
      }

    }

    G4double distanceThis, distanceOther, distanceToRoof;
    // we start at one dome surface, so only one non-zero solution exists
    do {
      distanceToRoof = (fRoofPos[pmtId] - fPhotonPosition.z()) / fPhotonMomentum.z();
      if (distanceToRoof < 0)
        distanceToRoof = 1e99;
      if (flag == INWARD) {
        distanceThis = GetEllipsoidIntersect(fPhotonPosition, fDomeRadius, fDomeRadiusz);
        distanceOther = GetEllipsoidIntersect(fPhotonPosition, fInterfaceRadius, fInterfaceRadiusz);
        if (distanceThis < distanceOther) {
          if (distanceToRoof < distanceThis)
            return true;  // photon killed
          dead = (G4UniformRand() > exp(-distanceThis / fDomeAbsLength.Y(fSampledMomentum)));
          if (dead)
            return true;
          fPhotonPosition += distanceThis * fPhotonMomentum;
          fPhotonTime += distanceThis * fDomeRIndex.Y(fSampledMomentum) / c_light;
          flag = OUTWARD;
          return false;
        } else {
          if (distanceToRoof < distanceOther)
            return true;  // photon killed
          dead = (G4UniformRand() > exp(-distanceOther / fDomeAbsLength.Y(fSampledMomentum)));
          if (dead)
            return true;
          fPhotonPosition += distanceOther * fPhotonMomentum;
          fPhotonTime += distanceOther * fDomeRIndex.Y(fSampledMomentum) / c_light;
          dead = TransitionToInterface<pmtId>(flag);
          if (dead)
            return true;
        }
      } else {
        distanceOther = GetEllipsoidIntersect(fPhotonPosition, fDomeRadius, fDomeRadiusz);
        if (distanceToRoof < distanceOther)
          return true; // photon killed
        dead = (G4UniformRand() > exp(-distanceOther / fDomeAbsLength.Y(fSampledMomentum)));
        if (dead)
          return true;
        fPhotonPosition += distanceOther * fPhotonMomentum;
        fPhotonTime += distanceOther * fDomeRIndex.Y(fSampledMomentum) / c_light;
        flag = OUTWARD;
        return false;
      }
    } while (true);
  }


  template<int pmtId>
  G4bool
  G4StationFastCerenkov::TransitionToInterface(G4int& flag)
  {
    G4ThreeVector alpha, beta, delta, normal;
    G4bool reflected, dead;
    G4double n, n2, transmission, vDotN;
    G4double paraComp, perpComp, paraMag, perpMag, x;
    G4ThreeVector paraPolarization, perpPolarization;

    if (flag == IN_INTERFACE_1 ||
        flag == IN_INTERFACE_2 ||
        flag == IN_INTERFACE_3 ||
        flag == IN_INTERFACE_4) {
      dead = PropagateInInterface<pmtId>(flag);  // flag needs to be altered in PiD
      if (dead)
        return true;
    }

    do {
      if (flag == INWARD) {
        n = fDomeRIndex.Y(fSampledMomentum);
        n2 = fInterfaceRIndex.Y(fSampledMomentum);
      } else {
        n = fInterfaceRIndex.Y(fSampledMomentum);
        n2 = fDomeRIndex.Y(fSampledMomentum);
      }
      normal.set(
        2 * fPhotonPosition.x() / fInterfaceRadiusSq,
        2 * fPhotonPosition.y() / fInterfaceRadiusSq,
        2 * fPhotonPosition.z() / fInterfaceRadiuszSq
      );
      normal *= 1/normal.mag();
      if (pmtId == 4)
        normal.set(0, 0, -1);
      vDotN = -normal.dot(fPhotonMomentum);  // is actually cos(theta)
      if (vDotN > 0.999999) {
        transmission = 4 * n * n2 / pow(n + n2, 2);
        reflected = (G4UniformRand() > transmission);
        if (reflected) {
          fPhotonMomentum = -fPhotonMomentum;
          fPhotonPolarization = -fPhotonPolarization;
          if (flag == INWARD) {
            flag = OUTWARD;
            return false;
          }
          flag = INWARD;
          dead = PropagateInInterface<pmtId>(flag);
          if (dead)
            return true;
        } else {
          if (flag == OUTWARD)
            return true;
          dead = PropagateInInterface<pmtId>(flag);
          if (dead)
            return true;
        }
      } else {
        alpha = normal.cross(fPhotonMomentum);
        alpha *= 1/alpha.mag();
        x = Sqr(n2) - Sqr(n) * (1 - Sqr(vDotN));
        if (x <= 0)
          x = 0;  // total reflection
        else
          x = (vDotN > 0) ? -sqrt(x) : sqrt(x);
        paraMag = fPhotonPolarization.dot(alpha);
        beta = fPhotonMomentum.cross(alpha); // should be unit vector
        perpMag = fPhotonPolarization.dot(beta);
        perpComp = perpMag / (n * vDotN - x);
        paraComp = paraMag / (n2 * vDotN - n * x / n2);
        transmission = -4 * n * vDotN * x * (Sqr(perpComp) + Sqr(paraComp));
        reflected = (G4UniformRand() >= transmission);
        if (x == 0)
          reflected = true;
        if (flag == INWARD) {
          if (reflected) {
            // reflected back into dome
            fPhotonMomentum += 2 * vDotN * normal;
            fPhotonMomentum *= 1/fPhotonMomentum.mag();
            beta = fPhotonMomentum.cross(alpha);
            fPhotonPolarization = paraComp * alpha + perpComp * beta;
            fPhotonPolarization *= 1/fPhotonPolarization.mag();
            flag = OUTWARD;
            return false;
          } else {
            // transmitted back into tank
            fPhotonMomentum *= n;
            fPhotonMomentum += (n * vDotN + x) * normal;
            fPhotonMomentum *= 1/fPhotonMomentum.mag();
            beta = fPhotonMomentum.cross(alpha);
            fPhotonPolarization = paraComp * alpha + perpComp * beta;
            fPhotonPolarization *= 1/fPhotonPolarization.mag();
            dead = PropagateInInterface<pmtId>(flag);
            if (dead)
              return true;
          }
        } else {
          // reflected off inside of dome - remains in interface
          if (reflected) {
            fPhotonMomentum += 2 * vDotN * normal;
            fPhotonMomentum *= 1/fPhotonMomentum.mag();
            beta = fPhotonMomentum.cross(alpha);
            fPhotonPolarization = paraComp * alpha + perpComp * beta;
            fPhotonPolarization *= 1/fPhotonPolarization.mag();
            flag = INWARD;
            dead = PropagateInInterface<pmtId>(flag);
            if (dead)
              return true;
          } else {
            // transmitted from interface back into dome
            fPhotonMomentum *= n;
            fPhotonMomentum += (n * vDotN + x) * normal;
            fPhotonMomentum *= 1/fPhotonMomentum.mag();
            beta = fPhotonMomentum.cross(alpha);
            fPhotonPolarization = paraComp * alpha + perpComp * beta;
            fPhotonPolarization *= 1/fPhotonPolarization.mag();
            return false;
          }
        }
      }
    } while (true);
  }


  template<int pmtId>
  G4bool
  G4StationFastCerenkov::PropagateInInterface(G4int& flag)
  {
    G4bool dead;

    while (pmtId == 4) {
      // We are inside the SPMT interface, that is above the dome level
      // So the z-component of the photon position relative to 
      // the upper part of the SPMT dome is > 0
      G4double distanceToDome = (fRoofPos[pmtId] + fDomeThick_SPMT - fPhotonPosition.z()) / fPhotonMomentum.z();
      G4double distanceToPMT = (fRoofPos[pmtId] + fDomeThick_SPMT + fInterfaceThick_SPMT - fPhotonPosition.z()) / fPhotonMomentum.z();
      
      G4ThreeVector displacement;
      G4double distance;

      if (flag == INWARD) {

        if (distanceToPMT < 0) {
          const std::string err = "Photon tracking inside SPMT interface run into an impossible condition : INWARD & distToPMT<0";
          ERROR(err);
          throw utl::OutOfBoundException(err);
        }

        displacement = fPhotonPosition + distanceToPMT * fPhotonMomentum;
        distance = (Sqr(displacement.x()) + Sqr(displacement.y())) / fPMTRadius_SPMT_Sq;
        if (distance < 1) {
          dead = (G4UniformRand() > exp(-distanceToPMT / fInterfaceAbsLength.Y(fSampledMomentum)));
          if (dead)
            return true;
          fPhotonPosition += distanceToPMT * fPhotonMomentum;
          fPhotonTime += distanceToPMT * fInterfaceRIndex.Y(fSampledMomentum) / c_light;
          
          distance = (Sqr(displacement.x()) + Sqr(displacement.y())) / fPMTActiveRadius_SPMT_Sq;
          if (distance < 1) {
            const G4double eff = fQE_SPMT.Y(fSampledMomentum);
            if (fPhotonMaxQE * G4UniformRand() <= eff)
              fG4StationSimulator.AddPhoton(pmtId, fPhotonTime);
            return true;  // photon absorbed, so is now dead
          }
          
          G4double maxDistanceInsidePMT = fPMTThick_SPMT / fPhotonMomentum.z();
          displacement = fPhotonPosition + maxDistanceInsidePMT * fPhotonMomentum;
          distance = (Sqr(displacement.x()) + Sqr(displacement.y())) / fPMTActiveRadius_SPMT_Sq;          
          if (distance < 1) {
            const G4double eff = fQE_SPMT.Y(fSampledMomentum);
            if (fPhotonMaxQE * G4UniformRand() <= eff)
              fG4StationSimulator.AddPhoton(pmtId, fPhotonTime);
          }
          return true;  // photon killed
          
        } else
          return true;  // photon killed

      } else {

        if (distanceToDome < 0) {
          const std::string err = "Photon tracking inside SPMT interface run into an impossible condition : OUTWARD & distToDome<0";
          ERROR(err);
          throw utl::OutOfBoundException(err);
        }

        displacement = fPhotonPosition + distanceToDome * fPhotonMomentum;
        distance = (Sqr(displacement.x()) + Sqr(displacement.y())) / fInterfaceRadius_SPMT_Sq;
        if (distance < 1) {
          dead = (G4UniformRand() > exp(-distanceToDome / fInterfaceAbsLength.Y(fSampledMomentum)));
          if (dead)
            return true;
          fPhotonPosition += distanceToDome * fPhotonMomentum;
          fPhotonTime += distanceToDome * fInterfaceRIndex.Y(fSampledMomentum) / c_light;
          return false;
        } else
          return true;  // photon killed
      }
    }

    G4double distanceToRoof = (fRoofPos[pmtId] - fPhotonPosition.z()) / fPhotonMomentum.z();
    if (distanceToRoof < 0)
      distanceToRoof = 1e99;

    // we start at interface surface, so only one non-zero solution exists
    //DistanceThis = -2.0 * fPhotonPosition.dot(fPhotonMomentum);
    const G4double distanceThis = GetEllipsoidIntersect(fPhotonPosition, fInterfaceRadius, fInterfaceRadiusz);
    const G4double distanceOther = GetEllipsoidIntersect(fPhotonPosition, fPMTRadius, fPMTRadiusz);

    if (distanceThis < distanceOther) {
      if (distanceToRoof < distanceThis)
        return true;  // photon killed
      dead = (G4UniformRand() > exp(-distanceThis / fInterfaceAbsLength.Y(fSampledMomentum)));
      if (dead)
        return true;
      fPhotonPosition += distanceThis * fPhotonMomentum;
      fPhotonTime += distanceThis * fInterfaceRIndex.Y(fSampledMomentum) / c_light;
      flag = OUTWARD;
      return false;
    } else {
      if (distanceToRoof < distanceOther)
        return true; // photon killed
      dead = (G4UniformRand() > exp(-distanceOther / fInterfaceAbsLength.Y(fSampledMomentum)));
      if (dead)
        return true;
      fPhotonTime += distanceOther * fInterfaceRIndex.Y(fSampledMomentum) / c_light;

      if ((-fPhotonPosition.z()) > fHeightz) {
        if (fHasSPMT) {
          const G4double eff = fQE.Y(fSampledMomentum);
          if (fPhotonMaxQE * G4UniformRand() <= eff)
            fG4StationSimulator.AddPhoton(pmtId, fPhotonTime);  
        } else
          fG4StationSimulator.AddPhoton(pmtId, fPhotonTime);
      }
      return true;  // photon absorbed, so is now dead
    }
  }


  // Scattering routines

  G4bool
  G4StationFastCerenkov::ScatterOffRoof()
  {
    const G4ThreeVector normal(0, 0, -1);
    const G4double reflectivity = fLinerReflectivity.Y(fSampledMomentum);
    G4double r = G4UniformRand();
    if (r > reflectivity)
      return true;  // absorbed
    r /= reflectivity;  // range 0..1
    const G4double lobe = fLinerSpecularLobe.Y(fSampledMomentum);
    if (r <= lobe) {
      LobeScatterHorizontal(normal);
      return false;
    }
    const G4double spike = lobe + fLinerSpecularSpike.Y(fSampledMomentum);
    if (r <= spike) {
      SpikeScatter(normal);
      return false;
    }
    const G4double back = spike + fLinerBackscatter.Y(fSampledMomentum);
    if (r <= back) {
      BackScatter();
      return false;
    }
    DiffuseScatterHorizontal(normal);
    return false;
  }


  G4bool
  G4StationFastCerenkov::ScatterOffWall()
  {
    const G4double reflectivity = fLinerReflectivity.Y(fSampledMomentum);
    G4double r = G4UniformRand();
    if (r > reflectivity)
      return true;  // absorbed
    r /= reflectivity;  // range 0..1
    const G4double x = fPhotonPosition.x();
    const G4double y = fPhotonPosition.y();
    const G4double w = 1 / sqrt(Sqr(x) + Sqr(y));
    const G4ThreeVector normal(-x * w, -y * w, 0);
    const G4double lobe = fLinerSpecularLobe.Y(fSampledMomentum);
    if (r <= lobe) {
      LobeScatterVertical(normal);
      return false;
    }
    const G4double spike = lobe + fLinerSpecularSpike.Y(fSampledMomentum);
    if (r <= spike) {
      SpikeScatter(normal);
      return false;
    }
    const G4double back = spike + fLinerBackscatter.Y(fSampledMomentum);
    if (r <= back) {
      BackScatter();
      return false;
    }
    DiffuseScatterVertical(normal);
    return false;
  }


  G4bool
  G4StationFastCerenkov::ScatterOffFloor()
  {
    const G4ThreeVector normal(0, 0, 1);
    const G4double reflectivity = fLinerReflectivity.Y(fSampledMomentum);
    G4double r = G4UniformRand();
    if (r > reflectivity)
      return true;  // absorbed
    r /= reflectivity;  // range 0..1
    const G4double lobe = fLinerSpecularLobe.Y(fSampledMomentum);
    if (r <= lobe) {
      LobeScatterHorizontal(normal);
      return false;
    }
    const G4double spike = lobe + fLinerSpecularSpike.Y(fSampledMomentum);
    if (r <= spike) {
      SpikeScatter(normal);
      return false;
    }
    const G4double back = spike + fLinerBackscatter.Y(fSampledMomentum);
    if (r <= back) {
      BackScatter();
      return false;
    }
    DiffuseScatterHorizontal(normal);
    return false;
  }


  void
  G4StationFastCerenkov::BackScatter()
  {
    //fPhotonMomentum.set(-fPhotonMomentum.x(), -fPhotonMomentum.y(), -fPhotonMomentum.z());
    fPhotonMomentum *= -1;
    // can leave polarisation as it is
  }


  void
  G4StationFastCerenkov::DiffuseScatterVertical(const G4ThreeVector& normal)
  {
    // for scattering off the tank walls
    const G4double rand = G4UniformRand();
    const G4double cosTheta = sqrt(1 - rand);
    const G4double sinTheta = sqrt(rand);
    const G4double phi = 2 * utl::kPi * G4UniformRand();
    const G4double sinThetaCosPhi = sinTheta * cos(phi);
    const G4double sinThetaSinPhi = sinTheta * sin(phi);
    const G4double x = normal.x();
    const G4double y = normal.y();
    const G4double w = sqrt(Sqr(x) + Sqr(y));
    const G4double vx = x * cosTheta - y * sinThetaSinPhi / w;
    const G4double vy = y * cosTheta + x * sinThetaSinPhi / w;
    const G4double vz = -w * sinThetaCosPhi;
    fPhotonMomentum.set(vx, vy, vz);
    const G4double dotProduct = fPhotonMomentum.dot(fPhotonPolarization);
    fPhotonPolarization = fPhotonPolarization - dotProduct * fPhotonMomentum;
    fPhotonPolarization /= fPhotonPolarization.mag();
  }


  void
  G4StationFastCerenkov::DiffuseScatterHorizontal(const G4ThreeVector& normal)
  {
    // for scattering off tank floor or roof
    const G4double rand = G4UniformRand();
    const G4double cosTheta = sqrt(1 - rand);
    const G4double sinTheta = sqrt(rand);
    const G4double phi = 2 * utl::kPi * G4UniformRand();
    const G4double vx = sinTheta * cos(phi);
    const G4double vy = sinTheta * sin(phi);
    const G4double vz = normal.z() * cosTheta;  // Normal.z() = +/- 1
    fPhotonMomentum.set(vx, vy, vz);
    const G4double dotProduct = fPhotonMomentum.dot(fPhotonPolarization);
    fPhotonPolarization = fPhotonPolarization - dotProduct * fPhotonMomentum;
    fPhotonPolarization /= fPhotonPolarization.mag();
  }


  void
  G4StationFastCerenkov::SpikeScatter(const G4ThreeVector& normal)
  {
    G4double dotProduct = fPhotonMomentum.dot(normal);  // less than zero;
    fPhotonMomentum = fPhotonMomentum - (2 * dotProduct) * normal;
    dotProduct = fPhotonMomentum.dot(fPhotonPolarization);
    fPhotonPolarization = fPhotonPolarization - dotProduct * fPhotonMomentum;
    fPhotonMomentum /= fPhotonMomentum.mag();
    fPhotonPolarization /= fPhotonPolarization.mag();
  }


  void
  G4StationFastCerenkov::LobeScatterVertical(const G4ThreeVector& normal)
  {
    const G4double fMax = (fSigmaAlpha < 0.25) ? 4 * fSigmaAlpha : 1;
    const G4double x = normal.x();
    const G4double y = normal.y();
    const G4double w = sqrt(Sqr(x) + Sqr(y));
    G4ThreeVector trialMomentum;
    do {
      G4ThreeVector facetNormal;
      G4double dotProduct;
      do {
        G4double alpha;
        G4double sinAlpha;
        do {
          alpha = G4RandGauss::shoot(0, fSigmaAlpha);
          sinAlpha = sin(alpha);
        } while (G4UniformRand()*fMax > sinAlpha || alpha >= 0.5*utl::kPi);
        const G4double phi = 2*utl::kPi * G4UniformRand();
        const G4double cosAlpha = cos(alpha);
        const G4double sinAlphaSinPhi = sinAlpha * sin(phi);
        const G4double sinAlphaCosPhi = sinAlpha * cos(phi);
        const G4double vx = x * cosAlpha - y * sinAlphaSinPhi / w;
        const G4double vy = y * cosAlpha + x * sinAlphaSinPhi / w;
        const G4double vz = -w * sinAlphaCosPhi;
        facetNormal.set(vx, vy, vz);
        dotProduct = fPhotonMomentum.dot(facetNormal);
      } while (dotProduct >= 0);
      trialMomentum = fPhotonMomentum - (2 * dotProduct) * facetNormal;
    } while (trialMomentum.dot(normal) <= 0);

    fPhotonMomentum = trialMomentum;
    const G4double dotProduct = fPhotonMomentum.dot(fPhotonPolarization);
    fPhotonPolarization = fPhotonPolarization - dotProduct * fPhotonMomentum;
    fPhotonMomentum /= fPhotonMomentum.mag();
    fPhotonPolarization /= fPhotonPolarization.mag();
  }


  void
  G4StationFastCerenkov::LobeScatterHorizontal(const G4ThreeVector& normal)
  {
    const G4double fMax = (fSigmaAlpha < 0.25) ? 4 * fSigmaAlpha : 1;
    G4ThreeVector trialMomentum;
    do {
      G4ThreeVector facetNormal;
      G4double dotProduct;
      do {
        G4double alpha;
        G4double sinAlpha;
        do {
          alpha = G4RandGauss::shoot(0, fSigmaAlpha);
          sinAlpha = sin(alpha);
        } while (G4UniformRand()*fMax > sinAlpha || alpha >= 0.5*utl::kPi);
        const G4double phi = 2*utl::kPi * G4UniformRand();
        const G4double cosAlpha = cos(alpha);
        const G4double vx = sinAlpha * cos(phi);
        const G4double vy = sinAlpha * sin(phi);
        const G4double vz = normal.z() * cosAlpha;  // Normal.z() = +/- 1
        facetNormal.set(vx, vy, vz);
        dotProduct = fPhotonMomentum.dot(facetNormal);
      } while (dotProduct >= 0);
      trialMomentum = fPhotonMomentum - (2 * dotProduct) * facetNormal;
    } while (trialMomentum.dot(normal) <= 0);

    fPhotonMomentum = trialMomentum;
    const G4double dotProduct = fPhotonMomentum.dot(fPhotonPolarization);
    fPhotonPolarization = fPhotonPolarization - dotProduct * fPhotonMomentum;
    fPhotonMomentum /= fPhotonMomentum.mag();
    fPhotonPolarization /= fPhotonPolarization.mag();
  }


  // Utility routines

  G4double
  G4StationFastCerenkov::GetSphereIntersect(const G4ThreeVector& pos,
                                         const G4double r)
    const
  {
    const G4double b = -pos.dot(fPhotonMomentum);
    const G4double c = pos.dot(pos) - Sqr(r);
    G4double d = Sqr(b) - c;
    if (d < 0)
      return 1e99;  // no intersection with sphere at all

    d = sqrt(d);
    const G4double distA = b + d;
    const G4double distB = b - d;
    if (distA > 0) {
      if (distB < distA) {
        if (distB > 0)
          return distB;
        return distA;  // inside sphere - B is behind us, A ahead
      }
      return distA;
    } else {
      if (distB > 0)
        return distB;  // inside sphere - A is behind us, B ahead
    }
    // if we get here, both are < 0 - intersections are both behind us already
    return 1e99;
  }


  G4double
  G4StationFastCerenkov::GetEllipsoidIntersect(const G4ThreeVector& pos,
                                            const G4double r1, const G4double r2)
    const
  {
    // The equation of the line with the photon direction is given by Pos+a*fPhotonMomentum being a the
    // distance from the photon position since fPhotonMomentum is a unitary vector in the direction of
    // the photon

    const double b[3] = { Sqr(fPhotonMomentum.x()), Sqr(fPhotonMomentum.y()), Sqr(fPhotonMomentum.z()) };
    const double c[3] = { Sqr(pos.x()), Sqr(pos.y()), Sqr(pos.z()) };
    const double d[3] = { fPhotonMomentum.x() * pos.x(), fPhotonMomentum.y() * pos.y(), fPhotonMomentum.z() * pos.z() };

    const double r12 = Sqr(r1);
    const double r22 = Sqr(r2);
    const double e = (b[0] + b[1]) / r12 + b[2] / r22;
    const double f = 2*((d[0] + d[1]) / r12 + d[2] / r22);
    const double g = (c[0] + c[1]) / r12 + c[2] / r22 - 1;

    const double dis = Sqr(f) - 4*e*g;
    // if no solution of the quadratic equation there is no intersection
    if (dis < 0)
      return 1e99;

    const double sq = sqrt(dis);
    // take smallest positive value of the solution
    const double i2e = 1 / (2*e);
    const double a1 = i2e * (-f + sq);
    const double a2 = i2e * (-f - sq);
    if (a2 > 1e-4) // to avoid almost 0 but positive a2 for current photon position
      return a2;

    if (a1 > 0)
      return a1;

    return 1e99;
  }

}