G4XTankFastCerenkov.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: Propogate.icc

#include "G4XTankFastCerenkov.h"
#include "G4XTankConstruction.h"
#include "G4XTankPMT.h"
using namespace G4XTankSimulatorAG;

#include <G4Version.hh>

#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 <iostream>
#include <vector>
#include <string>

double G4XTankFastCerenkov::fTankHeight;
double G4XTankFastCerenkov::fTankHalfHeight;
double G4XTankFastCerenkov::fTankThickness;
double G4XTankFastCerenkov::fTankRadius;
double G4XTankFastCerenkov::fTankRadiusSq;
double G4XTankFastCerenkov::fDomeRadius;
double G4XTankFastCerenkov::fDomeRadiusz;
double G4XTankFastCerenkov::fInterfaceRadius;
double G4XTankFastCerenkov::fInterfaceRadiusz;
double G4XTankFastCerenkov::fPMTRadius;
double G4XTankFastCerenkov::fPMTRadiusz;
double G4XTankFastCerenkov::fDomeRadiusSq;
double G4XTankFastCerenkov::fDomeRadiuszSq;
double G4XTankFastCerenkov::fInterfaceRadiusSq;
double G4XTankFastCerenkov::fInterfaceRadiuszSq;
double G4XTankFastCerenkov::fPMTRadiusSq;
double G4XTankFastCerenkov::fPMTRadiuszSq;
double G4XTankFastCerenkov::fHeightz;


double G4XTankFastCerenkov::fSigmaAlpha;

utl::TabulatedFunction G4XTankFastCerenkov::fQE;
utl::TabulatedFunction G4XTankFastCerenkov::fWaterAbsLength;

utl::TabulatedFunction G4XTankFastCerenkov::fInterfaceRIndex;
utl::TabulatedFunction G4XTankFastCerenkov::fInterfaceAbsLength;

utl::TabulatedFunction G4XTankFastCerenkov::fDomeRIndex;
utl::TabulatedFunction G4XTankFastCerenkov::fDomeAbsLength;

utl::TabulatedFunction G4XTankFastCerenkov::fLinerReflectivity;
utl::TabulatedFunction G4XTankFastCerenkov::fLinerSpecularLobe;
utl::TabulatedFunction G4XTankFastCerenkov::fLinerSpecularSpike;
utl::TabulatedFunction G4XTankFastCerenkov::fLinerBackscatter;
G4ThreeVector G4XTankFastCerenkov::fPMTPos[4];


double G4XTankFastCerenkov::fRoofPos[4];


G4XTankFastCerenkov::G4XTankFastCerenkov(const G4String& processName)
  : G4VContinuousProcess(processName)
{
  fTrackSecondariesFirst = false;
  fMaxPhotons = 0;
  fThePhysicsTable = NULL;
  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 CLHEP header PhysicalConstants.h - included from globals.hh
  fRayleighConst = 27.0 * hc4 / (8.0 * pow(utl::kPi, 3) * betat *
                                 temp * kboltz);

}

G4XTankFastCerenkov::~G4XTankFastCerenkov() {

  if (fThePhysicsTable != NULL) {

    fThePhysicsTable->clearAndDestroy();
    delete fThePhysicsTable;
    fThePhysicsTable = 0;
  
  }

}

void G4XTankFastCerenkov::GetDataFromConstruction() {

  fTankRadius = G4XTankConstruction::fTankRadius;
  fTankHalfHeight = G4XTankConstruction::fTankHalfHeight;
  fTankHeight = fTankHalfHeight * 2.0;
  fTankThickness = G4XTankConstruction::fTankThickness;
 
  fPMTRadius = G4XTankConstruction::fPmtRmax;
  fPMTRadiusz = G4XTankConstruction::fPmtRzmax;
  
  fInterfaceRadius = G4XTankConstruction::fInterfaceRmax;
  fInterfaceRadiusz = G4XTankConstruction::fInterfaceRzmax;
  fInterfaceRIndex = G4XTankConstruction::fInterfaceRINDEX;
  fInterfaceAbsLength = G4XTankConstruction::fInterfaceABSORPTION;

  fDomeRadius = G4XTankConstruction::fDomeRmax;
  fDomeRadiusz = G4XTankConstruction::fDomeRzmax;
  fDomeRIndex = G4XTankConstruction::fPmtdomeRINDEX;
  fDomeAbsLength = G4XTankConstruction::fPmtdomeABSORPTION;
  
  fHeightz = G4XTankConstruction::fHeightz;
  fWaterAbsLength = G4XTankConstruction::fWaterABSORPTION;

  fLinerReflectivity = G4XTankConstruction::fLinerREFLECTIVITY;
  fLinerSpecularLobe = G4XTankConstruction::fLinerSPECULARLOBECONSTANT;
  fLinerSpecularSpike = G4XTankConstruction::fLinerSPECULARSPIKECONSTANT;
  fLinerBackscatter = G4XTankConstruction::fLinerBACKSCATTERCONSTANT;
  fSigmaAlpha = G4XTankConstruction::fSIGMA_ALPHA;
  
  const sdet::Station *fCurrentDetectorStation = G4XTankSimulator::GetCurrentDetectorStation();
  
  fQE = fCurrentDetectorStation->GetPMT(1).GetQuantumEfficiency();
  
  const double& collectionEfficiency =
    fCurrentDetectorStation->GetPMT(1).GetCollectionEfficiency();
  


  std::vector<double>::iterator fIt;
  for (fIt = fQE.YBegin(); fIt != fQE.YEnd(); fIt++) 
    (*fIt) *= collectionEfficiency; // checked - this is right
  
  for (fIt = fQE.XBegin(); fIt != fQE.XEnd(); fIt++) 
    (*fIt) *= eV;
  
  // We need to check that values have been imported correctly, and that
  // e.g. water absorption length is correct and doesn't need scaling
  // by a MaxAbsLength parameter (likewise quantum efficiency).
  // NOW CHECKED - AbsLength is the actual absorption length, fQE is the
  // actual fQE (not including collection efficiency)
  //
  // Also need to see what the Pmt positions are in the utl::Point structures.
  // (need to specify coordinate system when accesing co-ordinate values)
  
  const utl::Point& pmt1 = G4XTankConstruction::fPmt1;
  const utl::Point& pmt2 = G4XTankConstruction::fPmt2;
  const utl::Point& pmt3 = G4XTankConstruction::fPmt3;

  const utl::CoordinateSystemPtr& pmt1Coord= pmt1.GetCoordinateSystem(); 
  const utl::CoordinateSystemPtr& pmt2Coord= pmt2.GetCoordinateSystem(); 
  const utl::CoordinateSystemPtr& pmt3Coord= pmt3.GetCoordinateSystem(); 

  fPMTPos[1].set(pmt1.GetX(pmt1Coord) * m,
                 pmt1.GetY(pmt1Coord) * m,
                 pmt1.GetZ(pmt1Coord) * m);

  fPMTPos[2].set(pmt2.GetX(pmt2Coord) * m,
                 pmt2.GetY(pmt2Coord) * m,
                 pmt2.GetZ(pmt2Coord) * m);

  fPMTPos[3].set(pmt3.GetX(pmt3Coord) * m,
                 pmt3.GetY(pmt3Coord) * m,
                 pmt3.GetZ(pmt3Coord) * m);

  for(unsigned int k=1; k<4; ++k)
    fRoofPos[k] = fTankHeight + fTankThickness - fPMTPos[k].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;
  return;
  
}

