Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LowEnergyRayleigh.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 //
28 // $Id$
29 // GEANT4 tag $Name: geant4-09-01-ref-00 $
30 //
31 // Author: A. Forti
32 // Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
33 //
34 // History:
35 // --------
36 // Added Livermore data table construction methods A. Forti
37 // Added BuildMeanFreePath A. Forti
38 // Added PostStepDoIt A. Forti
39 // Added SelectRandomAtom A. Forti
40 // Added map of the elements A.Forti
41 // 24.04.01 V.Ivanchenko remove RogueWave
42 // 11.08.2001 MGP - Major revision according to a design iteration
43 // 06.10.2001 MGP - Added strategy to test range for secondary generation
44 // 05.06.2002 F.Longo and G.Depaola - bug fixed in angular distribution
45 // 20.10.2002 G. Depaola - Change sampling method of theta
46 // 22.01.2003 V.Ivanchenko - Cut per region
47 // 24.04.2003 V.Ivanchenko - Cut per region mfpt
48 //
49 // --------------------------------------------------------------------
50 
51 #include "G4LowEnergyRayleigh.hh"
52 #include "Randomize.hh"
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4ParticleDefinition.hh"
56 #include "G4Track.hh"
57 #include "G4Step.hh"
58 #include "G4ForceCondition.hh"
59 #include "G4Gamma.hh"
60 #include "G4Electron.hh"
61 #include "G4DynamicParticle.hh"
62 #include "G4VParticleChange.hh"
63 #include "G4ThreeVector.hh"
66 #include "G4RDVEMDataSet.hh"
68 #include "G4RDVDataSetAlgorithm.hh"
70 
71 #include "G4MaterialCutsCouple.hh"
72 
74  : G4VDiscreteProcess(processName),
75  lowEnergyLimit(250*eV),
76  highEnergyLimit(100*GeV),
77  intrinsicLowEnergyLimit(10*eV),
78  intrinsicHighEnergyLimit(100*GeV)
79 {
80  if (lowEnergyLimit < intrinsicLowEnergyLimit ||
81  highEnergyLimit > intrinsicHighEnergyLimit)
82  {
83  G4Exception("G4LowEnergyRayleigh::G4LowEnergyRayleigh()",
84  "OutOfRange", FatalException,
85  "Energy limit outside intrinsic process validity range!");
86  }
87 
88  crossSectionHandler = new G4RDCrossSectionHandler();
89 
90  G4RDVDataSetAlgorithm* ffInterpolation = new G4RDLogLogInterpolation;
91  G4String formFactorFile = "rayl/re-ff-";
92  formFactorData = new G4RDCompositeEMDataSet(ffInterpolation,1.,1.);
93  formFactorData->LoadData(formFactorFile);
94 
95  meanFreePathTable = 0;
96 
97  if (verboseLevel > 0)
98  {
99  G4cout << GetProcessName() << " is created " << G4endl
100  << "Energy range: "
101  << lowEnergyLimit / keV << " keV - "
102  << highEnergyLimit / GeV << " GeV"
103  << G4endl;
104  }
105 }
106 
108 {
109  delete meanFreePathTable;
110  delete crossSectionHandler;
111  delete formFactorData;
112 }
113 
115 {
116 
117  crossSectionHandler->Clear();
118  G4String crossSectionFile = "rayl/re-cs-";
119  crossSectionHandler->LoadData(crossSectionFile);
120 
121  delete meanFreePathTable;
122  meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
123 }
124 
126  const G4Step& aStep)
127 {
128 
129  aParticleChange.Initialize(aTrack);
130 
131  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
132  G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
133 
134  if (photonEnergy0 <= lowEnergyLimit)
135  {
139  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
140  }
141 
142  // G4double e0m = photonEnergy0 / electron_mass_c2 ;
143  G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
144 
145  // Select randomly one element in the current material
146  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
147  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
148 
149  // Sample the angle of the scattered photon
150 
151  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
152 
153  G4double gReject,x,dataFormFactor;
154  G4double randomFormFactor;
155  G4double cosTheta;
156  G4double sinTheta;
157  G4double fcostheta;
158 
159  do
160  {
161  do
162  {
163  cosTheta = 2. * G4UniformRand() - 1.;
164  fcostheta = ( 1. + cosTheta*cosTheta)/2.;
165  } while (fcostheta < G4UniformRand());
166 
167  G4double sinThetaHalf = std::sqrt((1. - cosTheta) / 2.);
168  x = sinThetaHalf / (wlPhoton/cm);
169  if (x > 1.e+005)
170  dataFormFactor = formFactorData->FindValue(x,Z-1);
171  else
172  dataFormFactor = formFactorData->FindValue(0.,Z-1);
173  randomFormFactor = G4UniformRand() * Z * Z;
174  sinTheta = std::sqrt(1. - cosTheta*cosTheta);
175  gReject = dataFormFactor * dataFormFactor;
176 
177  } while( gReject < randomFormFactor);
178 
179  // Scattered photon angles. ( Z - axis along the parent photon)
180  G4double phi = twopi * G4UniformRand() ;
181  G4double dirX = sinTheta*std::cos(phi);
182  G4double dirY = sinTheta*std::sin(phi);
183  G4double dirZ = cosTheta;
184 
185  // Update G4VParticleChange for the scattered photon
186  G4ThreeVector photonDirection1(dirX, dirY, dirZ);
187 
188  photonDirection1.rotateUz(photonDirection0);
189  aParticleChange.ProposeEnergy(photonEnergy0);
190  aParticleChange.ProposeMomentumDirection(photonDirection1);
191 
193 
194  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
195 }
196 
198 {
199  return ( &particle == G4Gamma::Gamma() );
200 }
201 
203  G4double, // previousStepSize
205 {
207  G4double energy = photon->GetKineticEnergy();
208  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
209  size_t materialIndex = couple->GetIndex();
210 
211  G4double meanFreePath;
212  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
213  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
214  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
215  return meanFreePath;
216 }