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

Last change on this file was 1340, checked in by garnier, 15 years ago

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