G4bool G4XTankFastCerenkov::IsApplicable(const G4ParticleDefinition& aParticle) {
  
  return (aParticle.GetPDGCharge() != 0);

}

void G4XTankFastCerenkov::SetTrackSecondariesFirst(const G4bool state) {
  
  fTrackSecondariesFirst = state;
  
}

void G4XTankFastCerenkov::SetMaxNumPhotonsPerStep(const G4int numPhotons) {

  fMaxPhotons = numPhotons;

}

G4PhysicsTable *G4XTankFastCerenkov::GetPhysicsTable() const {

  return fThePhysicsTable;

}

void G4XTankFastCerenkov::DumpPhysicsTable() const {

  G4int PhysicsTableSize = fThePhysicsTable->entries();
  G4PhysicsOrderedFreeVector *v;
  //G4int i;
  
  for (G4int i = 0; i < PhysicsTableSize; ++i) {

    v = static_cast<G4PhysicsOrderedFreeVector *>((*fThePhysicsTable)[i]);
    v->DumpValues();

  }

}

G4VParticleChange* G4XTankFastCerenkov::AlongStepDoIt(const G4Track& aTrack,
                                                     const G4Step& aStep) {
    
  //  std::vector<double>::iterator VecItPP;

  aParticleChange.Initialize(aTrack);  
  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
  const G4Material* aMaterial = aTrack.GetMaterial();
  G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();

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

  if (t0 > 1.0e6) return G4VContinuousProcess::AlongStepDoIt(aTrack, aStep);

  G4ThreeVector aPrimaryDirection = pPreStepPoint->GetMomentumDirection();
  
  const G4double cosAlpha = aPrimaryDirection.z();
  const G4double sinAlpha = sqrt(1.0 - (cosAlpha * cosAlpha));
  const G4double cosBeta = aPrimaryDirection.x() / sinAlpha;
  const G4double sinBeta = aPrimaryDirection.y() / sinAlpha;

  G4MaterialPropertiesTable* aMaterialPropertiesTable =
    aMaterial->GetMaterialPropertiesTable();
  
  if (!aMaterialPropertiesTable)
    return G4VContinuousProcess::AlongStepDoIt(aTrack, aStep);

  const G4MaterialPropertyVector* Rindex = 
    aMaterialPropertiesTable->GetProperty("RINDEX"); 
  
  if (!Rindex) 
    return G4VContinuousProcess::AlongStepDoIt(aTrack, aStep);

  G4double MeanNumPhotons = GetAverageNumberOfPhotons(aParticle, aMaterial,
                                                      Rindex);
  if (MeanNumPhotons <= 0.0) {
    aParticleChange.SetNumberOfSecondaries(0);
    return G4VContinuousProcess::AlongStepDoIt(aTrack, aStep);
  }

  MeanNumPhotons *= aStep.GetStepLength();
  
  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(aTrack, aStep);
  
#if (G4VERSION_NUMBER < 920) 
    const G4double Pmin = Rindex->GetMinPhotonMomentum();
    const G4double Pmax = Rindex->GetMaxPhotonMomentum();
#else
    const G4double Pmin = Rindex->GetMinPhotonEnergy();
    const G4double Pmax = Rindex->GetMaxPhotonEnergy();
#endif

  const G4double dp = Pmax - Pmin;  
  const G4double nMax = Rindex->GetMaxProperty();
  
  const G4double betaInverse = aParticle->GetTotalEnergy() /
    aParticle->GetTotalMomentum();
  
  const G4double maxCos = betaInverse / nMax; 
  const G4double maxSin2 = (1.0 - (maxCos * maxCos));

  const double minQEMomentum = fQE[0].X();
  const double maxQEMomentum = fQE[fQE.GetNPoints() - 1].X();

  // leave in stuff for dealing with change of refractive index with
  // wavelength - even if the current setup ignores dispersion in water...

  // variables used in the loop later on
  G4double rand, cosTheta, sin2Theta;
  G4double phi, sinTheta, px, py, pz;
  G4double invSumMFP;
  G4double delta, deltaTime;
  G4double sinThetaSinPhi, sinThetaCosPhi;
  G4double radiusTest;
  double absorptionMFP; 
  double rayleighMFP; 
  for (G4int i = 0; i < NumPhotons; ++i) {
    do {
      rand = G4UniformRand();
      fSampledMomentum = Pmin + rand * dp; 
      fSampledRI = Rindex->GetProperty(fSampledMomentum);
      cosTheta = betaInverse / fSampledRI;  
      sin2Theta = (1.0 - (cosTheta * cosTheta));
      rand = G4UniformRand();
    } while ((rand * maxSin2) > sin2Theta);

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

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

    if (G4UniformRand() > fQE.Y(fSampledMomentum))
      continue; // Photon is ignored
    
    // If we get here, we want to track the photon.
    G4double SampledRI_Sq, RIConst;

    // calculate absorption and Rayleigh scattering constants for this photon
    int iFirstValidBin_Wat = 0;
    int iLastValidBin_Wat = fWaterAbsLength.GetNPoints() - 1;
    
    // If wrong... ignore it
    if(fSampledMomentum> fWaterAbsLength[iFirstValidBin_Wat].X() &&
       fSampledMomentum< fWaterAbsLength[iLastValidBin_Wat].X())
      absorptionMFP = fWaterAbsLength.Y(fSampledMomentum); // CHECK
    else continue;

    
    SampledRI_Sq = fSampledRI * fSampledRI;
    RIConst = (SampledRI_Sq + 2.0) * (SampledRI_Sq - 1.0);

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

    invSumMFP = 1.0 / (absorptionMFP + rayleighMFP);

    fInvMFP = (1.0 / absorptionMFP) + (1.0 / rayleighMFP);
    fRayleighFrac = absorptionMFP * invSumMFP;
    fAbsorbFrac = rayleighMFP * invSumMFP;

    // calculate photon vectors (direction, polarisation)
    sinTheta = sqrt(sin2Theta); 
    rand = G4UniformRand();
   
    phi = 2.0 * utl::kPi * rand;
    sinThetaSinPhi = sin(phi) * sinTheta;
    sinThetaCosPhi = cos(phi) * sinTheta;
    
    px = (cosAlpha * cosBeta * sinThetaCosPhi) - (sinBeta * sinThetaSinPhi)
      + (cosBeta * sinAlpha * cosTheta);
    py = (cosAlpha * sinBeta * sinThetaCosPhi) + (cosBeta * sinThetaSinPhi)
      + (sinBeta * sinAlpha * cosTheta);
    pz = (cosAlpha * cosTheta) - (sinAlpha * sinThetaCosPhi);
    
    fPhotonMomentum.set(px, py, pz); // This is a momentum direction
    fPhotonPolarization = (aPrimaryDirection - (cosTheta * fPhotonMomentum));
    fPhotonPolarization /= fPhotonPolarization.mag();
  
    rand = G4UniformRand();
    delta = rand * aStep.GetStepLength();
    deltaTime = delta / ((pPreStepPoint->GetVelocity() +
                          pPostStepPoint->GetVelocity()) / 2.0);
    fPhotonTime = t0 + deltaTime;
    fPhotonPosition = x0 + rand * aStep.GetDeltaPosition();
    
    radiusTest = fPhotonPosition.x() * fPhotonPosition.x() +
      fPhotonPosition.y() * fPhotonPosition.y() - fTankRadiusSq;
    
    // if Cerenkov photon is inside tank, then track it
    if ((radiusTest < 0.0) && (fPhotonPosition.z() > fTankThickness) &&
        (fPhotonPosition.z() < fTankHeight + fTankThickness)) 
      TrackPhotonInTank();
  }// for photons
  
  if (verboseLevel > 0) {
    G4cerr << "\n Exiting from G4Cerenkov::DoIt -- NumberOfSecondaries = " 
           << aParticleChange.GetNumberOfSecondaries() << G4endl;
  }

  return G4VContinuousProcess::AlongStepDoIt(aTrack, aStep);


}

