Geant4  10.02.p03
G4OpRayleigh Class Reference

#include <G4OpRayleigh.hh>

Inheritance diagram for G4OpRayleigh:
Collaboration diagram for G4OpRayleigh:

Public Member Functions

 G4OpRayleigh (const G4String &processName="OpRayleigh", G4ProcessType type=fOptical)
 
 ~G4OpRayleigh ()
 
G4bool IsApplicable (const G4ParticleDefinition &aParticleType)
 
void BuildPhysicsTable (const G4ParticleDefinition &aParticleType)
 
G4double GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *)
 
G4VParticleChange * PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
 
G4PhysicsTableGetPhysicsTable () const
 
void DumpPhysicsTable () const
 
- Public Member Functions inherited from G4VDiscreteProcess
 G4VDiscreteProcess (const G4String &, G4ProcessType aType=fNotDefined)
 
 G4VDiscreteProcess (G4VDiscreteProcess &)
 
virtual ~G4VDiscreteProcess ()
 
virtual G4double PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
 
virtual G4double AlongStepGetPhysicalInteractionLength (const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
 
virtual G4double AtRestGetPhysicalInteractionLength (const G4Track &, G4ForceCondition *)
 
virtual G4VParticleChange * AtRestDoIt (const G4Track &, const G4Step &)
 
virtual G4VParticleChange * AlongStepDoIt (const G4Track &, const G4Step &)
 
- Public Member Functions inherited from G4VProcess
 G4VProcess (const G4String &aName="NoName", G4ProcessType aType=fNotDefined)
 
 G4VProcess (const G4VProcess &right)
 
virtual ~G4VProcess ()
 
G4int operator== (const G4VProcess &right) const
 
G4int operator!= (const G4VProcess &right) const
 
G4double GetCurrentInteractionLength () const
 
void SetPILfactor (G4double value)
 
G4double GetPILfactor () const
 
G4double AlongStepGPIL (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)
 
G4double AtRestGPIL (const G4Track &track, G4ForceCondition *condition)
 
G4double PostStepGPIL (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
 
virtual void PreparePhysicsTable (const G4ParticleDefinition &)
 
virtual G4bool StorePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
 
virtual G4bool RetrievePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
 
const G4StringGetPhysicsTableFileName (const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
 
const G4StringGetProcessName () const
 
G4ProcessType GetProcessType () const
 
void SetProcessType (G4ProcessType)
 
G4int GetProcessSubType () const
 
void SetProcessSubType (G4int)
 
virtual void StartTracking (G4Track *)
 
virtual void EndTracking ()
 
virtual void SetProcessManager (const G4ProcessManager *)
 
virtual const G4ProcessManagerGetProcessManager ()
 
virtual void ResetNumberOfInteractionLengthLeft ()
 
G4double GetNumberOfInteractionLengthLeft () const
 
G4double GetTotalNumberOfInteractionLengthTraversed () const
 
G4bool isAtRestDoItIsEnabled () const
 
G4bool isAlongStepDoItIsEnabled () const
 
G4bool isPostStepDoItIsEnabled () const
 
virtual void DumpInfo () const
 
void SetVerboseLevel (G4int value)
 
G4int GetVerboseLevel () const
 
virtual void SetMasterProcess (G4VProcess *masterP)
 
const G4VProcessGetMasterProcess () const
 
virtual void BuildWorkerPhysicsTable (const G4ParticleDefinition &part)
 
virtual void PrepareWorkerPhysicsTable (const G4ParticleDefinition &)
 

Protected Attributes

G4PhysicsTablethePhysicsTable
 
- Protected Attributes inherited from G4VProcess
const G4ProcessManageraProcessManager
 
G4VParticleChange * pParticleChange
 
G4ParticleChange aParticleChange
 
G4double theNumberOfInteractionLengthLeft
 
G4double currentInteractionLength
 
G4double theInitialNumberOfInteractionLength
 
G4String theProcessName
 
G4String thePhysicsTableFileName
 
G4ProcessType theProcessType
 
G4int theProcessSubType
 
G4double thePILfactor
 
G4bool enableAtRestDoIt
 
G4bool enableAlongStepDoIt
 
G4bool enablePostStepDoIt
 
G4int verboseLevel
 

Private Member Functions

 G4OpRayleigh (const G4OpRayleigh &right)
 
G4OpRayleighoperator= (const G4OpRayleigh &right)
 
G4PhysicsOrderedFreeVectorCalculateRayleighMeanFreePaths (const G4Material *material) const
 

Additional Inherited Members

- Static Public Member Functions inherited from G4VProcess
static const G4StringGetProcessTypeName (G4ProcessType)
 
- Protected Member Functions inherited from G4VProcess
void SubtractNumberOfInteractionLengthLeft (G4double previousStepSize)
 
void ClearNumberOfInteractionLengthLeft ()
 

Detailed Description

Definition at line 77 of file G4OpRayleigh.hh.

Constructor & Destructor Documentation

◆ G4OpRayleigh() [1/2]

G4OpRayleigh::G4OpRayleigh ( const G4String processName = "OpRayleigh",
G4ProcessType  type = fOptical 
)

Definition at line 88 of file G4OpRayleigh.cc.

89  : G4VDiscreteProcess(processName, type)
90 {
92 
93  thePhysicsTable = NULL;
94 
95  if (verboseLevel>0) {
96  G4cout << GetProcessName() << " is created " << G4endl;
97  }
98 }
G4int verboseLevel
Definition: G4VProcess.hh:368
G4PhysicsTable * thePhysicsTable
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4GLOB_DLL std::ostream G4cout
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
#define G4endl
Definition: G4ios.hh:61
Here is the call graph for this function:

◆ ~G4OpRayleigh()

G4OpRayleigh::~G4OpRayleigh ( )

Definition at line 108 of file G4OpRayleigh.cc.

109 {
110  if (thePhysicsTable) {
112  delete thePhysicsTable;
113  }
114 }
G4PhysicsTable * thePhysicsTable
void clearAndDestroy()
Here is the call graph for this function:

◆ G4OpRayleigh() [2/2]

G4OpRayleigh::G4OpRayleigh ( const G4OpRayleigh right)
private

Member Function Documentation

◆ BuildPhysicsTable()

void G4OpRayleigh::BuildPhysicsTable ( const G4ParticleDefinition aParticleType)
virtual

Reimplemented from G4VProcess.

Definition at line 217 of file G4OpRayleigh.cc.

218 {
219  if (thePhysicsTable) {
221  delete thePhysicsTable;
222  thePhysicsTable = NULL;
223  }
224 
225  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
226  const G4int numOfMaterials = G4Material::GetNumberOfMaterials();
227 
228  thePhysicsTable = new G4PhysicsTable( numOfMaterials );
229 
230  for( G4int iMaterial = 0; iMaterial < numOfMaterials; iMaterial++ )
231  {
232  G4Material* material = (*theMaterialTable)[iMaterial];
233  G4MaterialPropertiesTable* materialProperties =
234  material->GetMaterialPropertiesTable();
235  G4PhysicsOrderedFreeVector* rayleigh = NULL;
236  if ( materialProperties != NULL ) {
237  rayleigh = materialProperties->GetProperty( "RAYLEIGH" );
238  if ( rayleigh == NULL ) rayleigh =
239  CalculateRayleighMeanFreePaths( material );
240  }
241  thePhysicsTable->insertAt( iMaterial, rayleigh );
242  }
243 }
G4MaterialPropertyVector * GetProperty(const char *key)
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
Definition: G4Material.hh:252
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:589
std::vector< G4Material * > G4MaterialTable
int G4int
Definition: G4Types.hh:78
G4PhysicsTable * thePhysicsTable
string material
Definition: eplot.py:19
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:596
void insertAt(size_t, G4PhysicsVector *)
G4PhysicsOrderedFreeVector * CalculateRayleighMeanFreePaths(const G4Material *material) const
void clearAndDestroy()
Here is the call graph for this function:

◆ CalculateRayleighMeanFreePaths()

G4PhysicsOrderedFreeVector * G4OpRayleigh::CalculateRayleighMeanFreePaths ( const G4Material material) const
private

Calculates the mean free paths for a material as a function of photon energy

Parameters
[in]materialinformation
Returns
the mean free path vector

Definition at line 269 of file G4OpRayleigh.cc.

270 {
271  G4MaterialPropertiesTable* materialProperties =
272  material->GetMaterialPropertiesTable();
273 
274  // Retrieve the beta_T or isothermal compressibility value. For backwards
275  // compatibility use a constant if the material is "Water". If the material
276  // doesn't have an ISOTHERMAL_COMPRESSIBILITY constant then return
277  G4double betat;
278  if ( material->GetName() == "Water" )
279  betat = 7.658e-23*m3/MeV;
280  else if(materialProperties->ConstPropertyExists("ISOTHERMAL_COMPRESSIBILITY"))
281  betat = materialProperties->GetConstProperty("ISOTHERMAL_COMPRESSIBILITY");
282  else
283  return NULL;
284 
285  // If the material doesn't have a RINDEX property vector then return
286  G4MaterialPropertyVector* rIndex = materialProperties->GetProperty("RINDEX");
287  if ( rIndex == NULL ) return NULL;
288 
289  // Retrieve the optional scale factor, (this just scales the scattering length
290  G4double scaleFactor = 1.0;
291  if( materialProperties->ConstPropertyExists( "RS_SCALE_FACTOR" ) )
292  scaleFactor= materialProperties->GetConstProperty("RS_SCALE_FACTOR" );
293 
294  // Retrieve the material temperature. For backwards compatibility use a
295  // constant if the material is "Water"
296  G4double temperature;
297  if( material->GetName() == "Water" )
298  temperature = 283.15*kelvin; // Temperature of water is 10 degrees celsius
299  else
300  temperature = material->GetTemperature();
301 
302  G4PhysicsOrderedFreeVector* rayleighMeanFreePaths =
304  // This calculates the meanFreePath via the Einstein-Smoluchowski formula
305  const G4double c1 = scaleFactor * betat * temperature * k_Boltzmann /
306  ( 6.0 * pi );
307 
308  for( size_t uRIndex = 0; uRIndex < rIndex->GetVectorLength(); uRIndex++ )
309  {
310  const G4double energy = rIndex->Energy( uRIndex );
311  const G4double rIndexSquared = (*rIndex)[uRIndex] * (*rIndex)[uRIndex];
312  const G4double xlambda = h_Planck * c_light / energy;
313  const G4double c2 = std::pow(twopi/xlambda,4);
314  const G4double c3 =
315  std::pow(((rIndexSquared-1.0)*(rIndexSquared+2.0 )/3.0),2);
316 
317  const G4double meanFreePath = 1.0 / ( c1 * c2 * c3 );
318 
319  if( verboseLevel>0 )
320  G4cout << energy << "MeV\t" << meanFreePath << "mm" << G4endl;
321 
322  rayleighMeanFreePaths->InsertValues( energy, meanFreePath );
323  }
324 
325  return rayleighMeanFreePaths;
326 }
static const double MeV
Definition: G4SIunits.hh:211
G4MaterialPropertyVector * GetProperty(const char *key)
G4int verboseLevel
Definition: G4VProcess.hh:368
float h_Planck
Definition: hepunit.py:263
static const double m3
Definition: G4SIunits.hh:130
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
Definition: G4Material.hh:252
void InsertValues(G4double energy, G4double value)
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
float k_Boltzmann
Definition: hepunit.py:299
static const G4double c3
static const double twopi
Definition: G4SIunits.hh:75
size_t GetVectorLength() const
G4double GetTemperature() const
Definition: G4Material.hh:182
static const double kelvin
Definition: G4SIunits.hh:278
TCanvas * c2
Definition: plot_hist.C:75
static const double pi
Definition: G4SIunits.hh:74
G4bool ConstPropertyExists(const char *key)
#define G4endl
Definition: G4ios.hh:61
G4double Energy(size_t index) const
double G4double
Definition: G4Types.hh:76
G4double GetConstProperty(const char *key)
const G4String & GetName() const
Definition: G4Material.hh:178
float c_light
Definition: hepunit.py:257
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DumpPhysicsTable()

void G4OpRayleigh::DumpPhysicsTable ( ) const
inline

Definition at line 167 of file G4OpRayleigh.hh.

169 {
170  G4int PhysicsTableSize = thePhysicsTable->entries();
172 
173  for (G4int i = 0 ; i < PhysicsTableSize ; i++ )
174  {
176  v->DumpValues();
177  }
178 }
int G4int
Definition: G4Types.hh:78
G4PhysicsTable * thePhysicsTable
size_t entries() const
Here is the call graph for this function:

◆ GetMeanFreePath()

G4double G4OpRayleigh::GetMeanFreePath ( const G4Track &  aTrack,
G4double  ,
G4ForceCondition *   
)
virtual

Implements G4VDiscreteProcess.

Definition at line 248 of file G4OpRayleigh.cc.

251 {
252  const G4DynamicParticle* particle = aTrack.GetDynamicParticle();
253  const G4double photonMomentum = particle->GetTotalMomentum();
254  const G4Material* material = aTrack.GetMaterial();
255 
256  G4PhysicsOrderedFreeVector* rayleigh =
257  static_cast<G4PhysicsOrderedFreeVector*>
258  ((*thePhysicsTable)(material->GetIndex()));
259 
260  G4double rsLength = DBL_MAX;
261  if( rayleigh != NULL ) rsLength = rayleigh->Value( photonMomentum );
262  return rsLength;
263 }
G4double GetTotalMomentum() const
size_t GetIndex() const
Definition: G4Material.hh:262
string material
Definition: eplot.py:19
G4double Value(G4double theEnergy, size_t &lastidx) const
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
Here is the call graph for this function:

◆ GetPhysicsTable()

G4PhysicsTable * G4OpRayleigh::GetPhysicsTable ( ) const
inline

Definition at line 180 of file G4OpRayleigh.hh.

181 {
182  return thePhysicsTable;
183 }
G4PhysicsTable * thePhysicsTable

◆ IsApplicable()

G4bool G4OpRayleigh::IsApplicable ( const G4ParticleDefinition aParticleType)
inlinevirtual

Reimplemented from G4VProcess.

Definition at line 161 of file G4OpRayleigh.hh.

162 {
163  return ( &aParticleType == G4OpticalPhoton::OpticalPhoton() );
164 }
static G4OpticalPhoton * OpticalPhoton()
Here is the call graph for this function:

◆ operator=()

G4OpRayleigh& G4OpRayleigh::operator= ( const G4OpRayleigh right)
private

◆ PostStepDoIt()

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

Reimplemented from G4VDiscreteProcess.

Definition at line 124 of file G4OpRayleigh.cc.

125 {
126  aParticleChange.Initialize(aTrack);
127 
128  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
129 
130  if (verboseLevel>0) {
131  G4cout << "Scattering Photon!" << G4endl;
132  G4cout << "Old Momentum Direction: "
133  << aParticle->GetMomentumDirection() << G4endl;
134  G4cout << "Old Polarization: "
135  << aParticle->GetPolarization() << G4endl;
136  }
137 
138  G4double cosTheta;
139  G4ThreeVector OldMomentumDirection, NewMomentumDirection;
140  G4ThreeVector OldPolarization, NewPolarization;
141 
142  G4double rand, constant;
143  G4double CosTheta, SinTheta, SinPhi, CosPhi, unit_x, unit_y, unit_z;
144 
145  do {
146  // Try to simulate the scattered photon momentum direction
147  // w.r.t. the initial photon momentum direction
148 
149  CosTheta = G4UniformRand();
150  SinTheta = std::sqrt(1.-CosTheta*CosTheta);
151  // consider for the angle 90-180 degrees
152  if (G4UniformRand() < 0.5) CosTheta = -CosTheta;
153 
154  // simulate the phi angle
155  rand = twopi*G4UniformRand();
156  SinPhi = std::sin(rand);
157  CosPhi = std::cos(rand);
158 
159  // start constructing the new momentum direction
160  unit_x = SinTheta * CosPhi;
161  unit_y = SinTheta * SinPhi;
162  unit_z = CosTheta;
163  NewMomentumDirection.set (unit_x,unit_y,unit_z);
164 
165  // Rotate the new momentum direction into global reference system
166  OldMomentumDirection = aParticle->GetMomentumDirection();
167  OldMomentumDirection = OldMomentumDirection.unit();
168  NewMomentumDirection.rotateUz(OldMomentumDirection);
169  NewMomentumDirection = NewMomentumDirection.unit();
170 
171  // calculate the new polarization direction
172  // The new polarization needs to be in the same plane as the new
173  // momentum direction and the old polarization direction
174  OldPolarization = aParticle->GetPolarization();
175  constant = -NewMomentumDirection.dot(OldPolarization);
176 
177  NewPolarization = OldPolarization + constant*NewMomentumDirection;
178  NewPolarization = NewPolarization.unit();
179 
180  // There is a corner case, where the Newmomentum direction
181  // is the same as oldpolariztion direction:
182  // random generate the azimuthal angle w.r.t. Newmomentum direction
183  if (NewPolarization.mag() == 0.) {
184  rand = G4UniformRand()*twopi;
185  NewPolarization.set(std::cos(rand),std::sin(rand),0.);
186  NewPolarization.rotateUz(NewMomentumDirection);
187  } else {
188  // There are two directions which are perpendicular
189  // to the new momentum direction
190  if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
191  }
192 
193  // simulate according to the distribution cos^2(theta)
194  cosTheta = NewPolarization.dot(OldPolarization);
195  // Loop checking, 13-Aug-2015, Peter Gumplinger
196  } while (std::pow(cosTheta,2) < G4UniformRand());
197 
198  aParticleChange.ProposePolarization(NewPolarization);
199  aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
200 
201  if (verboseLevel>0) {
202  G4cout << "New Polarization: "
203  << NewPolarization << G4endl;
204  G4cout << "Polarization Change: "
205  << *(aParticleChange.GetPolarization()) << G4endl;
206  G4cout << "New Momentum Direction: "
207  << NewMomentumDirection << G4endl;
208  G4cout << "Momentum Change: "
209  << *(aParticleChange.GetMomentumDirection()) << G4endl;
210  }
211 
212  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
213 }
void set(double x, double y, double z)
G4int verboseLevel
Definition: G4VProcess.hh:368
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double mag() const
Hep3Vector unit() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
static const double twopi
Definition: G4SIunits.hh:75
double dot(const Hep3Vector &) const
const G4ThreeVector & GetMomentumDirection() const
const G4ThreeVector & GetPolarization() const
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
Here is the call graph for this function:

Member Data Documentation

◆ thePhysicsTable

G4PhysicsTable* G4OpRayleigh::thePhysicsTable
protected

Definition at line 148 of file G4OpRayleigh.hh.


The documentation for this class was generated from the following files: