Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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