void G4XTankFastCerenkov::BuildThePhysicsTable() {

  if (fThePhysicsTable) 
    return;
  
  
  
  const G4MaterialTable* theMaterialTable= G4Material::GetMaterialTable();
  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
  fThePhysicsTable = new G4PhysicsTable(numOfMaterials);
  
  G4PhysicsOrderedFreeVector *aPhysicsOrderedFreeVector;
  G4Material *aMaterial;
  G4MaterialPropertiesTable *aMaterialPropertiesTable;
  G4MaterialPropertyVector *theRIndexVector;
  G4double currentRI, currentPM, currentCAI;
  G4double prevRI, prevPM, prevCAI;
  
  for (G4int i = 0; i < numOfMaterials; ++i) {
    
    aPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();
    aMaterial = (*theMaterialTable)[i];
    aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
    
    if (aMaterialPropertiesTable) {
      
      theRIndexVector =  aMaterialPropertiesTable->GetProperty("RINDEX");
      
      if (theRIndexVector) {
        
        theRIndexVector->ResetIterator();
        ++(*theRIndexVector);
        currentRI = theRIndexVector->GetProperty();
        
        if (currentRI > 1.0) {
          
#if (G4VERSION_NUMBER < 920) 
            currentPM = theRIndexVector->GetPhotonMomentum();
#else 
            currentPM = theRIndexVector->GetPhotonEnergy();
#endif

          currentCAI = 0.0;
          aPhysicsOrderedFreeVector->InsertValues(currentPM , currentCAI);
          prevRI  = currentRI;
          prevPM  = currentPM;
          prevCAI = currentCAI;
    
          while(++(*theRIndexVector)) {
                    
            currentRI = theRIndexVector->GetProperty();

#if (G4VERSION_NUMBER < 920)
            currentPM = theRIndexVector->GetPhotonMomentum();
#else
            currentPM = theRIndexVector->GetPhotonEnergy();
#endif

            currentCAI = 0.5 * ((1.0 / (prevRI * prevRI)) +
                                (1.0 / (currentRI * currentRI)));
            currentCAI = prevCAI + ((currentPM - prevPM) * currentCAI);
            aPhysicsOrderedFreeVector->InsertValues(currentPM, currentCAI);
            prevPM  = currentPM;
            prevCAI = currentCAI;
            prevRI  = currentRI;
                    
          }
                
        }
            
      }
        
    }

    fThePhysicsTable->insertAt(i, aPhysicsOrderedFreeVector); 

  }

}

G4double 
G4XTankFastCerenkov::GetContinuousStepLimit(const G4Track& aTrack,
                                           G4double, G4double, G4double &) {

  if (fMaxPhotons <= 0) 
    return DBL_MAX;
  
  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
  const G4Material* aMaterial = aTrack.GetMaterial();
  
  G4MaterialPropertiesTable* aMaterialPropertiesTable =
    aMaterial->GetMaterialPropertiesTable();
  
  if (!aMaterialPropertiesTable) 
    return DBL_MAX;
  
  const G4MaterialPropertyVector* Rindex =
    aMaterialPropertiesTable->GetProperty("RINDEX");
  
  if (!Rindex) 
    return DBL_MAX;
  
  G4double meanNumPhotons = GetAverageNumberOfPhotons(aParticle, aMaterial,
                                                      Rindex);
  if (meanNumPhotons <= 0.0) 
    return DBL_MAX;
  
  G4double StepLimit = fMaxPhotons / meanNumPhotons;
  return StepLimit;

}

G4double 
G4XTankFastCerenkov::GetAverageNumberOfPhotons(const G4DynamicParticle* aParticle,
                                              const G4Material* aMaterial,
                                              const G4MaterialPropertyVector* Rindex) const {  
  const G4double Rfact = 369.81 / (eV * cm);
  
  if (aParticle->GetTotalMomentum() <= 0.0) 
    return 0.0;
  
  G4double betaInverse = aParticle->GetTotalEnergy() /
    aParticle->GetTotalMomentum();
  
  G4int materialIndex = aMaterial->GetIndex();
  
  G4PhysicsOrderedFreeVector* CerenkovAngleIntegrals =
    (G4PhysicsOrderedFreeVector*)((*fThePhysicsTable)(materialIndex));
  
  if (!(CerenkovAngleIntegrals->IsFilledVectorExist())) 
    return 0.0;

#if (G4VERSION_NUMBER < 920) 
    G4double Pmin = Rindex->GetMinPhotonMomentum();
    G4double Pmax = Rindex->GetMaxPhotonMomentum();
#else
    G4double Pmin = Rindex->GetMinPhotonEnergy();
    G4double Pmax = Rindex->GetMaxPhotonEnergy();
#endif

  G4double nMin = Rindex->GetMinProperty();     
  G4double nMax = Rindex->GetMaxProperty();
  G4double CAImax = CerenkovAngleIntegrals->GetMaxValue();
  G4double dp, ge;
  
  if (nMax < betaInverse) {
    
    dp = 0.0;
    ge = 0.0;
  
  } else if (nMin > betaInverse) {
  
    dp = Pmax - Pmin;   
    ge = CAImax; 
  
  } else {

#if (G4VERSION_NUMBER < 920)  
    Pmin = Rindex->GetPhotonMomentum(betaInverse);
#else
    Pmin = Rindex->GetPhotonEnergy(betaInverse);
#endif
    dp = Pmax - Pmin;
    G4bool isOutRange;
    G4double CAImin = CerenkovAngleIntegrals->GetValue(Pmin, isOutRange);
    ge = CAImax - CAImin;
    
    if (verboseLevel>0) {
      G4cerr << "CAImin = " << CAImin << G4endl;
      G4cerr << "ge = " << ge << G4endl;
    }

  }

  G4double charge = aParticle->GetDefinition()->GetPDGCharge();
  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 G4XTankFastCerenkov::DoRayleighScatter() {

  const G4double rand = 2.0 * G4UniformRand() - 1.0;
  const G4double cosTheta= (rand < 0.0)? 
    (-pow(-rand, 1.0 / 3.0)): pow(rand, 1.0 / 3.0);
  const G4double sinTheta = sqrt(1.0 - cosTheta * cosTheta);
  const G4double phi = 2.0 * 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.0 / sqrt(1.0 - 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.0 / k;

  //    cerr << "Matrix\n";
  //    cerr << R11 << "\t" << R12 << "\t" << a << endl;
  //    cerr << R21 << "\t" << R22 << "\t" << b << endl;
  //    cerr << "0.000" << "\t" << R32 << "\t" << c << endl;
  G4double x, y, z;
  
  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));  
    //      cerr << "Rayleigh1 " << x << " " << y << " " << z << endl;
  } 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;
    //      cerr << "Rayleigh2 " << x << " " << y << " " << z << endl;
    fPhotonMomentum.set((R11 * x + R12 * y + a * z),
                        (R21 * x + R22 * y + b * z),
                        (R32 * y + c * z));  
  }

}

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

void G4XTankFastCerenkov::TrackPhotonInTank(void)
{

  // 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 PropogateInTank 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.0) &&
      (fPhotonPosition.z() < (fTankHeight + fTankThickness - fDomeRadiusz))) {
    PropogateInTank(DOWN_IN_TANK);
    return;
  }
    
  G4double closest; 
  G4ThreeVector displacement;

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

  return;
}

