Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
83 using namespace std;
84 
85 G4eBremsstrahlungModel::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 
105 G4eBremsstrahlungModel::~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 
117 void 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 
126 void 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 
166 G4double 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 
278 G4double 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 
380 G4double 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 
406 G4double 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 
484 G4double 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 
591 G4double 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 
614 G4DataVector* 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 
643 void 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 
865 const 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 
897 G4double 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 
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4