G4StationSimulator/G4TankOpBoundaryProcess.cc

Go to the documentation of this file.

#include <config.h>

// Optical Photon Boundary Process Class Implementation
//
// File:        G4TankOpBoundaryProcess.cc
// Description: Discrete Process -- reflection/refraction at
//                                  optical interfaces
// Version:     1.1
// Created:     1997-06-18
// Modified:    1998-05-25 - Correct parallel component of polarization
//                           (thanks to: Stefano Magni + Giovanni Pieri)
//              1998-05-28 - NULL Rindex pointer before reuse
//                           (thanks to: Stefano Magni)
//              1998-06-11 - delete *sint1 in oblique reflection
//                           (thanks to: Giovanni Pieri)
//              1998-06-19 - move from GetLocalExitNormal() to the new
//                           method: GetLocalExitNormal(&valid) to get
//                           the surface normal in all cases
//              1998-11-07 - NULL OpticalSurface pointer before use
//                           comparison not sharp for: abs(cost1) < 1.0
//                           remove sin1, sin2 in lines 556,567
//                           (thanks to Stefano Magni)
//              1999-10-10 - Accommodate changes done in DoAbsorption by
//                           changing logic in DielectricMetal
//
// Author:      Peter Gumplinger
//              adopted from work by Werner Keil - April 2/96
// mail:        gum@triumf.ca
//

#include "G4ios.hh"
#include "G4ProcessType.hh"

#include "G4TankOpBoundaryProcess.h"

#include <utl/Math.h>
#include <utl/config.h>

#include <iostream>

using namespace std;
using utl::Sqr;


namespace G4StationSimulatorOG {

  G4TankOpBoundaryProcess::G4TankOpBoundaryProcess(const G4String& processName,
                                                   const G4ProcessType /*type*/) :
    G4VDiscreteProcess(processName, fOptical)
  {
    if (verboseLevel > 0)
      G4cerr << GetProcessName() << " is created " << G4endl;
  }


  G4VParticleChange*
  G4TankOpBoundaryProcess::PostStepDoIt(const G4Track& track, const G4Step& step)
  {
    aParticleChange.Initialize(track);

    G4StepPoint* const preStepPoint = step.GetPreStepPoint();
    G4StepPoint* const postStepPoint = step.GetPostStepPoint();

    if (postStepPoint->GetStepStatus() != fGeomBoundary)
      return G4VDiscreteProcess::PostStepDoIt(track, step);

    if (track.GetStepLength() <= kCarTolerance / 2)
      return G4VDiscreteProcess::PostStepDoIt(track, step);

    fMaterial1 = preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial();
    fMaterial2 = postStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial();

    if (fMaterial1 == fMaterial2)
      return G4VDiscreteProcess::PostStepDoIt(track, step);

    const G4DynamicParticle* const particle = track.GetDynamicParticle();

    fPhotonMomentum = particle->GetTotalMomentum();
    fOldMomentum = particle->GetMomentumDirection();
    fOldPolarization = particle->GetPolarization();

    if (verboseLevel > 0) {
      G4cerr << " Photon at Boundary!\n"
                " Old Momentum Direction: " << fOldMomentum << "\n"
                " Old Polarization:       " << fOldPolarization << G4endl;
    }

    {
      const auto t = fMaterial1->GetMaterialPropertiesTable();
      const auto ri = t ? t->GetProperty("RINDEX") : nullptr;
      if (ri)
        fRIndex1 = ri->Value(fPhotonMomentum);
      else {
        fRIndex1 = 0;
        aParticleChange.ProposeTrackStatus(fStopAndKill);
        return G4VDiscreteProcess::PostStepDoIt(track, step);
      }
    }

    if (const auto s = G4LogicalBorderSurface::GetSurface(preStepPoint->GetPhysicalVolume(), postStepPoint->GetPhysicalVolume()))
      fOpticalSurface = (G4OpticalSurface*)s->GetSurfaceProperty();
    else if (const auto s = G4LogicalSkinSurface::GetSurface(preStepPoint->GetPhysicalVolume()->GetLogicalVolume()))
      fOpticalSurface = (G4OpticalSurface*)s->GetSurfaceProperty();
    else
      fOpticalSurface = nullptr;

    fModel = glisur;
    fFinish = polished;

    G4SurfaceType type = dielectric_dielectric;
    fRIndex2 = 0;
    {
      const auto t = fMaterial2->GetMaterialPropertiesTable();
      const auto ri = t ? t->GetProperty("RINDEX") : nullptr;
      if (ri) {
        fRIndex2 = ri->Value(fPhotonMomentum);
      } else if (fOpticalSurface) {
        type = fOpticalSurface->GetType();
      } else {
        fRIndex2 = 0;
        aParticleChange.ProposeTrackStatus(fStopAndKill);
        return G4VDiscreteProcess::PostStepDoIt(track, step);
      }
    }

    if (fOpticalSurface) {
      fModel = fOpticalSurface->GetModel();
      fFinish = fOpticalSurface->GetFinish();
      const auto t = fOpticalSurface->GetMaterialPropertiesTable();
      if (t) {
        if (!fRIndex2) {
          if (const auto p = t->GetProperty("RINDEX"))
            fRIndex2 = p->Value(fPhotonMomentum);
        }
        if (const auto p = t->GetProperty("REFLECTIVITY"))
          fReflectivity = p->Value(fPhotonMomentum);
        if (const auto p = t->GetProperty("EFFICIENCY"))
          fEfficiency = p->Value(fPhotonMomentum);
        if (fModel == unified) {
          if (const auto p = t->GetProperty("SPECULARLOBECONSTANT"))
            fProb_sl = p->Value(fPhotonMomentum);
          if (const auto p = t->GetProperty("SPECULARSPIKECONSTANT"))
            fProb_ss = p->Value(fPhotonMomentum);
          if (const auto p = t->GetProperty("BACKSCATTERCONSTANT"))
            fProb_bs = p->Value(fPhotonMomentum);
        }
      }
    }

    // No optical surface

    const G4ThreeVector& globalPoint = postStepPoint->GetPosition();

    G4Navigator* const navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();

    const G4ThreeVector localPoint = navigator->GetGlobalToLocalTransform().TransformPoint(globalPoint);

    G4bool valid = false;
    // Normal points back into volume
    const G4ThreeVector localNormal = -navigator->GetLocalExitNormal(&valid);
    if (!valid) {
      G4cerr << " G4TankOpBoundaryProcess/PostStepDoIt(): "
                " The Navigator reports that it returned an invalid normal"
             << G4endl;
    }

    fGlobalNormal = navigator->GetLocalToGlobalTransform().TransformAxis(localNormal);
    fStatus = Undefined;

    if (type == dielectric_metal)

      DielectricMetal();

    else if (type == dielectric_dielectric) {

      if (fFinish == polishedfrontpainted || fFinish == groundfrontpainted) {

        if (!G4BooleanRand(fReflectivity)) {

          DoAbsorption();

        } else {

          if (fFinish == groundfrontpainted)
            fStatus = LambertianReflection;
          DoReflection();

        }

      } else {

        // MY HACK (tp) (I added the DoAbsorption option so that it takes
        // effect even for the ground finish.  I had to do this because
        // groundbackpainted option causes a crash somewhere in G4Tubs..)
        // tm : Oct. 19, 2001
        // This hack creates a problem when the photons go directly from the
        // water to a surface where the OpticalSurface machinery has not been
        // used (like the dome) The reflectivity is then 0 and all the photon
        // are absorbed. Also, this hack screws up a normal optical interface
        // because it checks against the reflectivity instead of doing it the
        // way it should, in for example, Jackson. Only check against the
        // reflectivity for tyvek...

        if (fMaterial2->GetName() == "HDPE") {

          if (!G4BooleanRand(fReflectivity))
            DoAbsorption();
          else
            DielectricDielectric();

        } else {

          DielectricDielectric();

        }

      }

    } else {

      G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
      return G4VDiscreteProcess::PostStepDoIt(track, step);

    }

    fNewMomentum = fNewMomentum.unit();
    fNewPolarization = fNewPolarization.unit();

    aParticleChange.ProposeMomentumDirection(fNewMomentum);
    aParticleChange.ProposePolarization(fNewPolarization);

    return G4VDiscreteProcess::PostStepDoIt(track, step);
  }


