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

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

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