Geant4  10.01.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 84717 2014-10-20 07:39:47Z 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  do {
143  // Try to simulate the scattered photon momentum direction
144  // w.r.t. the initial photon momentum direction
145 
146  G4double CosTheta = G4UniformRand();
147  G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta);
148  // consider for the angle 90-180 degrees
149  if (G4UniformRand() < 0.5) CosTheta = -CosTheta;
150 
151  // simulate the phi angle
152  G4double rand = twopi*G4UniformRand();
153  G4double SinPhi = std::sin(rand);
154  G4double CosPhi = std::cos(rand);
155 
156  // start constructing the new momentum direction
157  G4double unit_x = SinTheta * CosPhi;
158  G4double unit_y = SinTheta * SinPhi;
159  G4double unit_z = CosTheta;
160  NewMomentumDirection.set (unit_x,unit_y,unit_z);
161 
162  // Rotate the new momentum direction into global reference system
163  OldMomentumDirection = aParticle->GetMomentumDirection();
164  OldMomentumDirection = OldMomentumDirection.unit();
165  NewMomentumDirection.rotateUz(OldMomentumDirection);
166  NewMomentumDirection = NewMomentumDirection.unit();
167 
168  // calculate the new polarization direction
169  // The new polarization needs to be in the same plane as the new
170  // momentum direction and the old polarization direction
171  OldPolarization = aParticle->GetPolarization();
172  G4double constant = -1./NewMomentumDirection.dot(OldPolarization);
173 
174  NewPolarization = NewMomentumDirection + constant*OldPolarization;
175  NewPolarization = NewPolarization.unit();
176 
177  // There is a corner case, where the Newmomentum direction
178  // is the same as oldpolariztion direction:
179  // random generate the azimuthal angle w.r.t. Newmomentum direction
180  if (NewPolarization.mag() == 0.) {
181  rand = G4UniformRand()*twopi;
182  NewPolarization.set(std::cos(rand),std::sin(rand),0.);
183  NewPolarization.rotateUz(NewMomentumDirection);
184  } else {
185  // There are two directions which are perpendicular
186  // to the new momentum direction
187  if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
188  }
189 
190  // simulate according to the distribution cos^2(theta)
191  cosTheta = NewPolarization.dot(OldPolarization);
192  } while (std::pow(cosTheta,2) < G4UniformRand());
193 
194  aParticleChange.ProposePolarization(NewPolarization);
195  aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
196 
197  if (verboseLevel>0) {
198  G4cout << "New Polarization: "
199  << NewPolarization << G4endl;
200  G4cout << "Polarization Change: "
201  << *(aParticleChange.GetPolarization()) << G4endl;
202  G4cout << "New Momentum Direction: "
203  << NewMomentumDirection << G4endl;
204  G4cout << "Momentum Change: "
205  << *(aParticleChange.GetMomentumDirection()) << G4endl;
206  }
207 
208  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
209 }
210 
211 // BuildPhysicsTable for the Rayleigh Scattering process
212 // --------------------------------------------------------
214 {
215  if (thePhysicsTable) {
217  delete thePhysicsTable;
218  thePhysicsTable = NULL;
219  }
220 
221  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
222  const G4int numOfMaterials = G4Material::GetNumberOfMaterials();
223 
224  thePhysicsTable = new G4PhysicsTable( numOfMaterials );
225 
226  for( G4int iMaterial = 0; iMaterial < numOfMaterials; iMaterial++ )
227  {
228  G4Material* material = (*theMaterialTable)[iMaterial];
229  G4MaterialPropertiesTable* materialProperties =
230  material->GetMaterialPropertiesTable();
231  G4PhysicsOrderedFreeVector* rayleigh = NULL;
232  if ( materialProperties != NULL ) {
233  rayleigh = materialProperties->GetProperty( "RAYLEIGH" );
234  if ( rayleigh == NULL ) rayleigh =
235  CalculateRayleighMeanFreePaths( material );
236  }
237  thePhysicsTable->insertAt( iMaterial, rayleigh );
238  }
239 }
240 
241 // GetMeanFreePath()
242 // -----------------
243 //
245  G4double ,
247 {
248  const G4DynamicParticle* particle = aTrack.GetDynamicParticle();
249  const G4double photonMomentum = particle->GetTotalMomentum();
250  const G4Material* material = aTrack.GetMaterial();
251 
252  G4PhysicsOrderedFreeVector* rayleigh =
253  static_cast<G4PhysicsOrderedFreeVector*>
254  ((*thePhysicsTable)(material->GetIndex()));
255 
256  G4double rsLength = DBL_MAX;
257  if( rayleigh != NULL ) rsLength = rayleigh->Value( photonMomentum );
258  return rsLength;
259 }
260 
261 // CalculateRayleighMeanFreePaths()
262 // --------------------------------
263 // Private method to compute Rayleigh Scattering Lengths
266 {
267  G4MaterialPropertiesTable* materialProperties =
268  material->GetMaterialPropertiesTable();
269 
270  // Retrieve the beta_T or isothermal compressibility value. For backwards
271  // compatibility use a constant if the material is "Water". If the material
272  // doesn't have an ISOTHERMAL_COMPRESSIBILITY constant then return
273  G4double betat;
274  if ( material->GetName() == "Water" )
275  betat = 7.658e-23*m3/MeV;
276  else if(materialProperties->ConstPropertyExists("ISOTHERMAL_COMPRESSIBILITY"))
277  betat = materialProperties->GetConstProperty("ISOTHERMAL_COMPRESSIBILITY");
278  else
279  return NULL;
280 
281  // If the material doesn't have a RINDEX property vector then return
282  G4MaterialPropertyVector* rIndex = materialProperties->GetProperty("RINDEX");
283  if ( rIndex == NULL ) return NULL;
284 
285  // Retrieve the optional scale factor, (this just scales the scattering length
286  G4double scaleFactor = 1.0;
287  if( materialProperties->ConstPropertyExists( "RS_SCALE_FACTOR" ) )
288  scaleFactor= materialProperties->GetConstProperty("RS_SCALE_FACTOR" );
289 
290  // Retrieve the material temperature. For backwards compatibility use a
291  // constant if the material is "Water"
292  G4double temperature;
293  if( material->GetName() == "Water" )
294  temperature = 283.15*kelvin; // Temperature of water is 10 degrees celsius
295  else
296  temperature = material->GetTemperature();
297 
298  G4PhysicsOrderedFreeVector* rayleighMeanFreePaths =
300  // This calculates the meanFreePath via the Einstein-Smoluchowski formula
301  const G4double c1 = scaleFactor * betat * temperature * k_Boltzmann /
302  ( 6.0 * pi );
303 
304  for( size_t uRIndex = 0; uRIndex < rIndex->GetVectorLength(); uRIndex++ )
305  {
306  const G4double energy = rIndex->Energy( uRIndex );
307  const G4double rIndexSquared = (*rIndex)[uRIndex] * (*rIndex)[uRIndex];
308  const G4double xlambda = h_Planck * c_light / energy;
309  const G4double c2 = std::pow(twopi/xlambda,4);
310  const G4double c3 =
311  std::pow(((rIndexSquared-1.0)*(rIndexSquared+2.0 )/3.0),2);
312 
313  const G4double meanFreePath = 1.0 / ( c1 * c2 * c3 );
314 
315  if( verboseLevel>0 )
316  G4cout << energy << "MeV\t" << meanFreePath << "mm" << G4endl;
317 
318  rayleighMeanFreePaths->InsertValues( energy, meanFreePath );
319  }
320 
321  return rayleighMeanFreePaths;
322 }
static c2_factory< G4double > c2
static const double MeV
Definition: G4SIunits.hh:193
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:112
size_t GetIndex() const
Definition: G4Material.hh:262
const G4String & GetName() const
Definition: G4Material.hh:178
const G4double pi
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:588
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:93
G4GLOB_DLL std::ostream G4cout
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
const G4ThreeVector & GetMomentumDirection() const
static const G4double c3
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:595
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4Material * GetMaterial() const
static const double kelvin
Definition: G4SIunits.hh:260
virtual void Initialize(const G4Track &)
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