  G4ThreeVector
  G4TankOpBoundaryProcess::G4IsotropicRand()
    const
  {
    // Returns a random isotropic unit vector
    G4ThreeVector vect;
    G4double len2;
    do {
      vect.set(G4UniformRand() - 0.5, G4UniformRand() - 0.5, G4UniformRand() - 0.5);
      len2 = vect.mag2();
    } while (len2 < 0.01 || len2 > 0.25);
    return vect.unit();
  }


  G4ThreeVector
  G4TankOpBoundaryProcess::G4LambertianRand(const G4ThreeVector& normal)
    const
  {
    // Returns a random lambertian unit vector
    G4ThreeVector vect;
    G4double ndotv;
    do {
      vect = G4IsotropicRand();
      ndotv = normal * vect;
      if (ndotv < 0) {
        vect *= -1;
        ndotv *= -1;
      }
    } while (!G4BooleanRand(ndotv));
    return vect;
  }


  G4ThreeVector
  G4TankOpBoundaryProcess::G4PlaneVectorRand(const G4ThreeVector& normal)
    const
  {
    // This function chooses a random vector within a plane given by the unit normal
    const G4ThreeVector vec1 = normal.orthogonal();
    const G4ThreeVector vec2 = vec1.cross(normal);
    const G4double phi = twopi * G4UniformRand();
    const G4double cosphi = std::cos(phi);
    const G4double sinphi = std::sin(phi);
    return cosphi * vec1 + sinphi * vec2;
  }


  G4ThreeVector
  G4TankOpBoundaryProcess::GetFacetNormal(const G4ThreeVector& momentum, const G4ThreeVector& normal)
    const
  {
    G4ThreeVector facetNormal;

    if (fModel == unified) {

      /* This function code alpha to a random value taken from the
       distribution p(alpha) = g(alpha; 0, sigma_alpha)*sin(alpha),
       for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
       is a gaussian distribution with mean 0 and standard deviation
       sigma_alpha.  */

      G4double sigmaAlpha = 0;

      if (fOpticalSurface)
        sigmaAlpha = fOpticalSurface->GetSigmaAlpha();

      const G4double maxF = std::min(1., 4 * sigmaAlpha);

      do {

        G4double alpha;
        do {
          alpha = G4RandGauss::shoot(0, sigmaAlpha);
        } while (G4UniformRand() * maxF > std::sin(alpha) || alpha >= halfpi);
        const G4double phi = G4UniformRand() * twopi;

        const G4double sinAlpha = std::sin(alpha);
        const G4double cosAlpha = std::cos(alpha);
        const G4double sinPhi = std::sin(phi);
        const G4double cosPhi = std::cos(phi);

        facetNormal.set(sinAlpha*cosPhi, sinAlpha*sinPhi, cosAlpha);
        facetNormal.rotateUz(normal);

      } while (momentum * facetNormal >= 0);

    } else {

      G4double polish = 1;

      if (fOpticalSurface)
        polish = fOpticalSurface->GetPolish();

      if (polish >= 1)
        facetNormal = normal;
      else {
        do {
          G4ThreeVector smear;
          do {
            smear.set(
              2 * G4UniformRand() - 1,
              2 * G4UniformRand() - 1,
              2 * G4UniformRand() - 1
            );
          } while (smear.mag() > 1);
          smear = (1 - polish) * smear;
          facetNormal = normal + smear;
        } while (momentum * facetNormal >= 0);
        facetNormal = facetNormal.unit();

      }

    }

    return facetNormal;
  }


  void
  G4TankOpBoundaryProcess::DielectricMetal()
  {
    do {
      if (!G4BooleanRand(fReflectivity)) {
        DoAbsorption();
        break;
      } else {
        DoReflection();
        fOldMomentum = fNewMomentum;
        fOldPolarization = fNewPolarization;
      }
    } while (fNewMomentum * fGlobalNormal < 0);
  }


