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