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