source: trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc@ 877

Last change on this file since 877 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 31.1 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: G4eBremsstrahlungModel.cc,v 1.39 2007/05/23 08:47:35 vnivanch Exp $
27// GEANT4 tag $Name: $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4eBremsstrahlungModel
35//
36// Author: Vladimir Ivanchenko on base of Laszlo Urban code
37//
38// Creation date: 03.01.2002
39//
40// Modifications:
41//
42// 11-11-02 Fix division by 0 (V.Ivanchenko)
43// 04-12-02 Change G4DynamicParticle constructor in PostStep (V.Ivanchenko)
44// 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
45// 24-01-03 Fix for compounds (V.Ivanchenko)
46// 27-01-03 Make models region aware (V.Ivanchenko)
47// 13-02-03 Add name (V.Ivanchenko)
48// 09-05-03 Fix problem of supression function + optimise sampling (V.Ivanchenko)
49// 20-05-04 Correction to ensure unit independence (L.Urban)
50// 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
51// 03-08-05 Add extra protection at initialisation (V.Ivantchenko)
52// 07-02-06 public function ComputeCrossSectionPerAtom() (mma)
53// 21-03-06 Fix problem of initialisation in case when cuts are not defined (VI)
54// 27-03-06 Fix calculation of fl parameter at low energy (energy loss) (VI)
55// 15-02-07 correct LPMconstant by a factor 2, thanks to G. Depaola (mma)
56//
57// Class Description:
58//
59//
60// -------------------------------------------------------------------
61//
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64
65#include "G4eBremsstrahlungModel.hh"
66#include "G4Electron.hh"
67#include "G4Positron.hh"
68#include "G4Gamma.hh"
69#include "Randomize.hh"
70#include "G4Material.hh"
71#include "G4Element.hh"
72#include "G4ElementVector.hh"
73#include "G4ProductionCutsTable.hh"
74#include "G4DataVector.hh"
75#include "G4ParticleChangeForLoss.hh"
76
77//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78
79using namespace std;
80
81G4eBremsstrahlungModel::G4eBremsstrahlungModel(const G4ParticleDefinition* p,
82 const G4String& nam)
83 : G4VEmModel(nam),
84 particle(0),
85 isElectron(true),
86 probsup(1.0),
87 MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length/pi),
88 LPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)),
89 theLPMflag(true),
90 isInitialised(false)
91{
92 if(p) SetParticle(p);
93 theGamma = G4Gamma::Gamma();
94 minThreshold = 1.0*keV;
95 highKinEnergy= 100.*TeV;
96 lowKinEnergy = 1.0*keV;
97 highEnergyTh = DBL_MAX;
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
102G4eBremsstrahlungModel::~G4eBremsstrahlungModel()
103{
104 size_t n = partialSumSigma.size();
105 if(n > 0) {
106 for(size_t i=0; i<n; i++) {
107 delete partialSumSigma[i];
108 }
109 }
110}
111
112//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113
114void G4eBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p)
115{
116 particle = p;
117 if(p == G4Electron::Electron()) isElectron = true;
118 else isElectron = false;
119}
120
121//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
122
123G4double G4eBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
124 const G4MaterialCutsCouple*)
125{
126 return minThreshold;
127}
128
129//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
130
131void G4eBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
132 const G4DataVector& cuts)
133{
134 if(p) SetParticle(p);
135 highKinEnergy = HighEnergyLimit();
136 lowKinEnergy = LowEnergyLimit();
137 const G4ProductionCutsTable* theCoupleTable=
138 G4ProductionCutsTable::GetProductionCutsTable();
139
140 if(theCoupleTable) {
141 G4int numOfCouples = theCoupleTable->GetTableSize();
142
143 G4int nn = partialSumSigma.size();
144 G4int nc = cuts.size();
145 if(nn > 0) {
146 for (G4int ii=0; ii<nn; ii++){
147 G4DataVector* a=partialSumSigma[ii];
148 if ( a ) delete a;
149 }
150 partialSumSigma.clear();
151 }
152 if(numOfCouples>0) {
153 for (G4int i=0; i<numOfCouples; i++) {
154 G4double cute = DBL_MAX;
155 if(i < nc) cute = cuts[i];
156 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
157 const G4Material* material = couple->GetMaterial();
158 G4DataVector* dv = ComputePartialSumSigma(material, 0.5*highKinEnergy,
159 std::min(cute, 0.25*highKinEnergy));
160 partialSumSigma.push_back(dv);
161 }
162 }
163 }
164 if(isInitialised) return;
165
166 if(pParticleChange)
167 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
168 else
169 fParticleChange = new G4ParticleChangeForLoss();
170
171 isInitialised = true;
172}
173
174//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
175
176G4double G4eBremsstrahlungModel::ComputeDEDXPerVolume(
177 const G4Material* material,
178 const G4ParticleDefinition* p,
179 G4double kineticEnergy,
180 G4double cutEnergy)
181{
182 if(!particle) SetParticle(p);
183 if(kineticEnergy < lowKinEnergy) return 0.0;
184
185 const G4double thigh = 100.*GeV;
186
187 G4double cut = std::min(cutEnergy, kineticEnergy);
188
189 G4double rate, loss;
190 const G4double factorHigh = 36./(1450.*GeV);
191 const G4double coef1 = -0.5;
192 const G4double coef2 = 2./9.;
193
194 const G4ElementVector* theElementVector = material->GetElementVector();
195 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
196
197 G4double totalEnergy = kineticEnergy + electron_mass_c2;
198 G4double dedx = 0.0;
199
200 // loop for elements in the material
201 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
202
203 G4double Z = (*theElementVector)[i]->GetZ();
204 G4double natom = theAtomicNumDensityVector[i];
205
206 // loss for MinKinEnergy<KineticEnergy<=100 GeV
207 if (kineticEnergy <= thigh) {
208
209 // x = log(totalEnergy/electron_mass_c2);
210 loss = ComputeBremLoss(Z, kineticEnergy, cut) ;
211 if (!isElectron) loss *= PositronCorrFactorLoss(Z, kineticEnergy, cut);
212
213 // extrapolation for KineticEnergy>100 GeV
214 } else {
215
216 // G4double xhigh = log(thigh/electron_mass_c2);
217 G4double cuthigh = thigh*0.5;
218
219 if (cut < thigh) {
220
221 loss = ComputeBremLoss(Z, thigh, cut) ;
222 if (!isElectron) loss *= PositronCorrFactorLoss(Z, thigh, cut) ;
223 rate = cut/totalEnergy;
224 loss *= (1. + coef1*rate + coef2*rate*rate);
225 rate = cut/thigh;
226 loss /= (1.+coef1*rate+coef2*rate*rate);
227
228 } else {
229
230 loss = ComputeBremLoss(Z, thigh, cuthigh) ;
231 if (!isElectron) loss *= PositronCorrFactorLoss(Z, thigh, cuthigh) ;
232 rate = cut/totalEnergy;
233 loss *= (1. + coef1*rate + coef2*rate*rate);
234 loss *= cut*factorHigh;
235 }
236 }
237 loss *= natom;
238
239 G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
240 * (material->GetElectronDensity()) ;
241
242 // now compute the correction due to the supression(s)
243 G4double kmin = 1.*eV;
244 G4double kmax = cut;
245
246 if (kmax > kmin) {
247
248 G4double floss = 0.;
249 G4int nmax = 100;
250
251 G4double vmin=log(kmin);
252 G4double vmax=log(kmax) ;
253 G4int nn = (G4int)(nmax*(vmax-vmin)/(log(highKinEnergy)-vmin)) ;
254 G4double u,fac,c,v,dv ;
255 if(nn > 0) {
256
257 dv = (vmax-vmin)/nn ;
258 v = vmin-dv ;
259
260 for(G4int n=0; n<=nn; n++) {
261
262 v += dv;
263 u = exp(v);
264 fac = u*SupressionFunction(material,kineticEnergy,u);
265 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
266 if ((n==0)||(n==nn)) c=0.5;
267 else c=1. ;
268 fac *= c ;
269 floss += fac ;
270 }
271 floss *=dv/(kmax-kmin);
272
273 } else {
274 floss = 1.;
275 }
276 if(floss > 1.) floss = 1.;
277 // correct the loss
278 loss *= floss;
279 }
280 dedx += loss;
281 }
282 if(dedx < 0.) dedx = 0.;
283 return dedx;
284}
285
286//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
287
288G4double G4eBremsstrahlungModel::ComputeBremLoss(G4double Z, G4double T,
289 G4double Cut)
290
291// compute loss due to soft brems
292{
293 static const G4double beta=1.0, ksi=2.0;
294 static const G4double clossh = 0.254 , closslow = 1./3. , alosslow = 1. ;
295 static const G4double Tlim= 10.*MeV ;
296
297 static const G4double xlim = 1.2 ;
298 static const G4int NZ = 8 ;
299 static const G4int Nloss = 11 ;
300 static const G4double ZZ[NZ] =
301 {2.,4.,6.,14.,26.,50.,82.,92.};
302 static const G4double coefloss[NZ][Nloss] = {
303 // Z=2
304 { 0.98916, 0.47564, -0.2505, -0.45186, 0.14462,
305 0.21307, -0.013738, -0.045689, -0.0042914, 0.0034429,
306 0.00064189},
307
308 // Z=4
309 { 1.0626, 0.37662, -0.23646, -0.45188, 0.14295,
310 0.22906, -0.011041, -0.051398, -0.0055123, 0.0039919,
311 0.00078003},
312 // Z=6
313 { 1.0954, 0.315, -0.24011, -0.43849, 0.15017,
314 0.23001, -0.012846, -0.052555, -0.0055114, 0.0041283,
315 0.00080318},
316
317 // Z=14
318 { 1.1649, 0.18976, -0.24972, -0.30124, 0.1555,
319 0.13565, -0.024765, -0.027047, -0.00059821, 0.0019373,
320 0.00027647},
321
322 // Z=26
323 { 1.2261, 0.14272, -0.25672, -0.28407, 0.13874,
324 0.13586, -0.020562, -0.026722, -0.00089557, 0.0018665,
325 0.00026981},
326
327 // Z=50
328 { 1.3147, 0.020049, -0.35543, -0.13927, 0.17666,
329 0.073746, -0.036076, -0.013407, 0.0025727, 0.00084005,
330 -1.4082e-05},
331
332 // Z=82
333 { 1.3986, -0.10586, -0.49187, -0.0048846, 0.23621,
334 0.031652, -0.052938, -0.0076639, 0.0048181, 0.00056486,
335 -0.00011995},
336
337 // Z=92
338 { 1.4217, -0.116, -0.55497, -0.044075, 0.27506,
339 0.081364, -0.058143, -0.023402, 0.0031322, 0.0020201,
340 0.00017519}
341
342 } ;
343 static G4double aaa = 0.414;
344 static G4double bbb = 0.345;
345 static G4double ccc = 0.460;
346
347 G4int iz = 0;
348 G4double delz = 1.e6;
349 for (G4int ii=0; ii<NZ; ii++)
350 {
351 G4double dz = std::abs(Z-ZZ[ii]);
352 if(dz < delz) {
353 iz = ii;
354 delz = dz;
355 }
356 }
357
358 G4double xx = log10(T/MeV);
359 G4double fl = 1.;
360
361 if (xx <= xlim)
362 {
363 xx /= xlim;
364 G4double yy = 1.0;
365 fl = 0.0;
366 for (G4int j=0; j<Nloss; j++) {
367 fl += yy+coefloss[iz][j];
368 yy *= xx;
369 }
370 if (fl < 0.00001) fl = 0.00001;
371 else if (fl > 1.0) fl = 1.0;
372 }
373
374 G4double loss;
375 G4double E = T+electron_mass_c2 ;
376
377 loss = Z*(Z+ksi)*E*E/(T+E)*exp(beta*log(Cut/T))*(2.-clossh*exp(log(Z)/4.));
378 if (T <= Tlim) loss /= exp(closslow*log(Tlim/T));
379 if( T <= Cut) loss *= exp(alosslow*log(T/Cut));
380 // correction
381 loss *= (aaa+bbb*T/Tlim)/(1.+ccc*T/Tlim);
382 loss *= fl;
383 loss /= Avogadro;
384
385 return loss;
386}
387
388//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389
390G4double G4eBremsstrahlungModel::PositronCorrFactorLoss(G4double Z,
391 G4double kineticEnergy, G4double cut)
392
393//calculates the correction factor for the energy loss due to bremsstrahlung for positrons
394//the same correction is in the (discrete) bremsstrahlung
395
396{
397 static const G4double K = 132.9416*eV ;
398 static const G4double a1=4.15e-1, a3=2.10e-3, a5=54.0e-5 ;
399
400 G4double x = log(kineticEnergy/(K*Z*Z)), x2 = x*x, x3 = x2*x;
401 G4double eta = 0.5+atan(a1*x+a3*x3+a5*x3*x2)/pi;
402 G4double e0 = cut/kineticEnergy;
403
404 G4double factor = 0.0;
405 if (e0 < 1.0) {
406 factor=log(1.-e0)/eta;
407 factor=exp(factor);
408 }
409 factor = eta*(1.-factor)/e0;
410
411 return factor;
412}
413
414//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
415
416G4double G4eBremsstrahlungModel::CrossSectionPerVolume(
417 const G4Material* material,
418 const G4ParticleDefinition* p,
419 G4double kineticEnergy,
420 G4double cutEnergy,
421 G4double maxEnergy)
422{
423 if(!particle) SetParticle(p);
424 G4double cross = 0.0;
425 G4double tmax = min(maxEnergy, kineticEnergy);
426 G4double cut = max(cutEnergy, minThreshold);
427 if(cut >= tmax) return cross;
428
429 const G4ElementVector* theElementVector = material->GetElementVector();
430 const G4double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
431 G4double dum=0.;
432
433 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
434
435 cross += theAtomNumDensityVector[i] * ComputeCrossSectionPerAtom(p,
436 kineticEnergy, (*theElementVector)[i]->GetZ(), dum, cut);
437 if (tmax < kineticEnergy) {
438 cross -= theAtomNumDensityVector[i] * ComputeCrossSectionPerAtom(p,
439 kineticEnergy, (*theElementVector)[i]->GetZ(), dum, tmax);
440 }
441 }
442
443 // now compute the correction due to the supression(s)
444
445 G4double kmax = tmax;
446 G4double kmin = cut;
447
448 G4double totalEnergy = kineticEnergy+electron_mass_c2 ;
449 G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
450 *(material->GetElectronDensity());
451
452 G4double fsig = 0.;
453 G4int nmax = 100;
454 G4double vmin=log(kmin);
455 G4double vmax=log(kmax) ;
456 G4int nn = (G4int)(nmax*(vmax-vmin)/(log(highKinEnergy)-vmin));
457 G4double u,fac,c,v,dv,y ;
458 if(nn > 0) {
459
460 dv = (vmax-vmin)/nn ;
461 v = vmin-dv ;
462 for(G4int n=0; n<=nn; n++) {
463
464 v += dv;
465 u = exp(v);
466 fac = SupressionFunction(material, kineticEnergy, u);
467 y = u/kmax;
468 fac *= (4.-4.*y+3.*y*y)/3.;
469 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
470
471 if ((n==0)||(n==nn)) c=0.5;
472 else c=1. ;
473
474 fac *= c;
475 fsig += fac;
476 }
477 y = kmin/kmax ;
478 fsig *=dv/(-4.*log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
479
480 } else {
481
482 fsig = 1.;
483 }
484 if (fsig > 1.) fsig = 1.;
485
486 // correct the cross section
487 cross *= fsig;
488
489 return cross;
490}
491
492//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
493
494G4double G4eBremsstrahlungModel::ComputeCrossSectionPerAtom(
495 const G4ParticleDefinition*,
496 G4double kineticEnergy,
497 G4double Z, G4double,
498 G4double cut, G4double)
499
500// Calculates the cross section per atom in GEANT4 internal units.
501//
502
503{
504 G4double cross = 0.0 ;
505 if ( kineticEnergy < 1*keV || kineticEnergy < cut) return cross;
506
507 static const G4double ksi=2.0, alfa=1.00;
508 static const G4double csigh = 0.127, csiglow = 0.25, asiglow = 0.020*MeV ;
509 static const G4double Tlim = 10.*MeV ;
510
511 static const G4double xlim = 1.2 ;
512 static const G4int NZ = 8 ;
513 static const G4int Nsig = 11 ;
514 static const G4double ZZ[NZ] =
515 {2.,4.,6.,14.,26.,50.,82.,92.} ;
516 static const G4double coefsig[NZ][Nsig] = {
517 // Z=2
518 { 0.4638, 0.37748, 0.32249, -0.060362, -0.065004,
519 -0.033457, -0.004583, 0.011954, 0.0030404, -0.0010077,
520 -0.00028131},
521
522 // Z=4
523 { 0.50008, 0.33483, 0.34364, -0.086262, -0.055361,
524 -0.028168, -0.0056172, 0.011129, 0.0027528, -0.00092265,
525 -0.00024348},
526
527 // Z=6
528 { 0.51587, 0.31095, 0.34996, -0.11623, -0.056167,
529 -0.0087154, 0.00053943, 0.0054092, 0.00077685, -0.00039635,
530 -6.7818e-05},
531
532 // Z=14
533 { 0.55058, 0.25629, 0.35854, -0.080656, -0.054308,
534 -0.049933, -0.00064246, 0.016597, 0.0021789, -0.001327,
535 -0.00025983},
536
537 // Z=26
538 { 0.5791, 0.26152, 0.38953, -0.17104, -0.099172,
539 0.024596, 0.023718, -0.0039205, -0.0036658, 0.00041749,
540 0.00023408},
541
542 // Z=50
543 { 0.62085, 0.27045, 0.39073, -0.37916, -0.18878,
544 0.23905, 0.095028, -0.068744, -0.023809, 0.0062408,
545 0.0020407},
546
547 // Z=82
548 { 0.66053, 0.24513, 0.35404, -0.47275, -0.22837,
549 0.35647, 0.13203, -0.1049, -0.034851, 0.0095046,
550 0.0030535},
551
552 // Z=92
553 { 0.67143, 0.23079, 0.32256, -0.46248, -0.20013,
554 0.3506, 0.11779, -0.1024, -0.032013, 0.0092279,
555 0.0028592}
556
557 } ;
558
559 G4int iz = 0 ;
560 G4double delz = 1.e6 ;
561 for (G4int ii=0; ii<NZ; ii++)
562 {
563 G4double absdelz = std::abs(Z-ZZ[ii]);
564 if(absdelz < delz)
565 {
566 iz = ii ;
567 delz = absdelz;
568 }
569 }
570
571 G4double xx = log10(kineticEnergy/MeV) ;
572 G4double fs = 1. ;
573
574 if (xx <= xlim) {
575
576 fs = coefsig[iz][Nsig-1] ;
577 for (G4int j=Nsig-2; j>=0; j--) {
578
579 fs = fs*xx+coefsig[iz][j] ;
580 }
581 if(fs < 0.) fs = 0.;
582 }
583
584 cross = Z*(Z+ksi)*(1.-csigh*exp(log(Z)/4.))*pow(log(kineticEnergy/cut),alfa);
585
586 if (kineticEnergy <= Tlim)
587 cross *= exp(csiglow*log(Tlim/kineticEnergy))
588 *(1.+asiglow/(sqrt(Z)*kineticEnergy));
589
590 if (!isElectron)
591 cross *= PositronCorrFactorSigma(Z, kineticEnergy, cut);
592
593 cross *= fs/Avogadro ;
594
595 if (cross < 0.) cross = 0.;
596 return cross;
597}
598
599//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
600
601G4double G4eBremsstrahlungModel::PositronCorrFactorSigma( G4double Z,
602 G4double kineticEnergy, G4double cut)
603
604//Calculates the correction factor for the total cross section of the positron
605// bremsstrahl.
606// Eta is the ratio of positron to electron energy loss by bremstrahlung.
607// A parametrized formula from L. Urban is used to estimate eta. It is a fit to
608// the results of L. Kim & al: Phys Rev. A33,3002 (1986)
609
610{
611 static const G4double K = 132.9416*eV;
612 static const G4double a1 = 4.15e-1, a3 = 2.10e-3, a5 = 54.0e-5;
613
614 G4double x = log(kineticEnergy/(K*Z*Z));
615 G4double x2 = x*x;
616 G4double x3 = x2*x;
617 G4double eta = 0.5 + atan(a1*x + a3*x3 + a5*x3*x2)/pi ;
618 G4double alfa = (1. - eta)/eta;
619 return eta*pow((1. - cut/kineticEnergy), alfa);
620}
621
622//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
623
624G4DataVector* G4eBremsstrahlungModel::ComputePartialSumSigma(
625 const G4Material* material,
626 G4double kineticEnergy,
627 G4double cut)
628
629// Build the table of cross section per element.
630//The table is built for MATERIALS.
631// This table is used by DoIt to select randomly an element in the material.
632{
633 G4int nElements = material->GetNumberOfElements();
634 const G4ElementVector* theElementVector = material->GetElementVector();
635 const G4double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
636 G4double dum = 0.;
637
638 G4DataVector* dv = new G4DataVector();
639
640 G4double cross = 0.0;
641
642 for (G4int i=0; i<nElements; i++ ) {
643
644 cross += theAtomNumDensityVector[i] * ComputeCrossSectionPerAtom( particle,
645 kineticEnergy, (*theElementVector)[i]->GetZ(), dum, cut);
646 dv->push_back(cross);
647 }
648 return dv;
649}
650
651//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
652
653void G4eBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
654 const G4MaterialCutsCouple* couple,
655 const G4DynamicParticle* dp,
656 G4double tmin,
657 G4double maxEnergy)
658// The emitted gamma energy is sampled using a parametrized formula
659// from L. Urban.
660// This parametrization is derived from :
661// cross-section values of Seltzer and Berger for electron energies
662// 1 keV - 10 GeV,
663// screened Bethe Heilter differential cross section above 10 GeV,
664// Migdal corrections in both case.
665// Seltzer & Berger: Nim B 12:95 (1985)
666// Nelson, Hirayama & Rogers: Technical report 265 SLAC (1985)
667// Migdal: Phys Rev 103:1811 (1956); Messel & Crawford: Pergamon Press (1970)
668//
669// A modified version of the random number techniques of Butcher&Messel is used
670// (Nuc Phys 20(1960),15).
671{
672 G4double kineticEnergy = dp->GetKineticEnergy();
673 G4double tmax = min(maxEnergy, kineticEnergy);
674 if(tmin >= tmax) return;
675
676//
677// GEANT4 internal units.
678//
679 static const G4double
680 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
681 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
682 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
683
684 static const G4double
685 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
686 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
687 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
688
689 static const G4double
690 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
691 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
692 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
693
694 static const G4double
695 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
696 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
697 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
698
699 static const G4double tlow = 1.*MeV;
700
701 G4double gammaEnergy;
702 G4bool LPMOK = false;
703 const G4Material* material = couple->GetMaterial();
704
705 // select randomly one element constituing the material
706 const G4Element* anElement = SelectRandomAtom(couple);
707
708 // Extract Z factors for this Element
709 G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3());
710 G4double FZ = lnZ* (4.- 0.55*lnZ);
711 G4double ZZ = anElement->GetIonisation()->GetZZ3();
712
713 // limits of the energy sampling
714 G4double totalEnergy = kineticEnergy + electron_mass_c2;
715 G4ThreeVector direction = dp->GetMomentumDirection();
716 G4double xmin = tmin/kineticEnergy;
717 G4double xmax = tmax/kineticEnergy;
718 G4double kappa = 0.0;
719 if(xmax >= 1.) xmax = 1.;
720 else kappa = log(xmax)/log(xmin);
721 G4double epsilmin = tmin/totalEnergy;
722 G4double epsilmax = tmax/totalEnergy;
723
724 // Migdal factor
725 G4double MigdalFactor = (material->GetElectronDensity())*MigdalConstant
726 / (epsilmax*epsilmax);
727
728 G4double x, epsil, greject, migdal, grejmax, q;
729 G4double U = log(kineticEnergy/electron_mass_c2);
730 G4double U2 = U*U;
731
732 // precalculated parameters
733 G4double ah, bh;
734 G4double screenfac = 0.0;
735
736 if (kineticEnergy > tlow) {
737
738 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12);
739 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22);
740 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
741
742 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12);
743 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22);
744 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
745
746 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
747 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
748
749 // limit of the screening variable
750 screenfac =
751 136.*electron_mass_c2/((anElement->GetIonisation()->GetZ3())*totalEnergy);
752 G4double screenmin = screenfac*epsilmin/(1.-epsilmin);
753
754 // Compute the maximum of the rejection function
755 G4double F1 = max(ScreenFunction1(screenmin) - FZ ,0.);
756 G4double F2 = max(ScreenFunction2(screenmin) - FZ ,0.);
757 grejmax = (F1 - epsilmin* (F1*ah - bh*epsilmin*F2))/(42.392 - FZ);
758
759 } else {
760
761 G4double al0 = al00 + ZZ* (al01 + ZZ* al02);
762 G4double al1 = al10 + ZZ* (al11 + ZZ* al12);
763 G4double al2 = al20 + ZZ* (al21 + ZZ* al22);
764
765 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02);
766 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12);
767 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
768
769 ah = al0 + al1*U + al2*U2;
770 bh = bl0 + bl1*U + bl2*U2;
771
772 // Compute the maximum of the rejection function
773 grejmax = max(1. + xmin* (ah + bh*xmin), 1.+ah+bh);
774 G4double xm = -ah/(2.*bh);
775 if ( xmin < xm && xm < xmax) grejmax = max(grejmax, 1.+ xm* (ah + bh*xm));
776 }
777
778 //
779 // sample the energy rate of the emitted gamma for e- kin energy > 1 MeV
780 //
781
782 do {
783 if (kineticEnergy > tlow) {
784 do {
785 q = G4UniformRand();
786 x = pow(xmin, q + kappa*(1.0 - q));
787 epsil = x*kineticEnergy/totalEnergy;
788 G4double screenvar = screenfac*epsil/(1.0-epsil);
789 G4double F1 = max(ScreenFunction1(screenvar) - FZ ,0.);
790 G4double F2 = max(ScreenFunction2(screenvar) - FZ ,0.);
791 migdal = (1. + MigdalFactor)/(1. + MigdalFactor/(x*x));
792 greject = migdal*(F1 - epsil* (ah*F1 - bh*epsil*F2))/(42.392 - FZ);
793 /*
794 if ( greject > grejmax ) {
795 G4cout << "### G4eBremsstrahlungModel Warning: Majoranta exceeded! "
796 << greject << " > " << grejmax
797 << " x= " << x
798 << " e= " << kineticEnergy
799 << G4endl;
800 }
801 */
802 } while( greject < G4UniformRand()*grejmax );
803
804 } else {
805
806 do {
807 q = G4UniformRand();
808 x = pow(xmin, q + kappa*(1.0 - q));
809 migdal = (1. + MigdalFactor)/(1. + MigdalFactor/(x*x));
810 greject = migdal*(1. + x* (ah + bh*x));
811 /*
812 if ( greject > grejmax ) {
813 G4cout << "### G4eBremsstrahlungModel Warning: Majoranta exceeded! "
814 << greject << " > " << grejmax
815 << " x= " << x
816 << " e= " << kineticEnergy
817 << G4endl;
818 }
819 */
820 } while( greject < G4UniformRand()*grejmax );
821 }
822 /*
823 if(x > 0.999) {
824 G4cout << "### G4eBremsstrahlungModel Warning: e= " << kineticEnergy
825 << " tlow= " << tlow
826 << " x= " << x
827 << " greject= " << greject
828 << " grejmax= " << grejmax
829 << " migdal= " << migdal
830 << G4endl;
831 // if(x >= 1.0) G4Exception("X=1");
832 }
833 */
834 gammaEnergy = x*kineticEnergy;
835
836 if (theLPMflag) {
837 // take into account the supression due to the LPM effect
838 if (G4UniformRand() <= SupressionFunction(material,kineticEnergy,
839 gammaEnergy))
840 LPMOK = true;
841 }
842 else LPMOK = true;
843
844 } while (!LPMOK);
845
846 //
847 // angles of the emitted gamma. ( Z - axis along the parent particle)
848 //
849 // universal distribution suggested by L. Urban
850 // (Geant3 manual (1993) Phys211),
851 // derived from Tsai distribution (Rev Mod Phys 49,421(1977))
852
853 G4double u;
854 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
855
856 if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1;
857 else u = - log(G4UniformRand()*G4UniformRand())/a2;
858
859 G4double theta = u*electron_mass_c2/totalEnergy;
860
861 G4double sint = sin(theta);
862
863 G4double phi = twopi * G4UniformRand() ;
864
865 G4ThreeVector gammaDirection(sint*cos(phi),sint*sin(phi), cos(theta));
866 gammaDirection.rotateUz(direction);
867
868 // create G4DynamicParticle object for the Gamma
869 G4DynamicParticle* g = new G4DynamicParticle(theGamma,gammaDirection,
870 gammaEnergy);
871 vdp->push_back(g);
872
873 G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
874 G4ThreeVector dir = totMomentum*direction - gammaEnergy*gammaDirection;
875 direction = dir.unit();
876
877 // energy of primary
878 G4double finalE = kineticEnergy - gammaEnergy;
879
880 // stop tracking and create new secondary instead of primary
881 if(gammaEnergy > highEnergyTh) {
882 fParticleChange->ProposeTrackStatus(fStopAndKill);
883 fParticleChange->SetProposedKineticEnergy(0.0);
884 G4DynamicParticle* el =
885 new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
886 direction, finalE);
887 vdp->push_back(el);
888
889 // continue tracking
890 } else {
891 fParticleChange->SetProposedMomentumDirection(direction);
892 fParticleChange->SetProposedKineticEnergy(finalE);
893 }
894}
895
896//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897
898const G4Element* G4eBremsstrahlungModel::SelectRandomAtom(
899 const G4MaterialCutsCouple* couple)
900{
901 // select randomly 1 element within the material
902
903 const G4Material* material = couple->GetMaterial();
904 G4int nElements = material->GetNumberOfElements();
905 const G4ElementVector* theElementVector = material->GetElementVector();
906
907 const G4Element* elm = 0;
908
909 if(1 < nElements) {
910
911 G4DataVector* dv = partialSumSigma[couple->GetIndex()];
912 G4double rval = G4UniformRand()*((*dv)[nElements-1]);
913
914 for (G4int i=0; i<nElements; i++) {
915 if (rval <= (*dv)[i]) elm = (*theElementVector)[i];
916 }
917 if(!elm) {
918 G4cout << "G4eBremsstrahlungModel::SelectRandomAtom: Warning -"
919 << " no elements found in "
920 << material->GetName()
921 << G4endl;
922 elm = (*theElementVector)[0];
923 }
924 } else elm = (*theElementVector)[0];
925
926 SetCurrentElement(elm);
927 return elm;
928}
929
930//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
931
932G4double G4eBremsstrahlungModel::SupressionFunction(const G4Material* material,
933 G4double kineticEnergy, G4double gammaEnergy)
934{
935 // supression due to the LPM effect+polarisation of the medium/
936 // supression due to the polarisation alone
937
938
939 G4double totEnergy = kineticEnergy+electron_mass_c2 ;
940 G4double totEnergySquare = totEnergy*totEnergy ;
941
942 G4double LPMEnergy = LPMconstant*(material->GetRadlen()) ;
943
944 G4double gammaEnergySquare = gammaEnergy*gammaEnergy ;
945
946 G4double electronDensity = material->GetElectronDensity();
947
948 G4double sp = gammaEnergySquare/
949 (gammaEnergySquare+MigdalConstant*totEnergySquare*electronDensity);
950
951 G4double supr = 1.0;
952
953 if (theLPMflag) {
954
955 G4double s2lpm = LPMEnergy*gammaEnergy/totEnergySquare;
956
957 if (s2lpm < 1.) {
958
959 G4double LPMgEnergyLimit = totEnergySquare/LPMEnergy ;
960 G4double LPMgEnergyLimit2 = LPMgEnergyLimit*LPMgEnergyLimit;
961 G4double splim = LPMgEnergyLimit2/
962 (LPMgEnergyLimit2+MigdalConstant*totEnergySquare*electronDensity);
963 G4double w = 1.+1./splim ;
964
965 if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp);
966 else w = s2lpm*(1.+1./sp);
967
968 supr = (sqrt(w*w+4.*s2lpm)-w)/(sqrt(w*w+4.)-w) ;
969 supr /= sp;
970 }
971
972 }
973 return supr;
974}
975
976//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
977
978
Note: See TracBrowser for help on using the repository browser.