void G4XTankFastCerenkov::PropogateInTank(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)) {
      //      G4cerr << "IN_DOME_2 or IN_INTERFACE_2" << G4endl;
      Dead = TransitionToDome<2>(flag);
      if (Dead) return;
    }
    else { // IN_DOME_3 or IN_INTERFACE_3
      //      G4cerr << "IN_DOME_3 or IN_INTERFACE_3" << G4endl;
      Dead = TransitionToDome<3>(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;
  G4double alpha, beta, gamma;
  G4double TestVal;
  G4bool repeat;
  G4double closest1, closest2, closest3;
  G4ThreeVector closest, Displacement1, Displacement2, Displacement3;
  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.0 / (fPhotonMomentum.x() * fPhotonMomentum.x() +
                     fPhotonMomentum.y() * fPhotonMomentum.y());
      alpha = (fPhotonMomentum.x() * fPhotonPosition.x() +
               fPhotonMomentum.y() * fPhotonPosition.y()) * gamma;
      beta = (fPhotonPosition.x() * fPhotonPosition.x() +
              fPhotonPosition.y() * fPhotonPosition.y() -
              fTankRadiusSq) * gamma;
     
      const double distance=alpha * 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=(DistanceToWall < DistanceToFloor)? 
        (1.0 - exp(-DistanceToWall * fInvMFP)):(1.0 - exp(-DistanceToFloor * fInvMFP));

      if (rand < TestVal) {
        // some interaction has occurred before reaching destination
        if (rand < (TestVal * fRayleighFrac)) {
          // we have a Rayleigh scatter
          DistanceTravelled = -log(1.0 - (rand / fRayleighFrac)) / fInvMFP;
          fPhotonPosition += DistanceTravelled * fPhotonMomentum;
          fPhotonTime += DistanceTravelled * fSampledRI / c_light;
          DoRayleighScatter(); // to pick new polarisation and direction
          if ((fPhotonMomentum.z() < 0.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;
      if (fPhotonMomentum.z() > 0)
        DistanceToRoof = (fTankHeight + fTankThickness - fPhotonPosition.z()) /
          fPhotonMomentum.z();
      else DistanceToRoof = -(fPhotonPosition.z()-fTankThickness) / fPhotonMomentum.z();// floor
        
      gamma = 1.0 / (fPhotonMomentum.x() * fPhotonMomentum.x() +
                     fPhotonMomentum.y() * fPhotonMomentum.y());
      alpha = (fPhotonMomentum.x() * fPhotonPosition.x() +
               fPhotonMomentum.y() * fPhotonPosition.y()) * gamma;
      beta = (fPhotonPosition.x() * fPhotonPosition.x() +
              fPhotonPosition.y() * fPhotonPosition.y() -
              fTankRadiusSq) * gamma;
      DistanceToWall = sqrt(alpha * alpha - beta) - alpha;
        
      Displacement1 = fPhotonPosition - fPMTPos[1];
      DotProduct = Displacement1.dot(fPhotonMomentum);
      closest = Displacement1 - (DotProduct * fPhotonMomentum);
      closest1 = (closest.x()*closest.x()+closest.y()*closest.y())/fDomeRadiusSq +
        (closest.z()*closest.z())/fDomeRadiuszSq;
      if (closest1 < 1. ) HitPMT += 1;
        
      Displacement2 = fPhotonPosition - fPMTPos[2];
      DotProduct = Displacement2.dot(fPhotonMomentum);
      closest = Displacement2 - (DotProduct * fPhotonMomentum);
      closest2 = (closest.x()*closest.x()+closest.y()*closest.y())/fDomeRadiusSq +
        (closest.z()*closest.z())/fDomeRadiuszSq;
      if (closest2 < 1. ) HitPMT += 2;
        
      Displacement3 = fPhotonPosition - fPMTPos[3];
      DotProduct = Displacement3.dot(fPhotonMomentum);
      closest = Displacement3 - (DotProduct * fPhotonMomentum);
      closest3 = (closest.x()*closest.x()+closest.y()*closest.y())/fDomeRadiusSq +
        (closest.z()*closest.z())/fDomeRadiuszSq;
      if (closest3 < 1. ) HitPMT += 4;
        
      switch (HitPMT) {
      case 0 :
        DistanceToPMT = 1.0e99;
        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;
      default :
        // hit all 3 pmts - error
        WARNING("Error determining next photon intersection in PropagateInTank()");
        return;
        break;
        //cerr << "Error determining next intersection of photon in" << endl;
        //cerr << "G4XTankFastCerenkov::PropogateInTank()" << endl;
        //exit(-1);
      }
        
      if (DistanceToPMT < DistanceToWall) {
        if (DistanceToPMT < DistanceToRoof) {
          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.0 - 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.0 - (rand / fRayleighFrac)) / fInvMFP;
        fPhotonPosition += DistanceTravelled * fPhotonMomentum;
        fPhotonTime += DistanceTravelled * fSampledRI / c_light;
        DoRayleighScatter(); // to pick new polarisation and direction
        if ((fPhotonMomentum.z() < 0.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;
        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.0)) 
      flag = DOWN_IN_TANK;
    else 
      flag = IN_TANK;
  
  } while (repeat);
  
  // never get here - we either get detected or absorbed inside the loop
  std::string err("Tank photon tracking ran into an impossible condition");
  ERROR(err);
  throw utl::OutOfBoundException(err);
}


// *************
// PMT ROUTINES
// *************

template<int pmtId> 
G4bool G4XTankFastCerenkov::TransitionToDome(G4int flag)
{
  fPhotonPosition -= fPMTPos[pmtId];
  G4ThreeVector Alpha, Beta, Delta, Normal;
  G4bool reflected, Dead/*, NormalFlag*/;
  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) ) {
    Dead = PropogateInDome<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.setX(2.*fPhotonPosition.x()/fDomeRadiusSq); 
    Normal.setY(2.*fPhotonPosition.y()/fDomeRadiusSq); 
    Normal.setZ(2.*fPhotonPosition.z()/fDomeRadiuszSq); 
    Normal /= Normal.mag();
    VDotN = -Normal.dot(fPhotonMomentum); // is actually cos(theta)
    if (VDotN > 0.999999) {
      Transmission = 4.0 * 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 = PropogateInDome<pmtId>(flag);
        if (Dead) return true;    
      }
      else {
        if (flag == OUTWARD) {
          fPhotonPosition += fPMTPos[pmtId];
          return false;
        }
        Dead = PropogateInDome<pmtId>(flag);
        if (Dead) return true;
      }
    }
    else {
      Alpha = Normal.cross(fPhotonMomentum);
      Alpha /= Alpha.mag();
      X = (n2 * n2) - (n * n) * (1.0 - (VDotN * VDotN));
      if (X <= 0.0) X = 0.0; // total reflection
      else {
        if (VDotN > 0.0) X = -sqrt(X);
        else 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.0 * n * VDotN * X) * (PerpComp * PerpComp +
                                               ParaComp * ParaComp);
      reflected = (G4UniformRand() >= Transmission);
      if (X == 0.0) reflected = true;
      if (flag == INWARD) {
        if (reflected) {
          // reflected back into tank
          fPhotonMomentum += 2.0 * VDotN * Normal;
          fPhotonMomentum /= fPhotonMomentum.mag();
          Beta = fPhotonMomentum.cross(Alpha);
          fPhotonPolarization = ParaComp * Alpha + PerpComp * Beta;
          fPhotonPolarization /= 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 /= fPhotonMomentum.mag();
          Beta = fPhotonMomentum.cross(Alpha);
          fPhotonPolarization = ParaComp * Alpha + PerpComp * Beta;
          fPhotonPolarization /= fPhotonPolarization.mag();
          Dead = PropogateInDome<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.0 * VDotN * Normal;
          fPhotonMomentum /= fPhotonMomentum.mag();
          Beta = fPhotonMomentum.cross(Alpha);
          fPhotonPolarization = ParaComp * Alpha + PerpComp * Beta;
          fPhotonPolarization /= fPhotonPolarization.mag();
          flag = INWARD;
          Dead = PropogateInDome<pmtId>(flag);
          if (Dead) return true;
        }
        else {
          // transmitted from dome back into tank
          fPhotonMomentum *= n;
          fPhotonMomentum += ((n * VDotN) + X) * Normal;
          fPhotonMomentum /= fPhotonMomentum.mag();
          Beta = fPhotonMomentum.cross(Alpha);
          fPhotonPolarization = ParaComp * Alpha + PerpComp * Beta;
          fPhotonPolarization /= fPhotonPolarization.mag();
          fPhotonPosition *= 1.00001; // Fudge to avoid rounding errors
          fPhotonPosition += fPMTPos[pmtId];
          return false;
        }
      }
    }
  } while (true);
}

template<int pmtId>
G4bool G4XTankFastCerenkov::PropogateInDome(G4int& flag)
{
  G4bool Dead;
  if (flag == IN_INTERFACE_1 ||
      flag == IN_INTERFACE_2 ||
      flag == IN_INTERFACE_3
      ) {
    Dead = TransitionToInterface<pmtId>(flag);
    if (Dead)  return true;
  }

  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.0) DistanceToRoof = 1.0e99;
    if (flag == INWARD) {
      //      DistanceThis = -2.0 * fPhotonPosition.dot(fPhotonMomentum);
      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 G4XTankFastCerenkov::TransitionToInterface(G4int& flag)
{
  G4ThreeVector Alpha, Beta, Delta, Normal;
  G4bool reflected, Dead/*, NormalFlag*/;
  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 
      ) {
    Dead = PropogateInInterface<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.setX(2.*fPhotonPosition.x()/fInterfaceRadiusSq); 
    Normal.setY(2.*fPhotonPosition.y()/fInterfaceRadiusSq); 
    Normal.setZ(2.*fPhotonPosition.z()/fInterfaceRadiuszSq); 
    Normal /= Normal.mag();
    VDotN = -Normal.dot(fPhotonMomentum); // is actually cos(theta)
    if (VDotN > 0.999999) {
      Transmission = 4.0 * 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 = PropogateInInterface<pmtId>(flag);
        if (Dead) return true;    
      }
      else {
        if (flag == OUTWARD) return true;
        Dead = PropogateInInterface<pmtId>(flag);
        if (Dead) return true;
      }
    }
    else {
      Alpha = Normal.cross(fPhotonMomentum);
      Alpha /= Alpha.mag();
      X = (n2 * n2) - (n * n) * (1.0 - (VDotN * VDotN));
      if (X <= 0.0) X = 0.0; // total reflection
      else {
        if (VDotN > 0.0) X = -sqrt(X);
        else 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.0 * n * VDotN * X) * (PerpComp * PerpComp +
                                               ParaComp * ParaComp);
      reflected = (G4UniformRand() >= Transmission);
      if (X == 0.0) reflected = true;
      if (flag == INWARD) {
        if (reflected) {
          // reflected back into dome
          fPhotonMomentum += 2.0 * VDotN * Normal;
          fPhotonMomentum /= fPhotonMomentum.mag();
          Beta = fPhotonMomentum.cross(Alpha);
          fPhotonPolarization = ParaComp * Alpha + PerpComp * Beta;
          fPhotonPolarization /= fPhotonPolarization.mag();
          flag = OUTWARD;
          return false;
        }
        else {
          // transmitted back into tank
          fPhotonMomentum *= n;
          fPhotonMomentum += ((n * VDotN) + X) * Normal;
          fPhotonMomentum /= fPhotonMomentum.mag();
          Beta = fPhotonMomentum.cross(Alpha);
          fPhotonPolarization = ParaComp * Alpha + PerpComp * Beta;
          fPhotonPolarization /= fPhotonPolarization.mag();
          Dead = PropogateInInterface<pmtId>(flag);
          if (Dead) return true;
        }
      }
      else {
        // reflected off inside of dome - remains in interface
        if (reflected) {
          fPhotonMomentum += 2.0 * VDotN * Normal;
          fPhotonMomentum /= fPhotonMomentum.mag();
          Beta = fPhotonMomentum.cross(Alpha);
          fPhotonPolarization = ParaComp * Alpha + PerpComp * Beta;
          fPhotonPolarization /= fPhotonPolarization.mag();
          flag = INWARD;
          Dead = PropogateInInterface<pmtId>(flag);
          if (Dead) return true;
        }
        else {
          // transmitted from interface back into dome
          fPhotonMomentum *= n;
          fPhotonMomentum += ((n * VDotN) + X) * Normal;
          fPhotonMomentum /= fPhotonMomentum.mag();
          Beta = fPhotonMomentum.cross(Alpha);
          fPhotonPolarization = ParaComp * Alpha + PerpComp * Beta;
          fPhotonPolarization /= fPhotonPolarization.mag();
          return false;
        }
      }
    }
  } while (true);
}

template<int pmtId>
G4bool G4XTankFastCerenkov::PropogateInInterface(G4int& flag)
{
  G4bool Dead;
  G4double DistanceThis, DistanceOther, DistanceToRoof;
  DistanceToRoof = (fRoofPos[pmtId] - fPhotonPosition.z()) / fPhotonMomentum.z();
  if (DistanceToRoof < 0.0) DistanceToRoof = 1.0e99;

  // we start at interface surface, so only one non-zero solution exists
  //  DistanceThis = -2.0 * fPhotonPosition.dot(fPhotonMomentum);
  DistanceThis = GetEllipsoidIntersect(fPhotonPosition, fInterfaceRadius, fInterfaceRadiusz);
  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 )
      G4XTankSimulator::fgPMTAction->AddPhoton(pmtId, fPhotonTime);
    return true; // photon absorbed, so is now dead
  }
}

// *******************
// Scattering routines
// *******************

G4bool G4XTankFastCerenkov::ScatterOffRoof(void)
{
  G4ThreeVector Normal(0.0, 0.0, -1.0);
  G4double Reflectivity = fLinerReflectivity.Y(fSampledMomentum);
  G4double r = G4UniformRand();
  if (r > Reflectivity) return true; // absorbed
  r /= Reflectivity; // range 0..1
  G4double Lobe = fLinerSpecularLobe.Y(fSampledMomentum);
  if (r <= Lobe) {
    LobeScatterHorizontal(Normal);
    return false;
  }
  G4double Spike = Lobe + fLinerSpecularSpike.Y(fSampledMomentum);
  if (r <= Spike) {
    SpikeScatter(Normal);
    return false;
  }
  G4double Back = Spike + fLinerBackscatter.Y(fSampledMomentum);
  if (r <= Back) {
    BackScatter();
    return false;
  }
  DiffuseScatterHorizontal(Normal);
  return false;
}

G4bool G4XTankFastCerenkov::ScatterOffWall(void)
{
  G4ThreeVector Normal;
  G4double Reflectivity = fLinerReflectivity.Y(fSampledMomentum);
  G4double r = G4UniformRand();
  if (r > Reflectivity) return true; // absorbed
  r /= Reflectivity; // range 0..1
  G4double X = fPhotonPosition.x();
  G4double Y = fPhotonPosition.y();
  G4double W = 1.0 / sqrt((X * X) + (Y * Y));
  Normal.set(-X * W, -Y * W, 0.0);
  G4double Lobe = fLinerSpecularLobe.Y(fSampledMomentum);
  if (r <= Lobe) {
    LobeScatterVertical(Normal);
    return false;
  }
  G4double Spike = Lobe + fLinerSpecularSpike.Y(fSampledMomentum);
  if (r <= Spike) {
    SpikeScatter(Normal);
    return false;
  }
  G4double Back = Spike + fLinerBackscatter.Y(fSampledMomentum);
  if (r <= Back) {
    BackScatter();
    return false;
  }
  DiffuseScatterVertical(Normal);
  return false;
}

