source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreGammaConversionModel.cc@ 1058

Last change on this file since 1058 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
File size: 12.0 KB
RevLine 
[968]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//
[1055]26// $Id: G4LivermoreGammaConversionModel.cc,v 1.6 2009/05/02 09:14:43 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
[968]28//
[1055]29//
30// Author: Sebastien Inserti
31// 30 October 2008
32//
33// History:
34// --------
35// 12 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// - use CLHEP electron mass for low-enegry limit
39// - remove MeanFreePath method and table
[968]40
[1055]41
[968]42#include "G4LivermoreGammaConversionModel.hh"
43
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45
46using namespace std;
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50G4LivermoreGammaConversionModel::G4LivermoreGammaConversionModel(const G4ParticleDefinition*,
[1055]51 const G4String& nam)
52 :G4VEmModel(nam),smallEnergy(2.*MeV),isInitialised(false),
53 crossSectionHandler(0),meanFreePathTable(0)
[968]54{
[1055]55 lowEnergyLimit = 2.0*electron_mass_c2;
[968]56 highEnergyLimit = 100 * GeV;
[1055]57 SetHighEnergyLimit(highEnergyLimit);
58
[968]59 verboseLevel= 0;
60 // Verbosity scale:
61 // 0 = nothing
62 // 1 = warning for energy non-conservation
63 // 2 = details of energy budget
64 // 3 = calculation of cross sections, file openings, sampling of atoms
65 // 4 = entering in methods
66
[1055]67 if(verboseLevel > 0) {
68 G4cout << "Livermore Gamma conversion is constructed " << G4endl
69 << "Energy range: "
70 << lowEnergyLimit / MeV << " MeV - "
71 << highEnergyLimit / GeV << " GeV"
72 << G4endl;
73 }
[968]74}
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78G4LivermoreGammaConversionModel::~G4LivermoreGammaConversionModel()
79{
[1055]80 if (crossSectionHandler) delete crossSectionHandler;
[968]81}
82
83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84
[1055]85void
86G4LivermoreGammaConversionModel::Initialise(const G4ParticleDefinition*,
87 const G4DataVector&)
[968]88{
89 if (verboseLevel > 3)
90 G4cout << "Calling G4LivermoreGammaConversionModel::Initialise()" << G4endl;
91
[1055]92 if (crossSectionHandler)
[968]93 {
[1055]94 crossSectionHandler->Clear();
95 delete crossSectionHandler;
[968]96 }
97
98 // Read data tables for all materials
99
100 crossSectionHandler = new G4CrossSectionHandler();
[1055]101 crossSectionHandler->Initialise(0,lowEnergyLimit,100.*GeV,400);
[968]102 G4String crossSectionFile = "pair/pp-cs-";
103 crossSectionHandler->LoadData(crossSectionFile);
104
105 //
106
107 if (verboseLevel > 2)
108 G4cout << "Loaded cross section files for PenelopeGammaConversion" << G4endl;
109
[1055]110 if (verboseLevel > 0) {
111 G4cout << "Livermore Gamma Conversion model is initialized " << G4endl
112 << "Energy range: "
113 << LowEnergyLimit() / MeV << " MeV - "
114 << HighEnergyLimit() / GeV << " GeV"
115 << G4endl;
116 }
[968]117
118 if(isInitialised) return;
[1055]119 fParticleChange = GetParticleChangeForGamma();
120 isInitialised = true;
121}
[968]122
123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124
[1055]125G4double
126G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
127 G4double GammaEnergy,
128 G4double Z, G4double,
129 G4double, G4double)
[968]130{
[1055]131 if (verboseLevel > 3) {
132 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel"
133 << G4endl;
134 }
135 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0;
[968]136
137 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
138 return cs;
139}
140
141//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142
143void G4LivermoreGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
144 const G4MaterialCutsCouple* couple,
145 const G4DynamicParticle* aDynamicGamma,
146 G4double,
147 G4double)
148{
149
150// The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
151// cross sections with Coulomb correction. A modified version of the random
152// number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
153
154// Note 1 : Effects due to the breakdown of the Born approximation at low
155// energy are ignored.
156// Note 2 : The differential cross section implicitly takes account of
157// pair creation in both nuclear and atomic electron fields. However triplet
158// prodution is not generated.
159
160 if (verboseLevel > 3)
161 G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModel" << G4endl;
162
163 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
164 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
165
166 G4double epsilon ;
167 G4double epsilon0 = electron_mass_c2 / photonEnergy ;
168
169 // Do it fast if photon energy < 2. MeV
170 if (photonEnergy < smallEnergy )
171 {
172 epsilon = epsilon0 + (0.5 - epsilon0) * G4UniformRand();
173 }
174 else
175 {
176 // Select randomly one element in the current material
[1055]177 //const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
178 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
179 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
[968]180
181 if (element == 0)
182 {
[1055]183 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - element = 0"
184 << G4endl;
185 return;
[968]186 }
187 G4IonisParamElm* ionisation = element->GetIonisation();
[1055]188 if (ionisation == 0)
[968]189 {
[1055]190 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0"
191 << G4endl;
192 return;
[968]193 }
194
195 // Extract Coulomb factor for this Element
196 G4double fZ = 8. * (ionisation->GetlogZ3());
197 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
198
199 // Limits of the screening variable
200 G4double screenFactor = 136. * epsilon0 / (element->GetIonisation()->GetZ3()) ;
201 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
202 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
203
204 // Limits of the energy sampling
205 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
206 G4double epsilonMin = std::max(epsilon0,epsilon1);
207 G4double epsilonRange = 0.5 - epsilonMin ;
208
209 // Sample the energy rate of the created electron (or positron)
210 G4double screen;
211 G4double gReject ;
212
213 G4double f10 = ScreenFunction1(screenMin) - fZ;
214 G4double f20 = ScreenFunction2(screenMin) - fZ;
215 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
216 G4double normF2 = std::max(1.5 * f20,0.);
217
218 do {
219 if (normF1 / (normF1 + normF2) > G4UniformRand() )
220 {
221 epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.3333) ;
222 screen = screenFactor / (epsilon * (1. - epsilon));
223 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
224 }
225 else
226 {
227 epsilon = epsilonMin + epsilonRange * G4UniformRand();
228 screen = screenFactor / (epsilon * (1 - epsilon));
229 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
230 }
231 } while ( gReject < G4UniformRand() );
232
233 } // End of epsilon sampling
234
235 // Fix charges randomly
236
237 G4double electronTotEnergy;
238 G4double positronTotEnergy;
239
240 if (CLHEP::RandBit::shootBit())
241 {
242 electronTotEnergy = (1. - epsilon) * photonEnergy;
243 positronTotEnergy = epsilon * photonEnergy;
244 }
245 else
246 {
247 positronTotEnergy = (1. - epsilon) * photonEnergy;
248 electronTotEnergy = epsilon * photonEnergy;
249 }
250
251 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
252 // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
253 // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
254
255 G4double u;
256 const G4double a1 = 0.625;
257 G4double a2 = 3. * a1;
258 // G4double d = 27. ;
259
260 // if (9. / (9. + d) > G4UniformRand())
261 if (0.25 > G4UniformRand())
262 {
263 u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ;
264 }
265 else
266 {
267 u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ;
268 }
269
270 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
271 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
272 G4double phi = twopi * G4UniformRand();
273
274 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
275 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
276
277
278 // Kinematics of the created pair:
279 // the electron and positron are assumed to have a symetric angular
280 // distribution with respect to the Z axis along the parent photon
281
282 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
283
[1055]284 // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
[968]285
286 G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
287 electronDirection.rotateUz(photonDirection);
288
289 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
290 electronDirection,
291 electronKineEnergy);
292
293 // The e+ is always created (even with kinetic energy = 0) for further annihilation
294 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
295
[1055]296 // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
[968]297
298 G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
299 positronDirection.rotateUz(photonDirection);
300
301 // Create G4DynamicParticle object for the particle2
302 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
303 positronDirection, positronKineEnergy);
304 // Fill output vector
305 fvect->push_back(particle1);
306 fvect->push_back(particle2);
307
308 // kill incident photon
309 fParticleChange->SetProposedKineticEnergy(0.);
310 fParticleChange->ProposeTrackStatus(fStopAndKill);
311
312}
313
314//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315
316G4double G4LivermoreGammaConversionModel::ScreenFunction1(G4double screenVariable)
317{
318 // Compute the value of the screening function 3*phi1 - phi2
319
320 G4double value;
321
322 if (screenVariable > 1.)
323 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
324 else
325 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
326
327 return value;
328}
329
330//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
331
332G4double G4LivermoreGammaConversionModel::ScreenFunction2(G4double screenVariable)
333{
334 // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
335
336 G4double value;
337
338 if (screenVariable > 1.)
339 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
340 else
341 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
342
343 return value;
344}
345
Note: See TracBrowser for help on using the repository browser.