Geant4  10.02.p02
G4OpRayleigh.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4OpRayleigh.cc 92045 2015-08-14 07:21:23Z gcosmo $
28 //
29 //
31 // Optical Photon Rayleigh Scattering Class Implementation
33 //
34 // File: G4OpRayleigh.cc
35 // Description: Discrete Process -- Rayleigh scattering of optical
36 // photons
37 // Version: 1.0
38 // Created: 1996-05-31
39 // Author: Juliet Armstrong
40 // Updated: 2014-10-10 - This version calculates the Rayleigh scattering
41 // length for more materials than just Water (although the Water
42 // default is kept). To do this the user would need to specify the
43 // ISOTHERMAL_COMPRESSIBILITY as a material property and
44 // optionally an RS_SCALE_LENGTH (useful for testing). Code comes
45 // from Philip Graham (Queen Mary University of London).
46 // 2010-06-11 - Fix Bug 207; Thanks to Xin Qian
47 // (Kellogg Radiation Lab of Caltech)
48 // 2005-07-28 - add G4ProcessType to constructor
49 // 2001-10-18 by Peter Gumplinger
50 // eliminate unused variable warning on Linux (gcc-2.95.2)
51 // 2001-09-18 by mma
52 // >numOfMaterials=G4Material::GetNumberOfMaterials() in BuildPhy
53 // 2001-01-30 by Peter Gumplinger
54 // > allow for positiv and negative CosTheta and force the
55 // > new momentum direction to be in the same plane as the
56 // > new and old polarization vectors
57 // 2001-01-29 by Peter Gumplinger
58 // > fix calculation of SinTheta (from CosTheta)
59 // 1997-04-09 by Peter Gumplinger
60 // > new physics/tracking scheme
61 // mail: gum@triumf.ca
62 //
64 
65 #include "G4OpRayleigh.hh"
66 
67 #include "G4ios.hh"
68 #include "G4PhysicalConstants.hh"
69 #include "G4SystemOfUnits.hh"
70 #include "G4OpProcessSubType.hh"
71 
73 // Class Implementation
75 
77  // Operators
79 
80 // G4OpRayleigh::operator=(const G4OpRayleigh &right)
81 // {
82 // }
83 
85  // Constructors
87 
89  : G4VDiscreteProcess(processName, type)
90 {
92 
93  thePhysicsTable = NULL;
94 
95  if (verboseLevel>0) {
96  G4cout << GetProcessName() << " is created " << G4endl;
97  }
98 }
99 
100 // G4OpRayleigh::G4OpRayleigh(const G4OpRayleigh &right)
101 // {
102 // }
103 
105  // Destructors
107 
109 {
110  if (thePhysicsTable) {
112  delete thePhysicsTable;
113  }
114 }
115 
117  // Methods
119 
120 // PostStepDoIt
121 // -------------
122 //
124 G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
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 }
214 
215 // BuildPhysicsTable for the Rayleigh Scattering process
216 // --------------------------------------------------------
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 }
244 
245 // GetMeanFreePath()
246 // -----------------
247 //
249  G4double ,
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 }
264 
265 // CalculateRayleighMeanFreePaths()
266 // --------------------------------
267 // Private method to compute Rayleigh Scattering Lengths
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 c2_factory< G4double > c2
static const double MeV
Definition: G4SIunits.hh:211
G4MaterialPropertyVector * GetProperty(const char *key)
G4int verboseLevel
Definition: G4VProcess.hh:368
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
static const double m3
Definition: G4SIunits.hh:130
size_t GetIndex() const
Definition: G4Material.hh:262
const G4String & GetName() const
Definition: G4Material.hh:178
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:589
std::vector< G4Material * > G4MaterialTable
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
void InsertValues(G4double energy, G4double value)
const G4ThreeVector * GetMomentumDirection() const
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
G4PhysicsTable * thePhysicsTable
G4double GetTotalMomentum() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
const G4ThreeVector & GetMomentumDirection() const
static const G4double c3
static const double twopi
Definition: G4SIunits.hh:75
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
Definition: G4Step.hh:76
static const G4double c1
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:596
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4Material * GetMaterial() const
static const double kelvin
Definition: G4SIunits.hh:278
virtual void Initialize(const G4Track &)
static const double pi
Definition: G4SIunits.hh:74
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
Definition: G4Material.hh:252
G4double energy(const ThreeVector &p, const G4double m)
G4OpRayleigh(const G4String &processName="OpRayleigh", G4ProcessType type=fOptical)
Definition: G4OpRayleigh.cc:88
const G4ThreeVector & GetPolarization() const
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
G4bool ConstPropertyExists(const char *key)
void insertAt(size_t, G4PhysicsVector *)
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
G4PhysicsOrderedFreeVector * CalculateRayleighMeanFreePaths(const G4Material *material) const
Calculates the mean free paths for a material as a function of photon energy.
G4ForceCondition
const G4ThreeVector * GetPolarization() const
G4double GetConstProperty(const char *key)
#define DBL_MAX
Definition: templates.hh:83
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void clearAndDestroy()
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
G4ProcessType