G4bool G4XTankFastCerenkov::ScatterOffFloor(void)
{
  G4ThreeVector Normal(0.0, 0.0, 1.0);
  G4double Reflectivity = fLinerReflectivity.Y(fSampledMomentum);
  G4double r = G4UniformRand();
  if (r > Reflectivity) return true; // absorbed
  r /= Reflectivity; // range 0..1
  G4double Lobe = fLinerSpecularLobe.Y(fSampledMomentum);
  if (r <= Lobe) {
    LobeScatterHorizontal(Normal);
    return false;
  }
  G4double Spike = Lobe + fLinerSpecularSpike.Y(fSampledMomentum);
  if (r <= Spike) {
    SpikeScatter(Normal);
    return false;
  }
  G4double Back = Spike + fLinerBackscatter.Y(fSampledMomentum);
  if (r <= Back) {
    BackScatter();
    return false;
  }
  DiffuseScatterHorizontal(Normal);
  return false;
}

void G4XTankFastCerenkov::BackScatter(void)
{
  fPhotonMomentum.set(-fPhotonMomentum.x(), -fPhotonMomentum.y(),
                      -fPhotonMomentum.z());
  // can leave polarisation as it is
  return;
}

void G4XTankFastCerenkov::DiffuseScatterVertical(const G4ThreeVector& Normal)
{
  // for scattering off the tank walls
  G4double rand = G4UniformRand();
  G4double cosTheta = sqrt(1.0 - rand);
  G4double sinTheta = sqrt(rand);
  G4double Phi = 2.0 * utl::kPi * G4UniformRand();
  G4double sinThetaCosPhi = sinTheta * cos(Phi);
  G4double sinThetaSinPhi = sinTheta * sin(Phi);
  G4double X = Normal.x();
  G4double Y = Normal.y();
  G4double W = sqrt((X * X) + (Y * Y));
  G4double Vx = (X * cosTheta) - (Y * sinThetaSinPhi / W);
  G4double Vy = (Y * cosTheta) + (X * sinThetaSinPhi / W);
  G4double Vz = -W * sinThetaCosPhi;
  fPhotonMomentum.set(Vx, Vy, Vz);
  G4double DotProduct = fPhotonMomentum.dot(fPhotonPolarization);
  fPhotonPolarization = fPhotonPolarization - DotProduct * fPhotonMomentum;
  fPhotonPolarization /= fPhotonPolarization.mag();
  return;
}

void G4XTankFastCerenkov::DiffuseScatterHorizontal(const G4ThreeVector& Normal)
{
  // for scattering off tank floor or roof
  G4double rand = G4UniformRand();
  G4double cosTheta = sqrt(1.0 - rand);
  G4double sinTheta = sqrt(rand);
  G4double Phi = 2.0 * utl::kPi * G4UniformRand();
  G4double Vx = sinTheta * cos(Phi);
  G4double Vy = sinTheta * sin(Phi);
  G4double Vz = Normal.z() * cosTheta; // Normal.z() = +/- 1
  fPhotonMomentum.set(Vx, Vy, Vz);
  G4double DotProduct = fPhotonMomentum.dot(fPhotonPolarization);
  fPhotonPolarization = fPhotonPolarization - DotProduct * fPhotonMomentum;
  fPhotonPolarization /= fPhotonPolarization.mag();
  return;
}

void G4XTankFastCerenkov::SpikeScatter(const G4ThreeVector& Normal)
{
  G4double DotProduct = fPhotonMomentum.dot(Normal); // less than zero;
  fPhotonMomentum = fPhotonMomentum - (2.0 * DotProduct) * Normal;
  DotProduct = fPhotonMomentum.dot(fPhotonPolarization);
  fPhotonPolarization = fPhotonPolarization - DotProduct * fPhotonMomentum;
  fPhotonMomentum /= fPhotonMomentum.mag();
  fPhotonPolarization /= fPhotonPolarization.mag();
  return;
}

