source: trunk/source/processes/optical/src/G4OpRayleigh.cc@ 899

Last change on this file since 899 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 11.4 KB
Line 
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,v 1.14 2006/06/29 21:08:54 gunter Exp $
28// GEANT4 tag $Name: $
29//
30//
31////////////////////////////////////////////////////////////////////////
32// Optical Photon Rayleigh Scattering Class Implementation
33////////////////////////////////////////////////////////////////////////
34//
35// File: G4OpRayleigh.cc
36// Description: Discrete Process -- Rayleigh scattering of optical
37// photons
38// Version: 1.0
39// Created: 1996-05-31
40// Author: Juliet Armstrong
41// Updated: 2005-07-28 - add G4ProcessType to constructor
42// 2001-10-18 by Peter Gumplinger
43// eliminate unused variable warning on Linux (gcc-2.95.2)
44// 2001-09-18 by mma
45// >numOfMaterials=G4Material::GetNumberOfMaterials() in BuildPhy
46// 2001-01-30 by Peter Gumplinger
47// > allow for positiv and negative CosTheta and force the
48// > new momentum direction to be in the same plane as the
49// > new and old polarization vectors
50// 2001-01-29 by Peter Gumplinger
51// > fix calculation of SinTheta (from CosTheta)
52// 1997-04-09 by Peter Gumplinger
53// > new physics/tracking scheme
54// mail: gum@triumf.ca
55//
56////////////////////////////////////////////////////////////////////////
57
58#include "G4ios.hh"
59#include "G4OpRayleigh.hh"
60
61/////////////////////////
62// Class Implementation
63/////////////////////////
64
65 //////////////
66 // Operators
67 //////////////
68
69// G4OpRayleigh::operator=(const G4OpRayleigh &right)
70// {
71// }
72
73 /////////////////
74 // Constructors
75 /////////////////
76
77G4OpRayleigh::G4OpRayleigh(const G4String& processName, G4ProcessType type)
78 : G4VDiscreteProcess(processName, type)
79{
80
81 thePhysicsTable = 0;
82
83 DefaultWater = false;
84
85 if (verboseLevel>0) {
86 G4cout << GetProcessName() << " is created " << G4endl;
87 }
88
89 BuildThePhysicsTable();
90}
91
92// G4OpRayleigh::G4OpRayleigh(const G4OpRayleigh &right)
93// {
94// }
95
96 ////////////////
97 // Destructors
98 ////////////////
99
100G4OpRayleigh::~G4OpRayleigh()
101{
102 if (thePhysicsTable!= 0) {
103 thePhysicsTable->clearAndDestroy();
104 delete thePhysicsTable;
105 }
106}
107
108 ////////////
109 // Methods
110 ////////////
111
112// PostStepDoIt
113// -------------
114//
115G4VParticleChange*
116G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
117{
118 aParticleChange.Initialize(aTrack);
119
120 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
121
122 if (verboseLevel>0) {
123 G4cout << "Scattering Photon!" << G4endl;
124 G4cout << "Old Momentum Direction: "
125 << aParticle->GetMomentumDirection() << G4endl;
126 G4cout << "Old Polarization: "
127 << aParticle->GetPolarization() << G4endl;
128 }
129
130 // find polar angle w.r.t. old polarization vector
131
132 G4double rand = G4UniformRand();
133
134 G4double CosTheta = std::pow(rand, 1./3.);
135 G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta);
136
137 if(G4UniformRand() < 0.5)CosTheta = -CosTheta;
138
139 // find azimuthal angle w.r.t old polarization vector
140
141 rand = G4UniformRand();
142
143 G4double Phi = twopi*rand;
144 G4double SinPhi = std::sin(Phi);
145 G4double CosPhi = std::cos(Phi);
146
147 G4double unit_x = SinTheta * CosPhi;
148 G4double unit_y = SinTheta * SinPhi;
149 G4double unit_z = CosTheta;
150
151 G4ThreeVector NewPolarization (unit_x,unit_y,unit_z);
152
153 // Rotate new polarization direction into global reference system
154
155 G4ThreeVector OldPolarization = aParticle->GetPolarization();
156 OldPolarization = OldPolarization.unit();
157
158 NewPolarization.rotateUz(OldPolarization);
159 NewPolarization = NewPolarization.unit();
160
161 // -- new momentum direction is normal to the new
162 // polarization vector and in the same plane as the
163 // old and new polarization vectors --
164
165 G4ThreeVector NewMomentumDirection =
166 OldPolarization - NewPolarization * CosTheta;
167
168 if(G4UniformRand() < 0.5)NewMomentumDirection = -NewMomentumDirection;
169 NewMomentumDirection = NewMomentumDirection.unit();
170
171 aParticleChange.ProposePolarization(NewPolarization);
172
173 aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
174
175 if (verboseLevel>0) {
176 G4cout << "New Polarization: "
177 << NewPolarization << G4endl;
178 G4cout << "Polarization Change: "
179 << *(aParticleChange.GetPolarization()) << G4endl;
180 G4cout << "New Momentum Direction: "
181 << NewMomentumDirection << G4endl;
182 G4cout << "Momentum Change: "
183 << *(aParticleChange.GetMomentumDirection()) << G4endl;
184 }
185
186 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
187}
188
189// BuildThePhysicsTable for the Rayleigh Scattering process
190// --------------------------------------------------------
191//
192void G4OpRayleigh::BuildThePhysicsTable()
193{
194// Builds a table of scattering lengths for each material
195
196 if (thePhysicsTable) return;
197
198 const G4MaterialTable* theMaterialTable=
199 G4Material::GetMaterialTable();
200 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
201
202 // create a new physics table
203
204 thePhysicsTable = new G4PhysicsTable(numOfMaterials);
205
206 // loop for materials
207
208 for (G4int i=0 ; i < numOfMaterials; i++)
209 {
210 G4PhysicsOrderedFreeVector* ScatteringLengths =
211 new G4PhysicsOrderedFreeVector();
212
213 G4MaterialPropertiesTable *aMaterialPropertiesTable =
214 (*theMaterialTable)[i]->GetMaterialPropertiesTable();
215
216 if(aMaterialPropertiesTable){
217
218 G4MaterialPropertyVector* AttenuationLengthVector =
219 aMaterialPropertiesTable->GetProperty("RAYLEIGH");
220
221 if(!AttenuationLengthVector){
222
223 if ((*theMaterialTable)[i]->GetName() == "Water")
224 {
225 // Call utility routine to Generate
226 // Rayleigh Scattering Lengths
227
228 DefaultWater = true;
229
230 ScatteringLengths =
231 RayleighAttenuationLengthGenerator(aMaterialPropertiesTable);
232 }
233 }
234 }
235
236 thePhysicsTable->insertAt(i,ScatteringLengths);
237 }
238}
239
240// GetMeanFreePath()
241// -----------------
242//
243G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack,
244 G4double ,
245 G4ForceCondition* )
246{
247 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
248 const G4Material* aMaterial = aTrack.GetMaterial();
249
250 G4double thePhotonMomentum = aParticle->GetTotalMomentum();
251
252 G4double AttenuationLength = DBL_MAX;
253
254 if (aMaterial->GetName() == "Water" && DefaultWater){
255
256 G4bool isOutRange;
257
258 AttenuationLength =
259 (*thePhysicsTable)(aMaterial->GetIndex())->
260 GetValue(thePhotonMomentum, isOutRange);
261 }
262 else {
263
264 G4MaterialPropertiesTable* aMaterialPropertyTable =
265 aMaterial->GetMaterialPropertiesTable();
266
267 if(aMaterialPropertyTable){
268 G4MaterialPropertyVector* AttenuationLengthVector =
269 aMaterialPropertyTable->GetProperty("RAYLEIGH");
270 if(AttenuationLengthVector){
271 AttenuationLength = AttenuationLengthVector ->
272 GetProperty(thePhotonMomentum);
273 }
274 else{
275// G4cout << "No Rayleigh scattering length specified" << G4endl;
276 }
277 }
278 else{
279// G4cout << "No Rayleigh scattering length specified" << G4endl;
280 }
281 }
282
283 return AttenuationLength;
284}
285
286// RayleighAttenuationLengthGenerator()
287// ------------------------------------
288// Private method to compute Rayleigh Scattering Lengths (for water)
289//
290G4PhysicsOrderedFreeVector*
291G4OpRayleigh::RayleighAttenuationLengthGenerator(G4MaterialPropertiesTable *aMPT)
292{
293 // Physical Constants
294
295 // isothermal compressibility of water
296 G4double betat = 7.658e-23*m3/MeV;
297
298 // K Boltzman
299 G4double kboltz = 8.61739e-11*MeV/kelvin;
300
301 // Temperature of water is 10 degrees celsius
302 // conversion to kelvin:
303 // TCelsius = TKelvin - 273.15 => 273.15 + 10 = 283.15
304 G4double temp = 283.15*kelvin;
305
306 // Retrieve vectors for refraction index
307 // and photon momentum from the material properties table
308
309 G4MaterialPropertyVector* Rindex = aMPT->GetProperty("RINDEX");
310
311 G4double refsq;
312 G4double e;
313 G4double xlambda;
314 G4double c1, c2, c3, c4;
315 G4double Dist;
316 G4double refraction_index;
317
318 G4PhysicsOrderedFreeVector *RayleighScatteringLengths =
319 new G4PhysicsOrderedFreeVector();
320
321 if (Rindex ) {
322
323 Rindex->ResetIterator();
324
325 while (++(*Rindex)) {
326
327 e = (Rindex->GetPhotonMomentum());
328
329 refraction_index = Rindex->GetProperty();
330 refsq = refraction_index*refraction_index;
331 xlambda = h_Planck*c_light/e;
332
333 if (verboseLevel>0) {
334 G4cout << Rindex->GetPhotonMomentum() << " MeV\t";
335 G4cout << xlambda << " mm\t";
336 }
337
338 c1 = 1 / (6.0 * pi);
339 c2 = std::pow((2.0 * pi / xlambda), 4);
340 c3 = std::pow( ( (refsq - 1.0) * (refsq + 2.0) / 3.0 ), 2);
341 c4 = betat * temp * kboltz;
342
343 Dist = 1.0 / (c1*c2*c3*c4);
344
345 if (verboseLevel>0) {
346 G4cout << Dist << " mm" << G4endl;
347 }
348 RayleighScatteringLengths->
349 InsertValues(Rindex->GetPhotonMomentum(), Dist);
350 }
351
352 }
353
354 return RayleighScatteringLengths;
355}
Note: See TracBrowser for help on using the repository browser.