source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreRayleighModel.cc@ 1197

Last change on this file since 1197 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 7.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// $Id: G4LivermoreRayleighModel.cc,v 1.8 2009/09/23 16:54:06 flongo Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29// Author: Sebastien Inserti
30// 30 October 2008
31//
32// History:
33// --------
34// 18 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
35// - apply internal high-energy limit only in constructor
36// - do not apply low-energy limit (default is 0)
37// - remove GetMeanFreePath method and table
38// - remove initialisation of element selector
39// - use G4ElementSelector
40
41#include "G4LivermoreRayleighModel.hh"
42
43//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44
45using namespace std;
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
49G4LivermoreRayleighModel::G4LivermoreRayleighModel(const G4ParticleDefinition*,
50 const G4String& nam)
51 :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),
52 formFactorData(0),crossSectionHandler(0)
53{
54 lowEnergyLimit = 250 * eV;
55 highEnergyLimit = 100 * GeV;
56
57 // SetLowEnergyLimit(lowEnergyLimit);
58 SetHighEnergyLimit(highEnergyLimit);
59 //
60 verboseLevel= 0;
61 // Verbosity scale:
62 // 0 = nothing
63 // 1 = warning for energy non-conservation
64 // 2 = details of energy budget
65 // 3 = calculation of cross sections, file openings, sampling of atoms
66 // 4 = entering in methods
67
68 if(verboseLevel > 0) {
69 G4cout << "Livermore Rayleigh is constructed " << G4endl
70 << "Energy range: "
71 << lowEnergyLimit / eV << " eV - "
72 << highEnergyLimit / GeV << " GeV"
73 << G4endl;
74 }
75}
76
77//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78
79G4LivermoreRayleighModel::~G4LivermoreRayleighModel()
80{
81 if (crossSectionHandler) delete crossSectionHandler;
82 if (formFactorData) delete formFactorData;
83}
84
85//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
87void G4LivermoreRayleighModel::Initialise(const G4ParticleDefinition* particle,
88 const G4DataVector& cuts)
89{
90 if (verboseLevel > 3)
91 G4cout << "Calling G4LivermoreRayleighModel::Initialise()" << G4endl;
92
93 if (crossSectionHandler)
94 {
95 crossSectionHandler->Clear();
96 delete crossSectionHandler;
97 }
98
99 // Data are read for all materials
100
101 crossSectionHandler = new G4CrossSectionHandler;
102 crossSectionHandler->Clear();
103 G4String crossSectionFile = "rayl/re-cs-";
104 crossSectionHandler->LoadData(crossSectionFile);
105
106 G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation;
107 G4String formFactorFile = "rayl/re-ff-";
108 formFactorData = new G4CompositeEMDataSet(ffInterpolation,1.,1.);
109 formFactorData->LoadData(formFactorFile);
110
111 InitialiseElementSelectors(particle,cuts);
112
113 //
114 if (verboseLevel > 2)
115 G4cout << "Loaded cross section files for Livermore Rayleigh model" << G4endl;
116
117 if (verboseLevel > 0) {
118 G4cout << "Livermore Rayleigh model is initialized " << G4endl
119 << "Energy range: "
120 << LowEnergyLimit() / eV << " eV - "
121 << HighEnergyLimit() / GeV << " GeV"
122 << G4endl;
123 }
124
125 if(isInitialised) return;
126 fParticleChange = GetParticleChangeForGamma();
127 isInitialised = true;
128
129}
130
131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132
133G4double G4LivermoreRayleighModel::ComputeCrossSectionPerAtom(
134 const G4ParticleDefinition*,
135 G4double GammaEnergy,
136 G4double Z, G4double,
137 G4double, G4double)
138{
139 if (verboseLevel > 3)
140 G4cout << "Calling CrossSectionPerAtom() of G4LivermoreRayleighModel" << G4endl;
141
142 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit)
143 return 0.0;
144
145 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
146 return cs;
147}
148
149//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150
151void G4LivermoreRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
152 const G4MaterialCutsCouple* couple,
153 const G4DynamicParticle* aDynamicGamma,
154 G4double,
155 G4double)
156{
157 if (verboseLevel > 3)
158 G4cout << "Calling SampleSecondaries() of G4LivermoreRayleighModel" << G4endl;
159
160 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
161
162 // absorption of low-energy gamma
163 if (photonEnergy0 <= lowEnergyLimit)
164 {
165 fParticleChange->ProposeTrackStatus(fStopAndKill);
166 fParticleChange->SetProposedKineticEnergy(0.);
167 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
168 return ;
169 }
170
171 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
172
173 // Select randomly one element in the current material
174 // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
175 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
176 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
177 G4int Z = (G4int)elm->GetZ();
178
179 // Sample the angle of the scattered photon
180
181 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
182
183 G4double gReject,x,dataFormFactor;
184 G4double randomFormFactor;
185 G4double cosTheta;
186 G4double sinTheta;
187 G4double fcostheta;
188
189 do
190 {
191 do
192 {
193 cosTheta = 2. * G4UniformRand() - 1.;
194 fcostheta = ( 1. + cosTheta*cosTheta)/2.;
195 } while (fcostheta < G4UniformRand());
196
197 if (photonEnergy0 > 5)
198 {
199 cosTheta = 1.;
200 }
201
202 G4double sinThetaHalf = std::sqrt((1. - cosTheta) / 2.);
203 x = sinThetaHalf / (wlPhoton/cm);
204 if (x > 1.e+005)
205 {
206 dataFormFactor = formFactorData->FindValue(x,Z-1);
207 }
208 else
209 {
210 dataFormFactor = formFactorData->FindValue(0.,Z-1);
211 }
212 randomFormFactor = G4UniformRand() * Z * Z;
213 sinTheta = std::sqrt(1. - cosTheta*cosTheta);
214 gReject = dataFormFactor * dataFormFactor;
215
216 } while( gReject < randomFormFactor);
217
218 // Scattered photon angles. ( Z - axis along the parent photon)
219 G4double phi = twopi * G4UniformRand() ;
220 G4double dirX = sinTheta*std::cos(phi);
221 G4double dirY = sinTheta*std::sin(phi);
222 G4double dirZ = cosTheta;
223
224 // Update G4VParticleChange for the scattered photon
225 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
226 photonDirection1.rotateUz(photonDirection0);
227 fParticleChange->ProposeMomentumDirection(photonDirection1);
228
229 fParticleChange->SetProposedKineticEnergy(photonEnergy0);
230}
231
232
Note: See TracBrowser for help on using the repository browser.