source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreComptonModel.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: 12.6 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: G4LivermoreComptonModel.cc,v 1.7 2009/06/10 13:32:36 mantero Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29//
30// Author: Sebastien Inserti
31// 30 October 2008
32//
33// History:
34// --------
35// 18 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
36// - apply internal high-energy limit only in constructor
37// - do not apply low-energy limit (default is 0)
38// - remove GetMeanFreePath method and table
39// - added protection against numerical problem in energy sampling
40// - use G4ElementSelector
41
42#include "G4LivermoreComptonModel.hh"
43
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45
46using namespace std;
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50G4LivermoreComptonModel::G4LivermoreComptonModel(const G4ParticleDefinition*,
51 const G4String& nam)
52 :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),
53 scatterFunctionData(0),crossSectionHandler(0)
54{
55 lowEnergyLimit = 250 * eV;
56 highEnergyLimit = 100 * GeV;
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 Compton model 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
79G4LivermoreComptonModel::~G4LivermoreComptonModel()
80{
81 if (crossSectionHandler) delete crossSectionHandler;
82 if (scatterFunctionData) delete scatterFunctionData;
83}
84
85//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
87void G4LivermoreComptonModel::Initialise(const G4ParticleDefinition* particle,
88 const G4DataVector& cuts)
89{
90 if (verboseLevel > 3)
91 G4cout << "Calling G4LivermoreComptonModel::Initialise()" << G4endl;
92
93 if (crossSectionHandler)
94 {
95 crossSectionHandler->Clear();
96 delete crossSectionHandler;
97 }
98
99 // Reading of data files - all materials are read
100
101 crossSectionHandler = new G4CrossSectionHandler;
102 crossSectionHandler->Clear();
103 G4String crossSectionFile = "comp/ce-cs-";
104 crossSectionHandler->LoadData(crossSectionFile);
105
106 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
107 G4String scatterFile = "comp/ce-sf-";
108 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
109 scatterFunctionData->LoadData(scatterFile);
110
111 // For Doppler broadening
112 shellData.SetOccupancyData();
113 G4String file = "/doppler/shell-doppler";
114 shellData.LoadData(file);
115
116 if (verboseLevel > 2)
117 G4cout << "Loaded cross section files for Livermore Compton model" << G4endl;
118
119 InitialiseElementSelectors(particle,cuts);
120
121 if( verboseLevel>0 ) {
122 G4cout << "Livermore Compton model is initialized " << G4endl
123 << "Energy range: "
124 << LowEnergyLimit() / eV << " eV - "
125 << HighEnergyLimit() / GeV << " GeV"
126 << G4endl;
127 }
128 //
129 if(isInitialised) return;
130 fParticleChange = GetParticleChangeForGamma();
131 isInitialised = true;
132}
133
134//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135
136G4double G4LivermoreComptonModel::ComputeCrossSectionPerAtom(
137 const G4ParticleDefinition*,
138 G4double GammaEnergy,
139 G4double Z, G4double,
140 G4double, G4double)
141{
142 if (verboseLevel > 3)
143 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModel" << G4endl;
144
145 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
146
147 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
148 return cs;
149}
150
151//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152
153void G4LivermoreComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
154 const G4MaterialCutsCouple* couple,
155 const G4DynamicParticle* aDynamicGamma,
156 G4double, G4double)
157{
158
159 // The scattered gamma energy is sampled according to Klein - Nishina formula.
160 // then accepted or rejected depending on the Scattering Function multiplied
161 // by factor from Klein - Nishina formula.
162 // Expression of the angular distribution as Klein Nishina
163 // angular and energy distribution and Scattering fuctions is taken from
164 // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
165 // Phys. Res. B 101 (1995). Method of sampling with form factors is different
166 // data are interpolated while in the article they are fitted.
167 // Reference to the article is from J. Stepanek New Photon, Positron
168 // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
169 // TeV (draft).
170 // The random number techniques of Butcher & Messel are used
171 // (Nucl Phys 20(1960),15).
172
173 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
174
175 if (verboseLevel > 3) {
176 G4cout << "G4LivermoreComptonModel::SampleSecondaries() E(MeV)= "
177 << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
178 << G4endl;
179 }
180
181 // low-energy gamma is absorpted by this process
182 if (photonEnergy0 <= lowEnergyLimit)
183 {
184 fParticleChange->ProposeTrackStatus(fStopAndKill);
185 fParticleChange->SetProposedKineticEnergy(0.);
186 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
187 return ;
188 }
189
190 G4double e0m = photonEnergy0 / electron_mass_c2 ;
191 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
192
193 // Select randomly one element in the current material
194 // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
195 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
196 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
197 G4int Z = (G4int)elm->GetZ();
198
199 G4double epsilon0 = 1. / (1. + 2. * e0m);
200 G4double epsilon0Sq = epsilon0 * epsilon0;
201 G4double alpha1 = -std::log(epsilon0);
202 G4double alpha2 = 0.5 * (1. - epsilon0Sq);
203
204 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
205
206 // Sample the energy of the scattered photon
207 G4double epsilon;
208 G4double epsilonSq;
209 G4double oneCosT;
210 G4double sinT2;
211 G4double gReject;
212
213 do
214 {
215 if ( alpha1/(alpha1+alpha2) > G4UniformRand())
216 {
217 // std::pow(epsilon0,G4UniformRand())
218 epsilon = std::exp(-alpha1 * G4UniformRand());
219 epsilonSq = epsilon * epsilon;
220 }
221 else
222 {
223 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand();
224 epsilon = std::sqrt(epsilonSq);
225 }
226
227 oneCosT = (1. - epsilon) / ( epsilon * e0m);
228 sinT2 = oneCosT * (2. - oneCosT);
229 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
230 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
231 gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
232
233 } while(gReject < G4UniformRand()*Z);
234
235 G4double cosTheta = 1. - oneCosT;
236 G4double sinTheta = std::sqrt (sinT2);
237 G4double phi = twopi * G4UniformRand() ;
238 G4double dirx = sinTheta * std::cos(phi);
239 G4double diry = sinTheta * std::sin(phi);
240 G4double dirz = cosTheta ;
241
242 // Doppler broadening - Method based on:
243 // Y. Namito, S. Ban and H. Hirayama,
244 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon
245 // into the EGS4 Code", NIM A 349, pp. 489-494, 1994
246
247 // Maximum number of sampling iterations
248 G4int maxDopplerIterations = 1000;
249 G4double bindingE = 0.;
250 G4double photonEoriginal = epsilon * photonEnergy0;
251 G4double photonE = -1.;
252 G4int iteration = 0;
253 G4double eMax = photonEnergy0;
254 do
255 {
256 iteration++;
257 // Select shell based on shell occupancy
258 G4int shell = shellData.SelectRandomShell(Z);
259 bindingE = shellData.BindingEnergy(Z,shell);
260
261 eMax = photonEnergy0 - bindingE;
262
263 // Randomly sample bound electron momentum
264 // (memento: the data set is in Atomic Units)
265 G4double pSample = profileData.RandomSelectMomentum(Z,shell);
266 // Rescale from atomic units
267 G4double pDoppler = pSample * fine_structure_const;
268 G4double pDoppler2 = pDoppler * pDoppler;
269 G4double var2 = 1. + oneCosT * e0m;
270 G4double var3 = var2*var2 - pDoppler2;
271 G4double var4 = var2 - pDoppler2 * cosTheta;
272 G4double var = var4*var4 - var3 + pDoppler2 * var3;
273 if (var > 0.)
274 {
275 G4double varSqrt = std::sqrt(var);
276 G4double scale = photonEnergy0 / var3;
277 // Random select either root
278 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
279 else photonE = (var4 + varSqrt) * scale;
280 }
281 else
282 {
283 photonE = -1.;
284 }
285 } while ( iteration <= maxDopplerIterations &&
286 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
287
288 // End of recalculation of photon energy with Doppler broadening
289 // Revert to original if maximum number of iterations threshold has been reached
290
291 if (iteration >= maxDopplerIterations)
292 {
293 photonE = photonEoriginal;
294 bindingE = 0.;
295 }
296
297 // Update G4VParticleChange for the scattered photon
298
299 G4ThreeVector photonDirection1(dirx,diry,dirz);
300 photonDirection1.rotateUz(photonDirection0);
301 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
302
303 G4double photonEnergy1 = photonE;
304
305 if (photonEnergy1 > 0.)
306 {
307 fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
308 }
309 else
310 {
311 photonEnergy1 = 0.;
312 fParticleChange->SetProposedKineticEnergy(0.) ;
313 fParticleChange->ProposeTrackStatus(fStopAndKill);
314 }
315
316 // Kinematics of the scattered electron
317 G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
318
319 // protection against negative final energy: no e- is created
320 if(eKineticEnergy < 0.0) {
321 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0 - photonEnergy1);
322 return;
323 }
324 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
325
326 G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2;
327 G4double electronP2 = electronE*electronE - electron_mass_c2*electron_mass_c2;
328 G4double sinThetaE = -1.;
329 G4double cosThetaE = 0.;
330 if (electronP2 > 0.)
331 {
332 cosThetaE = (eTotalEnergy + photonEnergy1 )* (1. - epsilon) / std::sqrt(electronP2);
333 sinThetaE = -1. * sqrt(1. - cosThetaE * cosThetaE);
334 }
335
336 G4double eDirX = sinThetaE * std::cos(phi);
337 G4double eDirY = sinThetaE * std::sin(phi);
338 G4double eDirZ = cosThetaE;
339
340 G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
341 eDirection.rotateUz(photonDirection0);
342
343 // SI - The range test has been removed wrt original G4LowEnergyCompton class
344
345 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
346
347 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
348 eDirection,eKineticEnergy) ;
349 fvect->push_back(dp);
350}
351
Note: See TracBrowser for help on using the repository browser.