  void
  G4TankOpBoundaryProcess::DielectricDielectric()
  {
    G4bool inside = false;
    G4bool swap = false;
    G4bool done;

    do {

      G4bool through = false;
      done = false;

      do {

        if (through) {
          swap = !swap;
          through = false;
          fGlobalNormal = -fGlobalNormal;
          std::swap(fMaterial1, fMaterial2);
          std::swap(fRIndex1, fRIndex2);
        }

        if (fFinish == ground || fFinish == groundbackpainted) {
          fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
        } else {
          fFacetNormal = fGlobalNormal;
        }

        G4double cost1 = -(fOldMomentum * fFacetNormal);
        G4double sint1 = 0;
        G4double sint2 = 0;

        if (std::abs(cost1) < 1 - kCarTolerance) {
          sint1 = std::sqrt(1 - Sqr(cost1));
          sint2 = sint1 * fRIndex1 / fRIndex2;  // Snell's Law
        }

        if (sint2 >= 1) {

          // Simulate total internal reflection
          if (swap)
            swap = false;

          fStatus = TotalInternalReflection;

          if (fModel == unified && fFinish != polished)
            ChooseReflection();

          if (fStatus == LambertianReflection) {
            DoReflection();
          } else if (fStatus == BackScattering) {
            fNewMomentum = -fOldMomentum;
            fNewPolarization = -fOldPolarization;
          } else {
            const G4double pDotN = fOldMomentum * fFacetNormal;
            fNewMomentum = fOldMomentum - (2 * pDotN) * fFacetNormal;
            const G4double eDotN = fOldPolarization * fFacetNormal;
            fNewPolarization = -fOldPolarization + (2 * eDotN) * fFacetNormal;
          }

        } else if (sint2 < 1) {

          // Calculate amplitude for transmission (Q = P x N)
          const G4double cost2 = (cost1 > 0) ? std::sqrt(1 - Sqr(sint2)) : -std::sqrt(1 - Sqr(sint2));

          G4ThreeVector a_trans;
          G4ThreeVector atrans;
          G4ThreeVector e1pp;
          G4ThreeVector e1pl;
          G4double e1_perp;
          G4double e1_parl;

          if (sint1 > 0) {
            a_trans = fOldMomentum.cross(fFacetNormal);
            atrans = a_trans.unit();
            e1_perp = fOldPolarization * atrans;
            e1pp = e1_perp * atrans;
            e1pl = fOldPolarization - e1pp;
            e1_parl = e1pl.mag();
          } else {
            a_trans = fOldPolarization;
            // Here we Follow Jackson's conventions and we set the parallel
            // component = 1 in case of a ray perpendicular to the surface
            e1_perp = 0;
            e1_parl = 1;
          }

          G4double e2_perp = 0;
          G4double e2_parl = 0;
          G4double e2_total = 0;
          G4double transCoeff = 0;

          if (cost1 != 0) {
            const G4double s1 = fRIndex1 * cost1;
            e2_perp = 2 * s1 * e1_perp / (fRIndex1 * cost1 + fRIndex2 * cost2);
            e2_parl = 2 * s1 * e1_parl / (fRIndex2 * cost1 + fRIndex1 * cost2);
            e2_total = e2_perp * e2_perp + e2_parl * e2_parl;
            const G4double s2 = fRIndex2 * cost2 * e2_total;
            transCoeff = s2 / s1;
          }

          G4ThreeVector refracted;
          G4ThreeVector deflected;
          G4double e2_abs;
          G4double c_parl;
          G4double c_perp;

          if (!G4BooleanRand(transCoeff)) {

            // Simulate reflection
            if (swap)
              swap = false;

            fStatus = FresnelReflection;

            if (fModel == unified && fFinish != polished)
              ChooseReflection();

            if (fStatus == LambertianReflection) {

              DoReflection();

            } else if (fStatus == BackScattering) {

              fNewMomentum = -fOldMomentum;
              fNewPolarization = -fOldPolarization;

            } else {

              const G4double pDotN = fOldMomentum * fFacetNormal;
              fNewMomentum = fOldMomentum - (2 * pDotN) * fFacetNormal;

              if (sint1 > 0) {  // incident ray oblique

                e2_parl = fRIndex2 * e2_parl / fRIndex1 - e1_parl;
                e2_perp = e2_perp - e1_perp;
                e2_total = e2_perp * e2_perp + e2_parl * e2_parl;
                refracted = fFacetNormal + pDotN * fNewMomentum;
                e2_abs = std::sqrt(e2_total);
                c_parl = e2_parl / e2_abs;
                c_perp = e2_perp / e2_abs;
                fNewPolarization = c_parl * refracted - c_perp * a_trans;

              } else if (fRIndex2 > fRIndex1) {  // incident ray perpendicular

                fNewPolarization = -fOldPolarization;

              }

            }

          } else {  // photon gets transmitted

            // Simulate transmission/refraction
            inside = !inside;
            through = true;
            fStatus = FresnelRefraction;

            if (sint1 > 0) {  // incident ray oblique

              const G4double alpha = cost1 - cost2 * fRIndex2 / fRIndex1;
              deflected = fOldMomentum + alpha * fFacetNormal;
              fNewMomentum = deflected.unit();
              const G4double pDotN = -cost2;
              refracted = fFacetNormal - pDotN * fNewMomentum;
              e2_abs = std::sqrt(e2_total);
              c_parl = e2_parl / e2_abs;
              c_perp = e2_perp / e2_abs;
              fNewPolarization = c_parl * refracted + c_perp * a_trans;

            } else {  // incident ray perpendicular

              fNewMomentum = fOldMomentum;
              fNewPolarization = fOldPolarization;

            }

          }

        }

        fOldMomentum = fNewMomentum;
        fOldPolarization = fNewPolarization;

        {
          const G4double proj = fNewMomentum * fGlobalNormal;
          done = (fStatus == FresnelRefraction) ? (proj <= 0) : (proj >= 0);
        }

      } while (!done);

      if (inside &&
          !swap &&
          (fFinish == polishedbackpainted || fFinish == groundbackpainted)) {

        if (!G4BooleanRand(fReflectivity))
          DoAbsorption();
        else {
          if (fStatus != FresnelRefraction)
            fGlobalNormal = -fGlobalNormal;
          else {
            swap = !swap;
            std::swap(fMaterial1, fMaterial2);
            std::swap(fRIndex1, fRIndex2);
          }

          if (fFinish == groundbackpainted)
            fStatus = LambertianReflection;

          DoReflection();

          fGlobalNormal = -fGlobalNormal;
          fOldMomentum = fNewMomentum;

          done = false;  // go again
        }

      }

    } while (!done);
  }


