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

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

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