Discrete process for reflection and diffusion at optical interfaces. More...

#include <OpBoundaryProcessSimple.hh>

Inheritance diagram for larg4::OpBoundaryProcessSimple:

Public Member Functions
	OpBoundaryProcessSimple (const G4String &processName="OpBoundary", G4ProcessType type=fOptical)

G4bool	IsApplicable (const G4ParticleDefinition &aParticleType)

G4double	GetMeanFreePath (const G4Track &, G4double, G4ForceCondition *condition)

G4VParticleChange *	PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)

OpBoundaryProcessSimpleStatus	GetStatus () const

Private Member Functions
G4bool	G4BooleanRand (const G4double prob) const

Private Attributes
OpBoundaryProcessSimpleStatus	fTheStatus

G4double	fCarTolerance

int	fVerbosity

Detailed Description

Discrete process for reflection and diffusion at optical interfaces.

See also: G4VDiscreteProcess

This class invokes a simplified model of optical reflections at boundaries between different materials. The relevant reflectivities are ultimately read from detinfo::LArProperties via larg4::MaterialPropertiesLoader.

The required parameters are total reflectance (detinfo::LArProperties::SurfaceReflectances()) and ratio of diffuse to specular reflectance (detinfo::LArProperties::SurfaceReflectanceDiffuseFractions()). Each photon crossing a boundary with a defined reflectance is randomly either reflected or absorbed and killed according to the supplied probability.

Every reflected photon with a defined diffuse reflection fraction is then randomly either diffusely or specularly reflected according to the supplied probability. All materials with no defined reflectance are assumed to be black and absorb all incident photons.

This physics process is loaded in larg4::OpticalPhysics physics constructor.

This class is based on the G4OpBoundaryProcess class in Geant4 and was adapted for LArSoft by Ben Jones, MIT, March 2010.

Definition at line 102 of file OpBoundaryProcessSimple.hh.

Constructor & Destructor Documentation

larg4::OpBoundaryProcessSimple::OpBoundaryProcessSimple	(	const G4String &	processName = `"OpBoundary"`,
		G4ProcessType	type = `fOptical`
	)

Definition at line 65 of file OpBoundaryProcessSimple.cxx.

     : G4VDiscreteProcess(processName, type)
   {
     G4cout << GetProcessName() << " is created " << G4endl;
 
     SetProcessSubType(fOpBoundary);
 
     fTheStatus = Undefined;
 
     fCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
   }

Member Function Documentation

G4bool larg4::OpBoundaryProcessSimple::G4BooleanRand ( const G4double prob ) const

inlineprivate

Definition at line 133 of file OpBoundaryProcessSimple.hh.

   {
     /* Returns a random boolean variable with the specified probability */
     return (G4UniformRand() < prob);
   }

G4double larg4::OpBoundaryProcessSimple::GetMeanFreePath	(	const G4Track &	,
		G4double	,
		G4ForceCondition *	condition
	)

Definition at line 233 of file OpBoundaryProcessSimple.cxx.

   {
     *condition = Forced;
 
     return DBL_MAX;
   }

OpBoundaryProcessSimpleStatus larg4::OpBoundaryProcessSimple::GetStatus ( void ) const

inline

Definition at line 146 of file OpBoundaryProcessSimple.hh.

   {
     return fTheStatus;
   }

G4bool larg4::OpBoundaryProcessSimple::IsApplicable ( const G4ParticleDefinition & aParticleType )

inline

Definition at line 140 of file OpBoundaryProcessSimple.hh.

   {
     return (&aParticleType == G4OpticalPhoton::OpticalPhoton());
   }

G4VParticleChange * larg4::OpBoundaryProcessSimple::PostStepDoIt	(	const G4Track &	aTrack,
		const G4Step &	aStep
	)