  void
  G4TankOpBoundaryProcess::ChooseReflection()
  {
    const G4double rand = G4UniformRand();
    if (rand >= 0 && rand < fProb_ss) {
      fStatus = SpikeReflection;
      fFacetNormal = fGlobalNormal;
    } else if (rand >= fProb_ss && rand <= fProb_ss + fProb_sl)
      fStatus = LobeReflection;
    else if (rand > fProb_ss + fProb_sl && rand < fProb_ss + fProb_sl + fProb_bs)
      fStatus = BackScattering;
    else
      fStatus = LambertianReflection;
  }


  void
  G4TankOpBoundaryProcess::DoAbsorption()
  {
    fStatus = Absorption;

    if (G4BooleanRand(fEfficiency)) {
      // EnergyDeposited =/= 0 means: photon has been detected
      fStatus = Detection;
      aParticleChange.ProposeLocalEnergyDeposit(fPhotonMomentum);
    } else
      aParticleChange.ProposeLocalEnergyDeposit(0);

    fNewMomentum = fOldMomentum;
    fNewPolarization = fOldPolarization;

    aParticleChange.ProposeTrackStatus(fStopAndKill);
  }


  void
  G4TankOpBoundaryProcess::DoReflection()
  {
    if (fStatus == LambertianReflection) {
      fNewMomentum = G4LambertianRand(fGlobalNormal);
      fFacetNormal = (fNewMomentum - fOldMomentum).unit();
    } else if (fFinish == ground) {
      fStatus = LobeReflection;
      fFacetNormal = GetFacetNormal(fOldMomentum, fGlobalNormal);
      const G4double pDotN = fOldMomentum * fFacetNormal;
      fNewMomentum = fOldMomentum - (2 * pDotN) * fFacetNormal;
    } else {
      fStatus = SpikeReflection;
      fFacetNormal = fGlobalNormal;
      const G4double pDotN = fOldMomentum * fFacetNormal;
      fNewMomentum = fOldMomentum - (2 * pDotN) * fFacetNormal;
    }
    const G4double eDotN = fOldPolarization * fFacetNormal;
    fNewPolarization = -fOldPolarization + (2 * eDotN) * fFacetNormal;
  }

}