| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // 27 //
28 // 28 //
29 ////////////////////////////////////////////// 29 ////////////////////////////////////////////////////////////////////////
30 // Optical Photon Rayleigh Scattering Class Im 30 // Optical Photon Rayleigh Scattering Class Implementation
31 ////////////////////////////////////////////// 31 ////////////////////////////////////////////////////////////////////////
32 // 32 //
33 // File: G4OpRayleigh.cc 33 // File: G4OpRayleigh.cc
34 // Description: Discrete Process -- Rayleigh s 34 // Description: Discrete Process -- Rayleigh scattering of optical
35 // photons 35 // photons
36 // Version: 1.0 36 // Version: 1.0
37 // Created: 1996-05-31 37 // Created: 1996-05-31
38 // Author: Juliet Armstrong 38 // Author: Juliet Armstrong
39 // Updated: 2014-10-10 - This version cal 39 // Updated: 2014-10-10 - This version calculates the Rayleigh scattering
40 // length for more materials than 40 // length for more materials than just Water (although the Water
41 // default is kept). To do this t 41 // default is kept). To do this the user would need to specify the
42 // ISOTHERMAL_COMPRESSIBILITY as 42 // ISOTHERMAL_COMPRESSIBILITY as a material property and
43 // optionally an RS_SCALE_LENGTH 43 // optionally an RS_SCALE_LENGTH (useful for testing). Code comes
44 // from Philip Graham (Queen Mary 44 // from Philip Graham (Queen Mary University of London).
45 // 2010-06-11 - Fix Bug 207; Than 45 // 2010-06-11 - Fix Bug 207; Thanks to Xin Qian
46 // (Kellogg Radiation Lab of Calt 46 // (Kellogg Radiation Lab of Caltech)
47 // 2005-07-28 - add G4ProcessType 47 // 2005-07-28 - add G4ProcessType to constructor
48 // 2001-10-18 by Peter Gumplinger 48 // 2001-10-18 by Peter Gumplinger
49 // eliminate unused variable warn 49 // eliminate unused variable warning on Linux (gcc-2.95.2)
50 // 2001-09-18 by mma 50 // 2001-09-18 by mma
51 // >numOfMaterials=G4Material::Ge 51 // >numOfMaterials=G4Material::GetNumberOfMaterials() in BuildPhy
52 // 2001-01-30 by Peter Gumplinger 52 // 2001-01-30 by Peter Gumplinger
53 // > allow for positiv and negati 53 // > allow for positiv and negative CosTheta and force the
54 // > new momentum direction to be 54 // > new momentum direction to be in the same plane as the
55 // > new and old polarization vec 55 // > new and old polarization vectors
56 // 2001-01-29 by Peter Gumplinger 56 // 2001-01-29 by Peter Gumplinger
57 // > fix calculation of SinTheta 57 // > fix calculation of SinTheta (from CosTheta)
58 // 1997-04-09 by Peter Gumplinger 58 // 1997-04-09 by Peter Gumplinger
59 // > new physics/tracking scheme 59 // > new physics/tracking scheme
60 // 60 //
61 ////////////////////////////////////////////// 61 ////////////////////////////////////////////////////////////////////////
62 62
63 #include "G4OpRayleigh.hh" 63 #include "G4OpRayleigh.hh"
64 #include "G4ios.hh" 64 #include "G4ios.hh"
65 #include "G4PhysicalConstants.hh" 65 #include "G4PhysicalConstants.hh"
66 #include "G4SystemOfUnits.hh" 66 #include "G4SystemOfUnits.hh"
67 #include "G4OpticalParameters.hh" 67 #include "G4OpticalParameters.hh"
68 #include "G4OpProcessSubType.hh" 68 #include "G4OpProcessSubType.hh"
69 69
70 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 G4OpRayleigh::G4OpRayleigh(const G4String& pro 71 G4OpRayleigh::G4OpRayleigh(const G4String& processName, G4ProcessType type)
72 : G4VDiscreteProcess(processName, type) 72 : G4VDiscreteProcess(processName, type)
73 { 73 {
74 Initialise(); 74 Initialise();
75 SetProcessSubType(fOpRayleigh); 75 SetProcessSubType(fOpRayleigh);
76 thePhysicsTable = nullptr; 76 thePhysicsTable = nullptr;
77 77
78 if(verboseLevel > 0) 78 if(verboseLevel > 0)
79 { 79 {
80 G4cout << GetProcessName() << " is created 80 G4cout << GetProcessName() << " is created " << G4endl;
81 } 81 }
82 } 82 }
83 83
84 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 G4OpRayleigh::~G4OpRayleigh() 85 G4OpRayleigh::~G4OpRayleigh()
86 { 86 {
87 // VI: inside this PhysicsTable all properti 87 // VI: inside this PhysicsTable all properties are unique
88 // it is not possible to destroy 88 // it is not possible to destroy
89 delete thePhysicsTable; 89 delete thePhysicsTable;
90 } 90 }
91 91
92 //....oooOO0OOooo........oooOO0OOooo........oo 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93 void G4OpRayleigh::PreparePhysicsTable(const G 93 void G4OpRayleigh::PreparePhysicsTable(const G4ParticleDefinition&)
94 { 94 {
95 Initialise(); 95 Initialise();
96 } 96 }
97 97
98 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 void G4OpRayleigh::Initialise() 99 void G4OpRayleigh::Initialise()
100 { 100 {
101 SetVerboseLevel(G4OpticalParameters::Instanc 101 SetVerboseLevel(G4OpticalParameters::Instance()->GetRayleighVerboseLevel());
102 } 102 }
103 103
104 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 G4VParticleChange* G4OpRayleigh::PostStepDoIt( 105 G4VParticleChange* G4OpRayleigh::PostStepDoIt(const G4Track& aTrack,
106 106 const G4Step& aStep)
107 { 107 {
108 aParticleChange.Initialize(aTrack); 108 aParticleChange.Initialize(aTrack);
109 const G4DynamicParticle* aParticle = aTrack. 109 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
110 110
111 if(verboseLevel > 1) 111 if(verboseLevel > 1)
112 { 112 {
113 G4cout << "OpRayleigh: Scattering Photon!" 113 G4cout << "OpRayleigh: Scattering Photon!" << G4endl
114 << "Old Momentum Direction: " << aP 114 << "Old Momentum Direction: " << aParticle->GetMomentumDirection()
115 << G4endl << "Old Polarization: " < 115 << G4endl << "Old Polarization: " << aParticle->GetPolarization()
116 << G4endl; 116 << G4endl;
117 } 117 }
118 118
119 G4double cosTheta; 119 G4double cosTheta;
120 G4ThreeVector oldMomDir, newMomDir; 120 G4ThreeVector oldMomDir, newMomDir;
121 G4ThreeVector oldPol, newPol; 121 G4ThreeVector oldPol, newPol;
122 G4double rand; 122 G4double rand;
123 G4double cost, sint, sinphi, cosphi; 123 G4double cost, sint, sinphi, cosphi;
124 124
125 do 125 do
126 { 126 {
127 // Try to simulate the scattered photon mo 127 // Try to simulate the scattered photon momentum direction
128 // w.r.t. the initial photon momentum dire 128 // w.r.t. the initial photon momentum direction
129 cost = G4UniformRand(); 129 cost = G4UniformRand();
130 sint = std::sqrt(1. - cost * cost); 130 sint = std::sqrt(1. - cost * cost);
131 // consider for the angle 90-180 degrees 131 // consider for the angle 90-180 degrees
132 if(G4UniformRand() < 0.5) 132 if(G4UniformRand() < 0.5)
133 cost = -cost; 133 cost = -cost;
134 134
135 // simulate the phi angle 135 // simulate the phi angle
136 rand = twopi * G4UniformRand(); 136 rand = twopi * G4UniformRand();
137 sinphi = std::sin(rand); 137 sinphi = std::sin(rand);
138 cosphi = std::cos(rand); 138 cosphi = std::cos(rand);
139 139
140 // construct the new momentum direction 140 // construct the new momentum direction
141 newMomDir.set(sint * cosphi, sint * sinphi 141 newMomDir.set(sint * cosphi, sint * sinphi, cost);
142 oldMomDir = aParticle->GetMomentumDirectio 142 oldMomDir = aParticle->GetMomentumDirection();
143 newMomDir.rotateUz(oldMomDir); 143 newMomDir.rotateUz(oldMomDir);
144 144
145 // calculate the new polarization directio 145 // calculate the new polarization direction
146 // The new polarization needs to be in the 146 // The new polarization needs to be in the same plane as the new
147 // momentum direction and the old polariza 147 // momentum direction and the old polarization direction
148 oldPol = aParticle->GetPolarization(); 148 oldPol = aParticle->GetPolarization();
149 newPol = (oldPol - newMomDir.dot(oldPol) * 149 newPol = (oldPol - newMomDir.dot(oldPol) * newMomDir).unit();
150 150
151 // There is a corner case, where the new m 151 // There is a corner case, where the new momentum direction
152 // is the same as old polarization directi 152 // is the same as old polarization direction:
153 // random generate the azimuthal angle w.r 153 // random generate the azimuthal angle w.r.t. new momentum direction
154 if(newPol.mag() == 0.) 154 if(newPol.mag() == 0.)
155 { 155 {
156 rand = G4UniformRand() * twopi; 156 rand = G4UniformRand() * twopi;
157 newPol.set(std::cos(rand), std::sin(rand 157 newPol.set(std::cos(rand), std::sin(rand), 0.);
158 newPol.rotateUz(newMomDir); 158 newPol.rotateUz(newMomDir);
159 } 159 }
160 else 160 else
161 { 161 {
162 // There are two directions perpendicula 162 // There are two directions perpendicular to the new momentum direction
163 if(G4UniformRand() < 0.5) 163 if(G4UniformRand() < 0.5)
164 newPol = -newPol; 164 newPol = -newPol;
165 } 165 }
166 166
167 // simulate according to the distribution 167 // simulate according to the distribution cos^2(theta)
168 cosTheta = newPol.dot(oldPol); 168 cosTheta = newPol.dot(oldPol);
169 // Loop checking, 13-Aug-2015, Peter Gumpl 169 // Loop checking, 13-Aug-2015, Peter Gumplinger
170 } while(std::pow(cosTheta, 2) < G4UniformRan 170 } while(std::pow(cosTheta, 2) < G4UniformRand());
171 171
172 aParticleChange.ProposePolarization(newPol); 172 aParticleChange.ProposePolarization(newPol);
173 aParticleChange.ProposeMomentumDirection(new 173 aParticleChange.ProposeMomentumDirection(newMomDir);
174 174
175 if(verboseLevel > 1) 175 if(verboseLevel > 1)
176 { 176 {
177 G4cout << "New Polarization: " << newPol < 177 G4cout << "New Polarization: " << newPol << G4endl
178 << "Polarization Change: " << *(aPa 178 << "Polarization Change: " << *(aParticleChange.GetPolarization())
179 << G4endl << "New Momentum Directio 179 << G4endl << "New Momentum Direction: " << newMomDir << G4endl
180 << "Momentum Change: " << *(aPartic 180 << "Momentum Change: " << *(aParticleChange.GetMomentumDirection())
181 << G4endl; 181 << G4endl;
182 } 182 }
183 183
184 return G4VDiscreteProcess::PostStepDoIt(aTra 184 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
185 } 185 }
186 186
187 //....oooOO0OOooo........oooOO0OOooo........oo 187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
188 void G4OpRayleigh::BuildPhysicsTable(const G4P 188 void G4OpRayleigh::BuildPhysicsTable(const G4ParticleDefinition&)
189 { 189 {
190 if(thePhysicsTable) 190 if(thePhysicsTable)
191 { 191 {
192 // thePhysicsTable->clearAndDestroy(); 192 // thePhysicsTable->clearAndDestroy();
193 delete thePhysicsTable; 193 delete thePhysicsTable;
194 thePhysicsTable = nullptr; 194 thePhysicsTable = nullptr;
195 } 195 }
196 196
197 const G4MaterialTable* theMaterialTable = G4 197 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
198 const size_t numOfMaterials = G4 198 const size_t numOfMaterials = G4Material::GetNumberOfMaterials();
199 thePhysicsTable = ne 199 thePhysicsTable = new G4PhysicsTable(numOfMaterials);
200 200
201 for(size_t i = 0; i < numOfMaterials; ++i) 201 for(size_t i = 0; i < numOfMaterials; ++i)
202 { 202 {
203 G4Material* material = (*the 203 G4Material* material = (*theMaterialTable)[i];
204 G4MaterialPropertiesTable* matProp = mater 204 G4MaterialPropertiesTable* matProp = material->GetMaterialPropertiesTable();
205 G4PhysicsFreeVector* rayleigh = nullptr; 205 G4PhysicsFreeVector* rayleigh = nullptr;
206 if(matProp) 206 if(matProp)
207 { 207 {
208 rayleigh = matProp->GetProperty(kRAYLEIG 208 rayleigh = matProp->GetProperty(kRAYLEIGH);
209 if(rayleigh == nullptr) 209 if(rayleigh == nullptr)
210 rayleigh = CalculateRayleighMeanFreePa 210 rayleigh = CalculateRayleighMeanFreePaths(material);
211 } 211 }
212 thePhysicsTable->insertAt(i, rayleigh); 212 thePhysicsTable->insertAt(i, rayleigh);
213 } 213 }
214 } 214 }
215 215
216 //....oooOO0OOooo........oooOO0OOooo........oo 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217 G4double G4OpRayleigh::GetMeanFreePath(const G 217 G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack, G4double,
218 G4Force 218 G4ForceCondition*)
219 { 219 {
220 auto rayleigh = static_cast<G4PhysicsFreeVec 220 auto rayleigh = static_cast<G4PhysicsFreeVector*>(
221 (*thePhysicsTable)(aTrack.GetMaterial()- 221 (*thePhysicsTable)(aTrack.GetMaterial()->GetIndex()));
222 222
223 G4double rsLength = DBL_MAX; 223 G4double rsLength = DBL_MAX;
224 if(rayleigh) 224 if(rayleigh)
225 { 225 {
226 rsLength = rayleigh->Value(aTrack.GetDynam 226 rsLength = rayleigh->Value(aTrack.GetDynamicParticle()->GetTotalMomentum(),
227 idx_rslength); 227 idx_rslength);
228 } 228 }
229 return rsLength; 229 return rsLength;
230 } 230 }
231 231
232 //....oooOO0OOooo........oooOO0OOooo........oo 232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
233 G4PhysicsFreeVector* G4OpRayleigh::CalculateRa 233 G4PhysicsFreeVector* G4OpRayleigh::CalculateRayleighMeanFreePaths(
234 const G4Material* material) const 234 const G4Material* material) const
235 { 235 {
236 G4MaterialPropertiesTable* MPT = material->G 236 G4MaterialPropertiesTable* MPT = material->GetMaterialPropertiesTable();
237 237
238 // Retrieve the beta_T or isothermal compres 238 // Retrieve the beta_T or isothermal compressibility value. For backwards
239 // compatibility use a constant if the mater 239 // compatibility use a constant if the material is "Water". If the material
240 // doesn't have an ISOTHERMAL_COMPRESSIBILIT 240 // doesn't have an ISOTHERMAL_COMPRESSIBILITY constant then return
241 G4double betat; 241 G4double betat;
242 if(material->GetName() == "Water") 242 if(material->GetName() == "Water")
243 { 243 {
244 betat = 7.658e-23 * m3 / MeV; 244 betat = 7.658e-23 * m3 / MeV;
245 } 245 }
246 else if(MPT->ConstPropertyExists(kISOTHERMAL 246 else if(MPT->ConstPropertyExists(kISOTHERMAL_COMPRESSIBILITY))
247 { 247 {
248 betat = MPT->GetConstProperty(kISOTHERMAL_ 248 betat = MPT->GetConstProperty(kISOTHERMAL_COMPRESSIBILITY);
249 } 249 }
250 else 250 else
251 { 251 {
252 return nullptr; 252 return nullptr;
253 } 253 }
254 254
255 // If the material doesn't have a RINDEX pro 255 // If the material doesn't have a RINDEX property vector then return
256 G4MaterialPropertyVector* rIndex = MPT->GetP 256 G4MaterialPropertyVector* rIndex = MPT->GetProperty(kRINDEX);
257 if(rIndex == nullptr) 257 if(rIndex == nullptr)
258 return nullptr; 258 return nullptr;
259 259
260 // Retrieve the optional scale factor (scale 260 // Retrieve the optional scale factor (scales the scattering length)
261 G4double scaleFactor = 1.0; 261 G4double scaleFactor = 1.0;
262 if(MPT->ConstPropertyExists(kRS_SCALE_FACTOR 262 if(MPT->ConstPropertyExists(kRS_SCALE_FACTOR))
263 { 263 {
264 scaleFactor = MPT->GetConstProperty(kRS_SC 264 scaleFactor = MPT->GetConstProperty(kRS_SCALE_FACTOR);
265 } 265 }
266 266
267 // Retrieve the material temperature. For ba 267 // Retrieve the material temperature. For backwards compatibility use a
268 // constant if the material is "Water" 268 // constant if the material is "Water"
269 G4double temperature; 269 G4double temperature;
270 if(material->GetName() == "Water") 270 if(material->GetName() == "Water")
271 { 271 {
272 temperature = 272 temperature =
273 283.15 * kelvin; // Temperature of wate 273 283.15 * kelvin; // Temperature of water is 10 degrees celsius
274 } 274 }
275 else 275 else
276 { 276 {
277 temperature = material->GetTemperature(); 277 temperature = material->GetTemperature();
278 } 278 }
279 279
280 auto rayleighMFPs = new G4PhysicsFreeVector( 280 auto rayleighMFPs = new G4PhysicsFreeVector();
281 // This calculates the meanFreePath via the 281 // This calculates the meanFreePath via the Einstein-Smoluchowski formula
282 const G4double c1 = 282 const G4double c1 =
283 scaleFactor * betat * temperature * k_Bolt 283 scaleFactor * betat * temperature * k_Boltzmann / (6.0 * pi);
284 284
285 for(size_t uRIndex = 0; uRIndex < rIndex->Ge 285 for(size_t uRIndex = 0; uRIndex < rIndex->GetVectorLength(); ++uRIndex)
286 { 286 {
287 const G4double energy = rIndex->Ene 287 const G4double energy = rIndex->Energy(uRIndex);
288 const G4double rIndexSquared = (*rIndex)[u 288 const G4double rIndexSquared = (*rIndex)[uRIndex] * (*rIndex)[uRIndex];
289 const G4double xlambda = h_Planck * 289 const G4double xlambda = h_Planck * c_light / energy;
290 const G4double c2 = std::pow(tw 290 const G4double c2 = std::pow(twopi / xlambda, 4);
291 const G4double c3 = 291 const G4double c3 =
292 std::pow(((rIndexSquared - 1.0) * (rInde 292 std::pow(((rIndexSquared - 1.0) * (rIndexSquared + 2.0) / 3.0), 2);
293 293
294 const G4double meanFreePath = 1.0 / (c1 * 294 const G4double meanFreePath = 1.0 / (c1 * c2 * c3);
295 295
296 if(verboseLevel > 0) 296 if(verboseLevel > 0)
297 { 297 {
298 G4cout << energy << "MeV\t" << meanFreeP 298 G4cout << energy << "MeV\t" << meanFreePath << "mm" << G4endl;
299 } 299 }
300 300
301 rayleighMFPs->InsertValues(energy, meanFre 301 rayleighMFPs->InsertValues(energy, meanFreePath);
302 } 302 }
303 303
304 return rayleighMFPs; 304 return rayleighMFPs;
305 } 305 }
306 306
307 //....oooOO0OOooo........oooOO0OOooo........oo 307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308 void G4OpRayleigh::SetVerboseLevel(G4int verbo 308 void G4OpRayleigh::SetVerboseLevel(G4int verbose)
309 { 309 {
310 verboseLevel = verbose; 310 verboseLevel = verbose;
311 G4OpticalParameters::Instance()->SetRayleigh 311 G4OpticalParameters::Instance()->SetRayleighVerboseLevel(verboseLevel);
312 } 312 }
313 313