void G4XTankFastCerenkov::LobeScatterVertical(const G4ThreeVector& Normal)
{
  G4double Alpha, Phi;
  G4ThreeVector TrialMomentum;
  G4ThreeVector FacetNormal;
  G4double DotProduct;

  const G4double f_max = (fSigmaAlpha < 0.25) ? 4.0 * fSigmaAlpha : 1.0;

  G4double sinAlpha, cosAlpha, sinAlphaSinPhi, sinAlphaCosPhi;
  G4double W, X, Y, Vx, Vy, Vz;
  do {
    do {
      do {
        Alpha = G4RandGauss::shoot(0.0, fSigmaAlpha);
        sinAlpha = sin(Alpha);
      } while (((G4UniformRand()*f_max) > sinAlpha) || (Alpha >= (0.5 *utl::kPi)));
      Phi = 2.0 * utl::kPi * G4UniformRand();
      cosAlpha = cos(Alpha);
      sinAlphaSinPhi = sinAlpha * sin(Phi);
      sinAlphaCosPhi = sinAlpha * cos(Phi);
      X = Normal.x();
      Y = Normal.y();
      W = sqrt((X * X) + (Y * Y));
      Vx = (X * cosAlpha) - (Y * sinAlphaSinPhi / W);
      Vy = (Y * cosAlpha) + (X * sinAlphaSinPhi / W);
      Vz = - W * sinAlphaCosPhi;
      FacetNormal.set(Vx, Vy, Vz);
      DotProduct = fPhotonMomentum.dot(FacetNormal);
    } while (DotProduct >= 0.0);
    TrialMomentum = fPhotonMomentum - (2.0 * DotProduct) * FacetNormal;
  } while (TrialMomentum.dot(Normal) <= 0.0);
 
  fPhotonMomentum = TrialMomentum;
  DotProduct = fPhotonMomentum.dot(fPhotonPolarization);
  fPhotonPolarization = fPhotonPolarization - DotProduct * fPhotonMomentum;
  fPhotonMomentum /= fPhotonMomentum.mag();
  fPhotonPolarization /= fPhotonPolarization.mag();
  return;
}

void G4XTankFastCerenkov::LobeScatterHorizontal(const G4ThreeVector& Normal)
{
  G4double Alpha, Phi;
  G4ThreeVector TrialMomentum;
  G4ThreeVector FacetNormal;
  G4double DotProduct;
  G4double f_max = (fSigmaAlpha < 0.25) ? 4.0 * fSigmaAlpha : 1.0;
  G4double sinAlpha, cosAlpha/*, sinAlphaSinPhi, sinAlphaCosPhi*/;
  G4double /*W, X, Y, Z,*/ Vx, Vy, Vz;
  do {
    do {
      do {
        Alpha = G4RandGauss::shoot(0.0, fSigmaAlpha);
        sinAlpha = sin(Alpha);
      } while (((G4UniformRand()*f_max) > sinAlpha) || (Alpha >= (0.5 *utl::kPi)));
          
      Phi = 2.0 * utl::kPi * G4UniformRand();
      cosAlpha = cos(Alpha);
      Vx = sinAlpha * cos(Phi);
      Vy = sinAlpha * sin(Phi);
      Vz = Normal.z() * cosAlpha; // Normal.z() = +/- 1
      FacetNormal.set(Vx, Vy, Vz);
      DotProduct = fPhotonMomentum.dot(FacetNormal);
    } while (DotProduct >= 0.0);
      
    TrialMomentum = fPhotonMomentum - (2.0 * DotProduct) * FacetNormal;
  } while (TrialMomentum.dot(Normal) <= 0.0);
  
  fPhotonMomentum = TrialMomentum;
  DotProduct = fPhotonMomentum.dot(fPhotonPolarization);
  fPhotonPolarization = fPhotonPolarization - DotProduct * fPhotonMomentum;
  fPhotonMomentum /= fPhotonMomentum.mag();
  fPhotonPolarization /= fPhotonPolarization.mag();
  return;
}

// ****************
// Utility routines
// ****************

inline G4double G4XTankFastCerenkov::GetSphereIntersect(const G4ThreeVector& Pos,
                                                       const G4double r)
  const
{
  G4double b, c, d;
  G4double DistA, DistB;
    
  b = -Pos.dot(fPhotonMomentum);
  c = Pos.dot(Pos) - (r * r);
  d = (b * b) - c;
  if (d < 0.0) return 1.0e99; // no intersection with sphere at all
  d = sqrt(d);
  DistA = b + d;
  DistB = b - d;
  if (DistA > 0.0) {
    if (DistB < DistA) {
      if(DistB > 0.0) return DistB;
      return DistA; // inside sphere - B is behind us, A ahead
    }
    return DistA;
  }
  else {
    if (DistB > 0.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 1.0e99;
}

inline G4double G4XTankFastCerenkov::GetEllipsoidIntersect(const G4ThreeVector& Pos,
                                                          const G4double r1, const G4double r2)
  const
{
  G4double a1, a2, b[3], c[3], d[3], e, f, g;
  //  G4double DistA, DistB;

  // 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
    
  b[0] = fPhotonMomentum.x()*fPhotonMomentum.x();
  b[1] = fPhotonMomentum.y()*fPhotonMomentum.y();
  b[2] = fPhotonMomentum.z()*fPhotonMomentum.z();
  c[0] = Pos.x()*Pos.x();
  c[1] = Pos.y()*Pos.y();
  c[2] = Pos.z()*Pos.z();
  d[0] = fPhotonMomentum.x()*Pos.x();
  d[1] = fPhotonMomentum.y()*Pos.y();
  d[2] = fPhotonMomentum.z()*Pos.z();

  e = (b[0]+b[1])/(r1*r1) + b[2]/(r2*r2);
  f = 2.*((d[0]+d[1])/(r1*r1) + d[2]/(r2*r2));
  g = (c[0]+c[1])/(r1*r1)+c[2]/(r2*r2) - 1.;

  // if no solution of the quadratic equation there is no intersection
  if ( (f*f-4*e*g) < 0 ) {
    return 1.0e99;
  }
  else { // take smallest positive value of the solution
    a1 = 1./(2*e) * ( -f + sqrt(f*f-4*e*g) );
    a2 = 1./(2*e) * ( -f - sqrt(f*f-4*e*g) );
    if ( a2 > 10e-5 ) // to avoid almost 0 but positive a2 for current photon position
      return a2;
    else if ( a1 > 0 )
      return a1;
  }
  
  return 1.0e99;
    
}


//not used
// G4double G4XTankFastCerenkov::RandGaussZero(const G4double Sigma)
// {
//   // generates Gaussian random numbers with given Sigma and mean of 0
//   static G4int i;
//   static G4double y1, y2;
//   G4double x1, x2, w;
  
//   i = 1 - i;
//   if (i == 0) return Sigma * y2;
//   do {
//     x1 = 2.0 * G4UniformRand() - 1.0;
//     x2 = 2.0 * G4UniformRand() - 1.0;
//     w = x1 * x1 + x2 * x2;
//   } while (w >= 1.0);
//   w = sqrt((-2.0 * log(w)) / w);
//   y1 = x1 * w;
//   y2 = x2 * w;
//   return Sigma * y1;
// }