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