| 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 // $Id: G4OpRayleigh.cc,v 1.19 2010-10-29 23:18:35 gum Exp $
>> 28 // GEANT4 tag $Name: not supported by cvs2svn $
27 // 29 //
28 // << 30 //
29 ////////////////////////////////////////////// 31 ////////////////////////////////////////////////////////////////////////
30 // Optical Photon Rayleigh Scattering Class Im 32 // Optical Photon Rayleigh Scattering Class Implementation
31 ////////////////////////////////////////////// 33 ////////////////////////////////////////////////////////////////////////
32 // 34 //
33 // File: G4OpRayleigh.cc << 35 // File: G4OpRayleigh.cc
34 // Description: Discrete Process -- Rayleigh s << 36 // Description: Discrete Process -- Rayleigh scattering of optical
35 // photons << 37 // photons
36 // Version: 1.0 38 // Version: 1.0
37 // Created: 1996-05-31 << 39 // Created: 1996-05-31
38 // Author: Juliet Armstrong 40 // Author: Juliet Armstrong
39 // Updated: 2014-10-10 - This version cal << 41 // Updated: 2010-06-11 - Fix Bug 207; Thanks to Xin Qian
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 42 // (Kellogg Radiation Lab of Caltech)
47 // 2005-07-28 - add G4ProcessType 43 // 2005-07-28 - add G4ProcessType to constructor
48 // 2001-10-18 by Peter Gumplinger 44 // 2001-10-18 by Peter Gumplinger
49 // eliminate unused variable warn 45 // eliminate unused variable warning on Linux (gcc-2.95.2)
50 // 2001-09-18 by mma 46 // 2001-09-18 by mma
51 // >numOfMaterials=G4Material::Ge << 47 // >numOfMaterials=G4Material::GetNumberOfMaterials() in BuildPhy
52 // 2001-01-30 by Peter Gumplinger 48 // 2001-01-30 by Peter Gumplinger
53 // > allow for positiv and negati 49 // > allow for positiv and negative CosTheta and force the
54 // > new momentum direction to be 50 // > new momentum direction to be in the same plane as the
55 // > new and old polarization vec 51 // > new and old polarization vectors
56 // 2001-01-29 by Peter Gumplinger 52 // 2001-01-29 by Peter Gumplinger
57 // > fix calculation of SinTheta 53 // > fix calculation of SinTheta (from CosTheta)
58 // 1997-04-09 by Peter Gumplinger 54 // 1997-04-09 by Peter Gumplinger
59 // > new physics/tracking scheme 55 // > new physics/tracking scheme
>> 56 // mail: gum@triumf.ca
60 // 57 //
61 ////////////////////////////////////////////// 58 ////////////////////////////////////////////////////////////////////////
62 59
63 #include "G4OpRayleigh.hh" 60 #include "G4OpRayleigh.hh"
>> 61
64 #include "G4ios.hh" 62 #include "G4ios.hh"
65 #include "G4PhysicalConstants.hh" 63 #include "G4PhysicalConstants.hh"
66 #include "G4SystemOfUnits.hh" 64 #include "G4SystemOfUnits.hh"
67 #include "G4OpticalParameters.hh" <<
68 #include "G4OpProcessSubType.hh" 65 #include "G4OpProcessSubType.hh"
69 66
70 //....oooOO0OOooo........oooOO0OOooo........oo << 67 /////////////////////////
>> 68 // Class Implementation
>> 69 /////////////////////////
>> 70
>> 71 //////////////
>> 72 // Operators
>> 73 //////////////
>> 74
>> 75 // G4OpRayleigh::operator=(const G4OpRayleigh &right)
>> 76 // {
>> 77 // }
>> 78
>> 79 /////////////////
>> 80 // Constructors
>> 81 /////////////////
>> 82
71 G4OpRayleigh::G4OpRayleigh(const G4String& pro 83 G4OpRayleigh::G4OpRayleigh(const G4String& processName, G4ProcessType type)
72 : G4VDiscreteProcess(processName, type) << 84 : G4VDiscreteProcess(processName, type)
73 { 85 {
74 Initialise(); << 86 SetProcessSubType(fOpRayleigh);
75 SetProcessSubType(fOpRayleigh); << 87
76 thePhysicsTable = nullptr; << 88 thePhysicsTable = 0;
77 << 89
78 if(verboseLevel > 0) << 90 DefaultWater = false;
79 { << 91
80 G4cout << GetProcessName() << " is created << 92 if (verboseLevel>0) {
81 } << 93 G4cout << GetProcessName() << " is created " << G4endl;
>> 94 }
>> 95
>> 96 BuildThePhysicsTable();
82 } 97 }
83 98
84 //....oooOO0OOooo........oooOO0OOooo........oo << 99 // G4OpRayleigh::G4OpRayleigh(const G4OpRayleigh &right)
>> 100 // {
>> 101 // }
>> 102
>> 103 ////////////////
>> 104 // Destructors
>> 105 ////////////////
>> 106
85 G4OpRayleigh::~G4OpRayleigh() 107 G4OpRayleigh::~G4OpRayleigh()
86 { 108 {
87 // VI: inside this PhysicsTable all properti << 109 if (thePhysicsTable!= 0) {
88 // it is not possible to destroy << 110 thePhysicsTable->clearAndDestroy();
89 delete thePhysicsTable; << 111 delete thePhysicsTable;
>> 112 }
90 } 113 }
91 114
92 //....oooOO0OOooo........oooOO0OOooo........oo << 115 ////////////
93 void G4OpRayleigh::PreparePhysicsTable(const G << 116 // Methods
94 { << 117 ////////////
95 Initialise(); <<
96 } <<
97 118
98 //....oooOO0OOooo........oooOO0OOooo........oo << 119 // PostStepDoIt
99 void G4OpRayleigh::Initialise() << 120 // -------------
>> 121 //
>> 122 G4VParticleChange*
>> 123 G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
100 { 124 {
101 SetVerboseLevel(G4OpticalParameters::Instanc << 125 aParticleChange.Initialize(aTrack);
102 } <<
103 126
104 //....oooOO0OOooo........oooOO0OOooo........oo << 127 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
105 G4VParticleChange* G4OpRayleigh::PostStepDoIt( <<
106 <<
107 { <<
108 aParticleChange.Initialize(aTrack); <<
109 const G4DynamicParticle* aParticle = aTrack. <<
110 128
111 if(verboseLevel > 1) << 129 if (verboseLevel>0) {
112 { << 130 G4cout << "Scattering Photon!" << G4endl;
113 G4cout << "OpRayleigh: Scattering Photon!" << 131 G4cout << "Old Momentum Direction: "
114 << "Old Momentum Direction: " << aP << 132 << aParticle->GetMomentumDirection() << G4endl;
115 << G4endl << "Old Polarization: " < << 133 G4cout << "Old Polarization: "
116 << G4endl; << 134 << aParticle->GetPolarization() << G4endl;
117 } << 135 }
118 << 136
119 G4double cosTheta; << 137 G4double cosTheta;
120 G4ThreeVector oldMomDir, newMomDir; << 138 G4ThreeVector OldMomentumDirection, NewMomentumDirection;
121 G4ThreeVector oldPol, newPol; << 139 G4ThreeVector OldPolarization, NewPolarization;
122 G4double rand; << 140
123 G4double cost, sint, sinphi, cosphi; << 141 do {
124 << 142 // Try to simulate the scattered photon momentum direction
125 do << 143 // w.r.t. the initial photon momentum direction
126 { << 144
127 // Try to simulate the scattered photon mo << 145 G4double CosTheta = G4UniformRand();
128 // w.r.t. the initial photon momentum dire << 146 G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta);
129 cost = G4UniformRand(); << 147 // consider for the angle 90-180 degrees
130 sint = std::sqrt(1. - cost * cost); << 148 if (G4UniformRand() < 0.5) CosTheta = -CosTheta;
131 // consider for the angle 90-180 degrees << 149
132 if(G4UniformRand() < 0.5) << 150 // simulate the phi angle
133 cost = -cost; << 151 G4double rand = twopi*G4UniformRand();
134 << 152 G4double SinPhi = std::sin(rand);
135 // simulate the phi angle << 153 G4double CosPhi = std::cos(rand);
136 rand = twopi * G4UniformRand(); << 154
137 sinphi = std::sin(rand); << 155 // start constructing the new momentum direction
138 cosphi = std::cos(rand); << 156 G4double unit_x = SinTheta * CosPhi;
139 << 157 G4double unit_y = SinTheta * SinPhi;
140 // construct the new momentum direction << 158 G4double unit_z = CosTheta;
141 newMomDir.set(sint * cosphi, sint * sinphi << 159 NewMomentumDirection.set (unit_x,unit_y,unit_z);
142 oldMomDir = aParticle->GetMomentumDirectio << 160
143 newMomDir.rotateUz(oldMomDir); << 161 // Rotate the new momentum direction into global reference system
144 << 162 OldMomentumDirection = aParticle->GetMomentumDirection();
145 // calculate the new polarization directio << 163 OldMomentumDirection = OldMomentumDirection.unit();
146 // The new polarization needs to be in the << 164 NewMomentumDirection.rotateUz(OldMomentumDirection);
147 // momentum direction and the old polariza << 165 NewMomentumDirection = NewMomentumDirection.unit();
148 oldPol = aParticle->GetPolarization(); << 166
149 newPol = (oldPol - newMomDir.dot(oldPol) * << 167 // calculate the new polarization direction
150 << 168 // The new polarization needs to be in the same plane as the new
151 // There is a corner case, where the new m << 169 // momentum direction and the old polarization direction
152 // is the same as old polarization directi << 170 OldPolarization = aParticle->GetPolarization();
153 // random generate the azimuthal angle w.r << 171 G4double constant = -1./NewMomentumDirection.dot(OldPolarization);
154 if(newPol.mag() == 0.) << 172
155 { << 173 NewPolarization = NewMomentumDirection + constant*OldPolarization;
156 rand = G4UniformRand() * twopi; << 174 NewPolarization = NewPolarization.unit();
157 newPol.set(std::cos(rand), std::sin(rand << 175
158 newPol.rotateUz(newMomDir); << 176 // There is a corner case, where the Newmomentum direction
159 } << 177 // is the same as oldpolariztion direction:
160 else << 178 // random generate the azimuthal angle w.r.t. Newmomentum direction
161 { << 179 if (NewPolarization.mag() == 0.) {
162 // There are two directions perpendicula << 180 rand = G4UniformRand()*twopi;
163 if(G4UniformRand() < 0.5) << 181 NewPolarization.set(std::cos(rand),std::sin(rand),0.);
164 newPol = -newPol; << 182 NewPolarization.rotateUz(NewMomentumDirection);
165 } << 183 } else {
166 << 184 // There are two directions which are perpendicular
167 // simulate according to the distribution << 185 // to the new momentum direction
168 cosTheta = newPol.dot(oldPol); << 186 if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
169 // Loop checking, 13-Aug-2015, Peter Gumpl << 187 }
170 } while(std::pow(cosTheta, 2) < G4UniformRan << 188
171 << 189 // simulate according to the distribution cos^2(theta)
172 aParticleChange.ProposePolarization(newPol); << 190 cosTheta = NewPolarization.dot(OldPolarization);
173 aParticleChange.ProposeMomentumDirection(new << 191 } while (std::pow(cosTheta,2) < G4UniformRand());
174 << 192
175 if(verboseLevel > 1) << 193 aParticleChange.ProposePolarization(NewPolarization);
176 { << 194 aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
177 G4cout << "New Polarization: " << newPol < << 195
178 << "Polarization Change: " << *(aPa << 196 if (verboseLevel>0) {
179 << G4endl << "New Momentum Directio << 197 G4cout << "New Polarization: "
180 << "Momentum Change: " << *(aPartic << 198 << NewPolarization << G4endl;
181 << G4endl; << 199 G4cout << "Polarization Change: "
182 } << 200 << *(aParticleChange.GetPolarization()) << G4endl;
>> 201 G4cout << "New Momentum Direction: "
>> 202 << NewMomentumDirection << G4endl;
>> 203 G4cout << "Momentum Change: "
>> 204 << *(aParticleChange.GetMomentumDirection()) << G4endl;
>> 205 }
183 206
184 return G4VDiscreteProcess::PostStepDoIt(aTra << 207 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
185 } 208 }
186 209
187 //....oooOO0OOooo........oooOO0OOooo........oo << 210 // BuildThePhysicsTable for the Rayleigh Scattering process
188 void G4OpRayleigh::BuildPhysicsTable(const G4P << 211 // --------------------------------------------------------
>> 212 //
>> 213 void G4OpRayleigh::BuildThePhysicsTable()
189 { 214 {
190 if(thePhysicsTable) << 215 // Builds a table of scattering lengths for each material
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
216 //....oooOO0OOooo........oooOO0OOooo........oo << 217 if (thePhysicsTable) return;
217 G4double G4OpRayleigh::GetMeanFreePath(const G << 218
218 G4Force << 219 const G4MaterialTable* theMaterialTable=
219 { << 220 G4Material::GetMaterialTable();
220 auto rayleigh = static_cast<G4PhysicsFreeVec << 221 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
221 (*thePhysicsTable)(aTrack.GetMaterial()- << 222
>> 223 // create a new physics table
>> 224
>> 225 thePhysicsTable = new G4PhysicsTable(numOfMaterials);
>> 226
>> 227 // loop for materials
>> 228
>> 229 for (G4int i=0 ; i < numOfMaterials; i++)
>> 230 {
>> 231 G4PhysicsOrderedFreeVector* ScatteringLengths = NULL;
>> 232
>> 233 G4MaterialPropertiesTable *aMaterialPropertiesTable =
>> 234 (*theMaterialTable)[i]->GetMaterialPropertiesTable();
>> 235
>> 236 if(aMaterialPropertiesTable){
>> 237
>> 238 G4MaterialPropertyVector* AttenuationLengthVector =
>> 239 aMaterialPropertiesTable->GetProperty("RAYLEIGH");
>> 240
>> 241 if(!AttenuationLengthVector){
>> 242
>> 243 if ((*theMaterialTable)[i]->GetName() == "Water")
>> 244 {
>> 245 // Call utility routine to Generate
>> 246 // Rayleigh Scattering Lengths
222 247
223 G4double rsLength = DBL_MAX; << 248 DefaultWater = true;
224 if(rayleigh) << 249
225 { << 250 ScatteringLengths =
226 rsLength = rayleigh->Value(aTrack.GetDynam << 251 RayleighAttenuationLengthGenerator(aMaterialPropertiesTable);
227 idx_rslength); << 252 }
228 } << 253 }
229 return rsLength; << 254 }
>> 255
>> 256 thePhysicsTable->insertAt(i,ScatteringLengths);
>> 257 }
230 } 258 }
231 259
232 //....oooOO0OOooo........oooOO0OOooo........oo << 260 // GetMeanFreePath()
233 G4PhysicsFreeVector* G4OpRayleigh::CalculateRa << 261 // -----------------
234 const G4Material* material) const << 262 //
>> 263 G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
>> 264 G4double ,
>> 265 G4ForceCondition* )
235 { 266 {
236 G4MaterialPropertiesTable* MPT = material->G << 267 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
>> 268 const G4Material* aMaterial = aTrack.GetMaterial();
237 269
238 // Retrieve the beta_T or isothermal compres << 270 G4double thePhotonEnergy = aParticle->GetTotalEnergy();
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 271
301 rayleighMFPs->InsertValues(energy, meanFre << 272 G4double AttenuationLength = DBL_MAX;
302 } <<
303 273
304 return rayleighMFPs; << 274 if (aMaterial->GetName() == "Water" && DefaultWater){
>> 275
>> 276 G4bool isOutRange;
>> 277
>> 278 AttenuationLength =
>> 279 (*thePhysicsTable)(aMaterial->GetIndex())->
>> 280 GetValue(thePhotonEnergy, isOutRange);
>> 281 }
>> 282 else {
>> 283
>> 284 G4MaterialPropertiesTable* aMaterialPropertyTable =
>> 285 aMaterial->GetMaterialPropertiesTable();
>> 286
>> 287 if(aMaterialPropertyTable){
>> 288 G4MaterialPropertyVector* AttenuationLengthVector =
>> 289 aMaterialPropertyTable->GetProperty("RAYLEIGH");
>> 290 if(AttenuationLengthVector){
>> 291 AttenuationLength = AttenuationLengthVector ->
>> 292 Value(thePhotonEnergy);
>> 293 }
>> 294 else{
>> 295 // G4cout << "No Rayleigh scattering length specified" << G4endl;
>> 296 }
>> 297 }
>> 298 else{
>> 299 // G4cout << "No Rayleigh scattering length specified" << G4endl;
>> 300 }
>> 301 }
>> 302
>> 303 return AttenuationLength;
305 } 304 }
306 305
307 //....oooOO0OOooo........oooOO0OOooo........oo << 306 // RayleighAttenuationLengthGenerator()
308 void G4OpRayleigh::SetVerboseLevel(G4int verbo << 307 // ------------------------------------
>> 308 // Private method to compute Rayleigh Scattering Lengths (for water)
>> 309 //
>> 310 G4PhysicsOrderedFreeVector*
>> 311 G4OpRayleigh::RayleighAttenuationLengthGenerator(G4MaterialPropertiesTable *aMPT)
309 { 312 {
310 verboseLevel = verbose; << 313 // Physical Constants
311 G4OpticalParameters::Instance()->SetRayleigh << 314
>> 315 // isothermal compressibility of water
>> 316 G4double betat = 7.658e-23*m3/MeV;
>> 317
>> 318 // K Boltzman
>> 319 G4double kboltz = 8.61739e-11*MeV/kelvin;
>> 320
>> 321 // Temperature of water is 10 degrees celsius
>> 322 // conversion to kelvin:
>> 323 // TCelsius = TKelvin - 273.15 => 273.15 + 10 = 283.15
>> 324 G4double temp = 283.15*kelvin;
>> 325
>> 326 // Retrieve vectors for refraction index
>> 327 // and photon energy from the material properties table
>> 328
>> 329 G4MaterialPropertyVector* Rindex = aMPT->GetProperty("RINDEX");
>> 330
>> 331 G4double refsq;
>> 332 G4double e;
>> 333 G4double xlambda;
>> 334 G4double c1, c2, c3, c4;
>> 335 G4double Dist;
>> 336 G4double refraction_index;
>> 337
>> 338 G4PhysicsOrderedFreeVector *RayleighScatteringLengths =
>> 339 new G4PhysicsOrderedFreeVector();
>> 340
>> 341 if (Rindex ) {
>> 342
>> 343 for (size_t i = 0; i < Rindex->GetVectorLength(); i++) {
>> 344
>> 345 e = Rindex->Energy(i);
>> 346
>> 347 refraction_index = (*Rindex)[i];
>> 348
>> 349 refsq = refraction_index*refraction_index;
>> 350 xlambda = h_Planck*c_light/e;
>> 351
>> 352 if (verboseLevel>0) {
>> 353 G4cout << Rindex->Energy(i) << " MeV\t";
>> 354 G4cout << xlambda << " mm\t";
>> 355 }
>> 356
>> 357 c1 = 1 / (6.0 * pi);
>> 358 c2 = std::pow((2.0 * pi / xlambda), 4);
>> 359 c3 = std::pow( ( (refsq - 1.0) * (refsq + 2.0) / 3.0 ), 2);
>> 360 c4 = betat * temp * kboltz;
>> 361
>> 362 Dist = 1.0 / (c1*c2*c3*c4);
>> 363
>> 364 if (verboseLevel>0) {
>> 365 G4cout << Dist << " mm" << G4endl;
>> 366 }
>> 367 RayleighScatteringLengths->
>> 368 InsertValues(Rindex->Energy(i), Dist);
>> 369 }
>> 370
>> 371 }
>> 372
>> 373 return RayleighScatteringLengths;
312 } 374 }
313 375