Definition at line 81 of file OpBoundaryProcessSimple.cxx.

   {
 
     // Note - these have to be loaded here, since this object
     // is constructed before LArG4, and hence before the
     // LArG4 parameters are read from config
 
     fTheStatus = Undefined;
     aParticleChange.Initialize(aTrack);
 
     G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
     G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
 
     if (pPostStepPoint->GetStepStatus() != fGeomBoundary) {
       fTheStatus = NotAtBoundary;
       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
     }
     if (aTrack.GetStepLength() <= fCarTolerance / 2) {
       fTheStatus = StepTooSmall;
       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
     }
 
     G4Material* Material1 = pPreStepPoint->GetMaterial();
     G4Material* Material2 = pPostStepPoint->GetMaterial();
 
     if (Material1 != Material2) {
 
       art::ServiceHandle<sim::LArG4Parameters const> lgp;
 
       fVerbosity = lgp->OpVerbosity();
 
       G4ThreeVector NewMomentum;
       G4ThreeVector NewPolarization;
 
       G4ThreeVector theGlobalNormal;
       G4ThreeVector theFacetNormal;
 
       std::string Material1Name = Material1->GetName();
       std::string Material2Name = Material2->GetName();
 
       const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
       G4double thePhotonMomentum = aParticle->GetTotalMomentum();
       G4ThreeVector OldMomentum = aParticle->GetMomentumDirection();
       G4ThreeVector OldPolarization = aParticle->GetPolarization();
 
       if (fVerbosity > 9)
         std::cout << "OpBoundaryProcessSimple Debug: Photon " << aTrack.GetTrackID() << " momentum "
                   << thePhotonMomentum << " Material1 " << Material1Name.c_str() << " Material 2 "
                   << Material2Name.c_str() << std::endl;
 
       G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();
 
       G4Navigator* theNavigator =
         G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
 
       G4ThreeVector theLocalPoint =
         theNavigator->GetGlobalToLocalTransform().TransformPoint(theGlobalPoint);
 
       // Normal points back into volume
       G4bool valid;
       G4ThreeVector theLocalNormal = theNavigator->GetLocalExitNormal(&valid);
 
       if (valid) { theLocalNormal = -theLocalNormal; }
       else {
         G4cerr << " OpBoundaryProcessSimple/PostStepDoIt(): "
                << " The Navigator reports that it returned an invalid normal" << G4endl;
       }
 
       theGlobalNormal = theNavigator->GetLocalToGlobalTransform().TransformAxis(theLocalNormal);
 
       if (OldMomentum * theGlobalNormal > 0.0) {
 #ifdef G4DEBUG_OPTICAL
         G4cerr << " OpBoundaryProcessSimple/PostStepDoIt(): "
                << " theGlobalNormal points the wrong direction " << G4endl;
 #endif
         theGlobalNormal = -theGlobalNormal;
       }
 
       G4MaterialPropertiesTable* aMaterialPropertiesTable;
       G4MaterialPropertyVector* Reflectance;
       G4MaterialPropertyVector* DiffuseReflectanceFraction;
 
       aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
 
       if (aMaterialPropertiesTable) {
 
         std::stringstream PropertyName;
         PropertyName << "REFLECTANCE_" << Material2Name.c_str();
         Reflectance = aMaterialPropertiesTable->GetProperty(PropertyName.str().c_str());
 
         PropertyName.str("");
         PropertyName << "DIFFUSE_REFLECTANCE_FRACTION_" << Material2Name.c_str();
         DiffuseReflectanceFraction =
           aMaterialPropertiesTable->GetProperty(PropertyName.str().c_str());
 
         if (Reflectance) {
           double theReflectance = Reflectance->Value(thePhotonMomentum);
           if (G4BooleanRand(theReflectance)) {
             if (DiffuseReflectanceFraction) {
               double theDiffuseReflectanceFraction =
                 DiffuseReflectanceFraction->Value(thePhotonMomentum);
               if (G4BooleanRand(theDiffuseReflectanceFraction)) { fTheStatus = SimpleDiffuse; }
               else {
                 fTheStatus = SimpleSpecular;
               }
             }
             else {
               fTheStatus = SimpleSpecular;
             }
           }
           else {
             fTheStatus = SimpleAbsorbed;
           }
         }
         else {
           fTheStatus = SimpleAbsorbedNoRefl;
         }
       }
 
       // Take action according to track status
 
       if (fTheStatus == SimpleDiffuse) {
         NewMomentum = G4LambertianRand(theGlobalNormal);
         theFacetNormal = (NewMomentum - OldMomentum).unit();
       }
 
       else if (fTheStatus == SimpleSpecular) {
         theFacetNormal = theGlobalNormal;
         G4double PdotN = OldMomentum * theFacetNormal;
         NewMomentum = OldMomentum - (2. * PdotN) * theFacetNormal;
       }
 
       else if ((fTheStatus == SimpleAbsorbed) || (fTheStatus == SimpleAbsorbedNoRefl)) {
         aParticleChange.ProposeTrackStatus(fStopAndKill);
       }
 
       else if ((fTheStatus == NotAtBoundary) || (fTheStatus == StepTooSmall)) {
         NewMomentum = OldMomentum;
       }
       else
         aParticleChange.ProposeTrackStatus(fStopAndKill);
 
       NewMomentum = NewMomentum.unit();
       aParticleChange.ProposeMomentumDirection(NewMomentum);
     }
     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
   }