| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
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 // 7 //
>> 8 // $Id: G4OpRayleigh.cc,v 1.1.10.1 1999/12/07 20:52:52 gunter Exp $
>> 9 // GEANT4 tag $Name: geant4-01-00 $
28 // 10 //
>> 11 //
29 ////////////////////////////////////////////// 12 ////////////////////////////////////////////////////////////////////////
30 // Optical Photon Rayleigh Scattering Class Im 13 // Optical Photon Rayleigh Scattering Class Implementation
31 ////////////////////////////////////////////// 14 ////////////////////////////////////////////////////////////////////////
32 // 15 //
33 // File: G4OpRayleigh.cc << 16 // File: G4OpRayleigh.cc
34 // Description: Discrete Process -- Rayleigh s << 17 // Description: Discrete Process -- Rayleigh scattering of optical
35 // photons << 18 // photons
36 // Version: 1.0 19 // Version: 1.0
37 // Created: 1996-05-31 << 20 // Created: 1996-05-31
38 // Author: Juliet Armstrong 21 // Author: Juliet Armstrong
39 // Updated: 2014-10-10 - This version cal << 22 // Updated: 1997-04-09 by Peter Gumplinger
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 23 // > new physics/tracking scheme
>> 24 // mail: gum@triumf.ca
60 // 25 //
61 ////////////////////////////////////////////// 26 ////////////////////////////////////////////////////////////////////////
62 27
63 #include "G4OpRayleigh.hh" <<
64 #include "G4ios.hh" 28 #include "G4ios.hh"
65 #include "G4PhysicalConstants.hh" << 29 #include "G4OpRayleigh.hh"
66 #include "G4SystemOfUnits.hh" << 30
67 #include "G4OpticalParameters.hh" << 31 /////////////////////////
68 #include "G4OpProcessSubType.hh" << 32 // Class Implementation
69 << 33 /////////////////////////
70 //....oooOO0OOooo........oooOO0OOooo........oo << 34
71 G4OpRayleigh::G4OpRayleigh(const G4String& pro << 35 //////////////
72 : G4VDiscreteProcess(processName, type) << 36 // Operators
>> 37 //////////////
>> 38
>> 39 // G4OpRayleigh::operator=(const G4OpRayleigh &right)
>> 40 // {
>> 41 // }
>> 42
>> 43 /////////////////
>> 44 // Constructors
>> 45 /////////////////
>> 46
>> 47 G4OpRayleigh::G4OpRayleigh(const G4String& processName)
>> 48 : G4VDiscreteProcess(processName)
73 { 49 {
74 Initialise(); << 50
75 SetProcessSubType(fOpRayleigh); << 51 thePhysicsTable = NULL;
76 thePhysicsTable = nullptr; << 52
77 << 53 if (verboseLevel>0) {
78 if(verboseLevel > 0) << 54 G4cout << GetProcessName() << " is created " << endl;
79 { << 55 }
80 G4cout << GetProcessName() << " is created << 56
81 } << 57 BuildThePhysicsTable();
82 } 58 }
83 59
84 //....oooOO0OOooo........oooOO0OOooo........oo << 60 // G4OpRayleigh::G4OpRayleigh(const G4OpRayleigh &right)
>> 61 // {
>> 62 // }
>> 63
>> 64 ////////////////
>> 65 // Destructors
>> 66 ////////////////
>> 67
85 G4OpRayleigh::~G4OpRayleigh() 68 G4OpRayleigh::~G4OpRayleigh()
86 { 69 {
87 // VI: inside this PhysicsTable all properti << 70 if (thePhysicsTable!= NULL) {
88 // it is not possible to destroy << 71 thePhysicsTable->clearAndDestroy();
89 delete thePhysicsTable; << 72 delete thePhysicsTable;
>> 73 }
90 } 74 }
91 75
92 //....oooOO0OOooo........oooOO0OOooo........oo << 76 ////////////
93 void G4OpRayleigh::PreparePhysicsTable(const G << 77 // Methods
94 { << 78 ////////////
95 Initialise(); <<
96 } <<
97 79
98 //....oooOO0OOooo........oooOO0OOooo........oo << 80 // PostStepDoIt
99 void G4OpRayleigh::Initialise() << 81 // -------------
>> 82 //
>> 83 G4VParticleChange*
>> 84 G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
100 { 85 {
101 SetVerboseLevel(G4OpticalParameters::Instanc << 86 aParticleChange.Initialize(aTrack);
102 } <<
103 87
104 //....oooOO0OOooo........oooOO0OOooo........oo << 88 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
105 G4VParticleChange* G4OpRayleigh::PostStepDoIt( << 89 const G4Material* aMaterial = aTrack.GetMaterial();
106 <<
107 { <<
108 aParticleChange.Initialize(aTrack); <<
109 const G4DynamicParticle* aParticle = aTrack. <<
110 90
111 if(verboseLevel > 1) << 91 if (verboseLevel>0) {
112 { << 92 G4cout << "Scattering Photon!" << endl;
113 G4cout << "OpRayleigh: Scattering Photon!" << 93 G4cout << "Old Momentum Direction: "
114 << "Old Momentum Direction: " << aP << 94 << aParticle->GetMomentumDirection() << endl;
115 << G4endl << "Old Polarization: " < << 95 G4cout << "Old Polarization: "
116 << G4endl; << 96 << aParticle->GetPolarization() << endl;
117 } << 97 }
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 98
184 return G4VDiscreteProcess::PostStepDoIt(aTra << 99 // find polar angle w.r.t. old polarization vector
185 } <<
186 100
187 //....oooOO0OOooo........oooOO0OOooo........oo << 101 G4double rand = G4UniformRand();
188 void G4OpRayleigh::BuildPhysicsTable(const G4P << 102
189 { << 103 G4double CosTheta = pow(rand, 1./3.);
190 if(thePhysicsTable) << 104 G4double SinTheta = 1.-CosTheta*CosTheta;
191 { << 105
192 // thePhysicsTable->clearAndDestroy(); << 106 // find azimuthal angle w.r.t old polarization vector
193 delete thePhysicsTable; << 107
194 thePhysicsTable = nullptr; << 108 rand = G4UniformRand();
195 } << 109
196 << 110 G4double Phi = twopi*rand;
197 const G4MaterialTable* theMaterialTable = G4 << 111 G4double SinPhi = sin(Phi);
198 const size_t numOfMaterials = G4 << 112 G4double CosPhi = cos(Phi);
199 thePhysicsTable = ne << 113
200 << 114 G4double unit_x = SinTheta * CosPhi;
201 for(size_t i = 0; i < numOfMaterials; ++i) << 115 G4double unit_y = SinTheta * SinPhi;
202 { << 116 G4double unit_z = CosTheta;
203 G4Material* material = (*the << 117
204 G4MaterialPropertiesTable* matProp = mater << 118 G4ThreeVector NewPolarization (unit_x,unit_y,unit_z);
205 G4PhysicsFreeVector* rayleigh = nullptr; << 119
206 if(matProp) << 120 // Rotate new polarization direction into global reference system
207 { << 121
208 rayleigh = matProp->GetProperty(kRAYLEIG << 122 G4ThreeVector OldPolarization = aParticle->GetPolarization();
209 if(rayleigh == nullptr) << 123
210 rayleigh = CalculateRayleighMeanFreePa << 124 NewPolarization.rotateUz(OldPolarization);
211 } << 125
212 thePhysicsTable->insertAt(i, rayleigh); << 126 // -- new momentum direction is normal to the new polarization
213 } << 127 // vector (components below expressed in reference system where
>> 128 // new polarization vector is aligned with the z axis)
>> 129
>> 130 SinTheta = 1.0;
>> 131 CosTheta = 0.0;
>> 132
>> 133 rand = G4UniformRand();
>> 134
>> 135 Phi = twopi * rand;
>> 136 SinPhi = sin(Phi);
>> 137 CosPhi = cos(Phi);
>> 138
>> 139 unit_x = SinTheta*CosPhi;
>> 140 unit_y = SinTheta*SinPhi;
>> 141 unit_z = CosTheta;
>> 142
>> 143 G4ThreeVector NewMomentumDirection(unit_x,unit_y,unit_z);
>> 144
>> 145 // Rotate New momentum direction vector into global reference system
>> 146
>> 147 NewMomentumDirection.rotateUz(NewPolarization);
>> 148
>> 149 aParticleChange.SetPolarizationChange(NewPolarization.unit());
>> 150
>> 151 aParticleChange.SetMomentumChange(NewMomentumDirection.unit());
>> 152
>> 153 if (verboseLevel>0) {
>> 154 G4cout << "New Polarization: "
>> 155 << NewPolarization << endl;
>> 156 G4cout << "Polarization Change: "
>> 157 << *(aParticleChange.GetPolarizationChange()) << endl;
>> 158 G4cout << "New Momentum Direction: "
>> 159 << NewMomentumDirection << endl;
>> 160 G4cout << "Momentum Change: "
>> 161 << *(aParticleChange.GetMomentumChange()) << endl;
>> 162 }
>> 163
>> 164 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
214 } 165 }
215 166
216 //....oooOO0OOooo........oooOO0OOooo........oo << 167 // BuildThePhysicsTable for the Rayleigh Scattering process
217 G4double G4OpRayleigh::GetMeanFreePath(const G << 168 // --------------------------------------------------------
218 G4Force << 169 //
>> 170 void G4OpRayleigh::BuildThePhysicsTable()
219 { 171 {
220 auto rayleigh = static_cast<G4PhysicsFreeVec << 172 // Builds a table of scattering lengths for each material
221 (*thePhysicsTable)(aTrack.GetMaterial()- << 173
>> 174 if (thePhysicsTable) return;
>> 175
>> 176 const G4MaterialTable* theMaterialTable=
>> 177 G4Material::GetMaterialTable();
>> 178 G4int numOfMaterials = theMaterialTable->length();
>> 179
>> 180 // create a new physics table
222 181
223 G4double rsLength = DBL_MAX; << 182 thePhysicsTable = new G4PhysicsTable(numOfMaterials);
224 if(rayleigh) << 183
225 { << 184 // loop for materials
226 rsLength = rayleigh->Value(aTrack.GetDynam << 185
227 idx_rslength); << 186 for (G4int i=0 ; i < numOfMaterials; i++)
228 } << 187 {
229 return rsLength; << 188 G4PhysicsOrderedFreeVector* ScatteringLengths =
>> 189 new G4PhysicsOrderedFreeVector();
>> 190
>> 191 if ((*theMaterialTable)[i]->GetName() == "Water")
>> 192 {
>> 193 G4MaterialPropertiesTable *MaterialPT =
>> 194 (*theMaterialTable)[i]->GetMaterialPropertiesTable();
>> 195 // Call utility routine to Generate
>> 196 // Rayleigh Scattering Lengths
>> 197 ScatteringLengths =
>> 198 RayleighAttenuationLengthGenerator(MaterialPT);
>> 199 }
>> 200
>> 201 thePhysicsTable->insertAt(i,ScatteringLengths);
>> 202 }
230 } 203 }
231 204
232 //....oooOO0OOooo........oooOO0OOooo........oo << 205 // GetMeanFreePath()
233 G4PhysicsFreeVector* G4OpRayleigh::CalculateRa << 206 // -----------------
234 const G4Material* material) const << 207 //
>> 208 G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
>> 209 G4double ,
>> 210 G4ForceCondition* )
235 { 211 {
236 G4MaterialPropertiesTable* MPT = material->G << 212 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
>> 213 const G4Material* aMaterial = aTrack.GetMaterial();
237 214
238 // Retrieve the beta_T or isothermal compres << 215 G4double thePhotonMomentum = aParticle->GetTotalMomentum();
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 216
301 rayleighMFPs->InsertValues(energy, meanFre << 217 G4double AttenuationLength = DBL_MAX;
302 } << 218 G4bool isOutRange;
303 219
304 return rayleighMFPs; << 220 if (aMaterial->GetName() == "Water")
>> 221 {
>> 222 AttenuationLength =
>> 223 (*thePhysicsTable)(aMaterial->GetIndex())->
>> 224 GetValue(thePhotonMomentum, isOutRange);
>> 225 }
>> 226 return AttenuationLength;
305 } 227 }
306 228
307 //....oooOO0OOooo........oooOO0OOooo........oo << 229 // RayleighAttenuationLengthGenerator()
308 void G4OpRayleigh::SetVerboseLevel(G4int verbo << 230 // ------------------------------------
>> 231 // Private method to compute Rayleigh Scattering Lengths (for water)
>> 232 //
>> 233 G4PhysicsOrderedFreeVector*
>> 234 G4OpRayleigh::RayleighAttenuationLengthGenerator(G4MaterialPropertiesTable *aMPT)
309 { 235 {
310 verboseLevel = verbose; << 236 // Physical Constants
311 G4OpticalParameters::Instance()->SetRayleigh << 237
>> 238 // isothermal compressibility of water
>> 239 G4double betat = 7.658e-23*m3/MeV;
>> 240
>> 241 // K Boltzman
>> 242 G4double kboltz = 8.61739e-11*MeV/kelvin;
>> 243
>> 244 // Temperature of water is 10 degrees celsius
>> 245 // conversion to kelvin:
>> 246 // TCelsius = TKelvin - 273.15 => 273.15 + 10 = 283.15
>> 247 G4double temp = 283.15*kelvin;
>> 248
>> 249 // Retrieve vectors for refraction index
>> 250 // and photon momentum from the material properties table
>> 251
>> 252 G4MaterialPropertyVector* Rindex = aMPT->GetProperty("RINDEX");
>> 253
>> 254 G4double refsq;
>> 255 G4double e;
>> 256 G4double xlambda;
>> 257 G4double c1, c2, c3, c4;
>> 258 G4double Dist;
>> 259 G4double refraction_index;
>> 260
>> 261 G4double no_unit = 1.0;
>> 262
>> 263 G4PhysicsOrderedFreeVector *RayleighScatteringLengths =
>> 264 new G4PhysicsOrderedFreeVector();
>> 265 Rindex->ResetIterator();
>> 266
>> 267 while (++(*Rindex)) {
>> 268
>> 269 e = (Rindex->GetPhotonMomentum());
>> 270
>> 271 refraction_index = Rindex->GetProperty();
>> 272 refsq = refraction_index*refraction_index;
>> 273 xlambda = h_Planck*c_light/e;
>> 274
>> 275 if (verboseLevel>0) {
>> 276 G4cout << Rindex->GetPhotonMomentum() << " MeV\t";
>> 277 G4cout << xlambda << " mm\t";
>> 278 }
>> 279
>> 280 c1 = 1 / (6.0 * pi);
>> 281 c2 = pow((2.0 * pi / xlambda), 4);
>> 282 c3 = pow( ( (refsq - 1.0) * (refsq + 2.0) / 3.0 ), 2);
>> 283 c4 = betat * temp * kboltz;
>> 284
>> 285 Dist = 1.0 / (c1*c2*c3*c4);
>> 286
>> 287 if (verboseLevel>0) {
>> 288 G4cout << Dist << " mm" << endl;
>> 289 }
>> 290 RayleighScatteringLengths->
>> 291 InsertValues(Rindex->GetPhotonMomentum(), Dist);
>> 292 }
>> 293
>> 294 return RayleighScatteringLengths;
312 } 295 }
313 296