Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eBremsstrahlungRelModel.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: G4eBremsstrahlungRelModel
34 //
35 // Author: Andreas Schaelicke
36 //
37 // Creation date: 12.08.2008
38 //
39 // Modifications:
40 //
41 // 13.11.08 add SetLPMflag and SetLPMconstant methods
42 // 13.11.08 change default LPMconstant value
43 // 13.10.10 add angular distributon interface (VI)
44 //
45 // Main References:
46 // Y.-S.Tsai, Rev. Mod. Phys. 46 (1974) 815; Rev. Mod. Phys. 49 (1977) 421.
47 // S.Klein, Rev. Mod. Phys. 71 (1999) 1501.
48 // T.Stanev et.al., Phys. Rev. D25 (1982) 1291.
49 // M.L.Ter-Mikaelian, High-energy Electromagnetic Processes in Condensed Media, Wiley, 1972.
50 //
51 // -------------------------------------------------------------------
52 //
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 
56 #include "G4eBremsstrahlungRelModel.hh"
57 #include "G4PhysicalConstants.hh"
58 #include "G4SystemOfUnits.hh"
59 #include "G4Electron.hh"
60 #include "G4Positron.hh"
61 #include "G4Gamma.hh"
62 #include "Randomize.hh"
63 #include "G4Material.hh"
64 #include "G4Element.hh"
65 #include "G4ElementVector.hh"
66 #include "G4ProductionCutsTable.hh"
67 #include "G4ParticleChangeForLoss.hh"
68 #include "G4LossTableManager.hh"
69 #include "G4ModifiedTsai.hh"
70 #include "G4DipBustGenerator.hh"
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
74 const G4double G4eBremsstrahlungRelModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
75  0.5917, 0.7628, 0.8983, 0.9801 };
76 const G4double G4eBremsstrahlungRelModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
77  0.1813, 0.1569, 0.1112, 0.0506 };
78 const G4double G4eBremsstrahlungRelModel::Fel_light[] = {0., 5.31 , 4.79 , 4.74 , 4.71} ;
79 const G4double G4eBremsstrahlungRelModel::Finel_light[] = {0., 6.144 , 5.621 , 5.805 , 5.924} ;
80 
81 using namespace std;
82 
83 G4eBremsstrahlungRelModel::G4eBremsstrahlungRelModel(const G4ParticleDefinition* p,
84  const G4String& nam)
85  : G4VEmModel(nam),
86  particle(0),
87  bremFactor(fine_structure_const*classic_electr_radius*classic_electr_radius*16./3.),
88  isElectron(true),
89  fMigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi),
90  fLPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)*0.5),
91  fXiLPM(0), fPhiLPM(0), fGLPM(0),
92  use_completescreening(false),isInitialised(false)
93 {
94  fParticleChange = 0;
95  theGamma = G4Gamma::Gamma();
96 
97  lowKinEnergy = 0.1*GeV;
98  SetLowEnergyLimit(lowKinEnergy);
99 
100  nist = G4NistManager::Instance();
101 
102  SetLPMFlag(true);
103  //SetAngularDistribution(new G4ModifiedTsai());
104  SetAngularDistribution(new G4DipBustGenerator());
105 
106  particleMass = kinEnergy = totalEnergy = currentZ = z13 = z23 = lnZ = Fel
107  = Finel = fCoulomb = fMax = densityFactor = densityCorr = lpmEnergy
108  = xiLPM = phiLPM = gLPM = klpm = kp = 0.0;
109 
110  energyThresholdLPM = 1.e39;
111 
112  InitialiseConstants();
113  if(p) { SetParticle(p); }
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 
118 void G4eBremsstrahlungRelModel::InitialiseConstants()
119 {
120  facFel = log(184.15);
121  facFinel = log(1194.);
122 
123  preS1 = 1./(184.15*184.15);
124  logTwo = log(2.);
125 }
126 
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
128 
129 G4eBremsstrahlungRelModel::~G4eBremsstrahlungRelModel()
130 {
131 }
132 
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134 
135 void G4eBremsstrahlungRelModel::SetParticle(const G4ParticleDefinition* p)
136 {
137  particle = p;
138  particleMass = p->GetPDGMass();
139  if(p == G4Electron::Electron()) { isElectron = true; }
140  else { isElectron = false;}
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144 
145 void G4eBremsstrahlungRelModel::SetupForMaterial(const G4ParticleDefinition*,
146  const G4Material* mat,
147  G4double kineticEnergy)
148 {
149  densityFactor = mat->GetElectronDensity()*fMigdalConstant;
150  lpmEnergy = mat->GetRadlen()*fLPMconstant;
151 
152  // Threshold for LPM effect (i.e. below which LPM hidden by density effect)
153  if (LPMFlag()) {
154  energyThresholdLPM=sqrt(densityFactor)*lpmEnergy;
155  } else {
156  energyThresholdLPM=1.e39; // i.e. do not use LPM effect
157  }
158  // calculate threshold for density effect
159  kinEnergy = kineticEnergy;
160  totalEnergy = kineticEnergy + particleMass;
161  densityCorr = densityFactor*totalEnergy*totalEnergy;
162 
163  // define critical gamma energies (important for integration/dicing)
164  klpm=totalEnergy*totalEnergy/lpmEnergy;
165  kp=sqrt(densityCorr);
166 
167 }
168 
169 
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171 
172 void G4eBremsstrahlungRelModel::Initialise(const G4ParticleDefinition* p,
173  const G4DataVector& cuts)
174 {
175  if(p) { SetParticle(p); }
176 
177  lowKinEnergy = LowEnergyLimit();
178 
179  currentZ = 0.;
180 
181  InitialiseElementSelectors(p, cuts);
182 
183  if(isInitialised) { return; }
184  fParticleChange = GetParticleChangeForLoss();
185  isInitialised = true;
186 }
187 
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189 
190 G4double G4eBremsstrahlungRelModel::ComputeDEDXPerVolume(
191  const G4Material* material,
192  const G4ParticleDefinition* p,
193  G4double kineticEnergy,
194  G4double cutEnergy)
195 {
196  if(!particle) { SetParticle(p); }
197  if(kineticEnergy < lowKinEnergy) { return 0.0; }
198  G4double cut = std::min(cutEnergy, kineticEnergy);
199  if(cut == 0.0) { return 0.0; }
200 
201  SetupForMaterial(particle, material,kineticEnergy);
202 
203  const G4ElementVector* theElementVector = material->GetElementVector();
204  const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
205 
206  G4double dedx = 0.0;
207 
208  // loop for elements in the material
209  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
210 
211  G4VEmModel::SetCurrentElement((*theElementVector)[i]);
212  SetCurrentElement((*theElementVector)[i]->GetZ());
213 
214  dedx += theAtomicNumDensityVector[i]*currentZ*currentZ*ComputeBremLoss(cut);
215  }
216  dedx *= bremFactor;
217 
218 
219  return dedx;
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
223 
224 G4double G4eBremsstrahlungRelModel::ComputeBremLoss(G4double cut)
225 {
226  G4double loss = 0.0;
227 
228  // number of intervals and integration step
229  G4double vcut = cut/totalEnergy;
230  G4int n = (G4int)(20*vcut) + 3;
231  G4double delta = vcut/G4double(n);
232 
233  G4double e0 = 0.0;
234  G4double xs;
235 
236  // integration
237  for(G4int l=0; l<n; l++) {
238 
239  for(G4int i=0; i<8; i++) {
240 
241  G4double eg = (e0 + xgi[i]*delta)*totalEnergy;
242 
243  if(totalEnergy > energyThresholdLPM) {
244  xs = ComputeRelDXSectionPerAtom(eg);
245  } else {
246  xs = ComputeDXSectionPerAtom(eg);
247  }
248  loss += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
249  }
250  e0 += delta;
251  }
252 
253  loss *= delta*totalEnergy;
254 
255  return loss;
256 }
257 
258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
259 
260 G4double G4eBremsstrahlungRelModel::ComputeCrossSectionPerAtom(
261  const G4ParticleDefinition* p,
262  G4double kineticEnergy,
263  G4double Z, G4double,
264  G4double cutEnergy,
265  G4double maxEnergy)
266 {
267  if(!particle) { SetParticle(p); }
268  if(kineticEnergy < lowKinEnergy) { return 0.0; }
269  G4double cut = std::min(cutEnergy, kineticEnergy);
270  G4double tmax = std::min(maxEnergy, kineticEnergy);
271 
272  if(cut >= tmax) { return 0.0; }
273 
274  SetCurrentElement(Z);
275 
276  G4double cross = ComputeXSectionPerAtom(cut);
277 
278  // allow partial integration
279  if(tmax < kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
280 
281  cross *= Z*Z*bremFactor;
282 
283  return cross;
284 }
285 
286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
287 
288 
289 G4double G4eBremsstrahlungRelModel::ComputeXSectionPerAtom(G4double cut)
290 {
291  G4double cross = 0.0;
292 
293  // number of intervals and integration step
294  G4double vcut = log(cut/totalEnergy);
295  G4double vmax = log(kinEnergy/totalEnergy);
296  G4int n = (G4int)(0.45*(vmax - vcut)) + 4;
297  // n=1; // integration test
298  G4double delta = (vmax - vcut)/G4double(n);
299 
300  G4double e0 = vcut;
301  G4double xs;
302 
303  // integration
304  for(G4int l=0; l<n; l++) {
305 
306  for(G4int i=0; i<8; i++) {
307 
308  G4double eg = exp(e0 + xgi[i]*delta)*totalEnergy;
309 
310  if(totalEnergy > energyThresholdLPM) {
311  xs = ComputeRelDXSectionPerAtom(eg);
312  } else {
313  xs = ComputeDXSectionPerAtom(eg);
314  }
315  cross += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
316  }
317  e0 += delta;
318  }
319 
320  cross *= delta;
321 
322  return cross;
323 }
324 
325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
326 void G4eBremsstrahlungRelModel::CalcLPMFunctions(G4double k)
327 {
328  // *** calculate lpm variable s & sprime ***
329  // Klein eqs. (78) & (79)
330  G4double sprime = sqrt(0.125*k*lpmEnergy/(totalEnergy*(totalEnergy-k)));
331 
332  G4double s1 = preS1*z23;
333  G4double logS1 = 2./3.*lnZ-2.*facFel;
334  G4double logTS1 = logTwo+logS1;
335 
336  xiLPM = 2.;
337 
338  if (sprime>1)
339  xiLPM = 1.;
340  else if (sprime>sqrt(2.)*s1) {
341  G4double h = log(sprime)/logTS1;
342  xiLPM = 1+h-0.08*(1-h)*(1-sqr(1-h))/logTS1;
343  }
344 
345  G4double s0 = sprime/sqrt(xiLPM);
346 
347  // *** merging with density effect*** should be only necessary in region "close to" kp, e.g. k<100*kp
348  // using Ter-Mikaelian eq. (20.9)
349  G4double k2 = k*k;
350  s0 *= (1 + (densityCorr/k2) );
351 
352  // recalculate Xi using modified s above
353  // Klein eq. (75)
354  xiLPM = 1.;
355  if (s0<=s1) xiLPM = 2.;
356  else if ( (s1<s0) && (s0<=1) ) xiLPM = 1. + log(s0)/logS1;
357 
358 
359  // *** calculate supression functions phi and G ***
360  // Klein eqs. (77)
361  G4double s2=s0*s0;
362  G4double s3=s0*s2;
363  G4double s4=s2*s2;
364 
365  if (s0<0.1) {
366  // high suppression limit
367  phiLPM = 6.*s0 - 18.84955592153876*s2 + 39.47841760435743*s3
368  - 57.69873135166053*s4;
369  gLPM = 37.69911184307752*s2 - 236.8705056261446*s3 + 807.7822389*s4;
370  }
371  else if (s0<1.9516) {
372  // intermediate suppression
373  // using eq.77 approxim. valid s<2.
374  phiLPM = 1.-exp(-6.*s0*(1.+(3.-pi)*s0)
375  +s3/(0.623+0.795*s0+0.658*s2));
376  if (s0<0.415827397755) {
377  // using eq.77 approxim. valid 0.07<s<2
378  G4double psiLPM = 1-exp(-4*s0-8*s2/(1+3.936*s0+4.97*s2-0.05*s3+7.50*s4));
379  gLPM = 3*psiLPM-2*phiLPM;
380  }
381  else {
382  // using alternative parametrisiation
383  G4double pre = -0.16072300849123999 + s0*3.7550300067531581 + s2*-1.7981383069010097
384  + s3*0.67282686077812381 + s4*-0.1207722909879257;
385  gLPM = tanh(pre);
386  }
387  }
388  else {
389  // low suppression limit valid s>2.
390  phiLPM = 1. - 0.0119048/s4;
391  gLPM = 1. - 0.0230655/s4;
392  }
393 
394  // *** make sure suppression is smaller than 1 ***
395  // *** caused by Migdal approximation in xi ***
396  if (xiLPM*phiLPM>1. || s0>0.57) { xiLPM=1./phiLPM; }
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 
401 
402 G4double G4eBremsstrahlungRelModel::ComputeRelDXSectionPerAtom(G4double gammaEnergy)
403 // Ultra relativistic model
404 // only valid for very high energies, but includes LPM suppression
405 // * complete screening
406 {
407  if(gammaEnergy < 0.0) { return 0.0; }
408 
409  G4double y = gammaEnergy/totalEnergy;
410  G4double y2 = y*y*.25;
411  G4double yone2 = (1.-y+2.*y2);
412 
413  // ** form factors complete screening case **
414 
415  // ** calc LPM functions -- include ter-mikaelian merging with density effect **
416  // G4double xiLPM, gLPM, phiLPM; // to be made member variables !!!
417  CalcLPMFunctions(gammaEnergy);
418 
419  G4double mainLPM = xiLPM*(y2 * gLPM + yone2*phiLPM) * ( (Fel-fCoulomb) + Finel/currentZ );
420  G4double secondTerm = (1.-y)/12.*(1.+1./currentZ);
421 
422  G4double cross = mainLPM+secondTerm;
423  return cross;
424 }
425 
426 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
427 
428 G4double G4eBremsstrahlungRelModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
429 // Relativistic model
430 // only valid for high energies (and if LPM suppression does not play a role)
431 // * screening according to thomas-fermi-Model (only valid for Z>5)
432 // * no LPM effect
433 {
434 
435  if(gammaEnergy < 0.0) { return 0.0; }
436 
437  G4double y = gammaEnergy/totalEnergy;
438 
439  G4double main=0.,secondTerm=0.;
440 
441  if (use_completescreening|| currentZ<5) {
442  // ** form factors complete screening case **
443  main = (3./4.*y*y - y + 1.) * ( (Fel-fCoulomb) + Finel/currentZ );
444  secondTerm = (1.-y)/12.*(1.+1./currentZ);
445  }
446  else {
447  // ** intermediate screening using Thomas-Fermi FF from Tsai only valid for Z>=5**
448  G4double dd=100.*electron_mass_c2*y/(totalEnergy-gammaEnergy);
449  G4double gg=dd/z13;
450  G4double eps=dd/z23;
451  G4double phi1=Phi1(gg,currentZ), phi1m2=Phi1M2(gg,currentZ);
452  G4double psi1=Psi1(eps,currentZ), psi1m2=Psi1M2(eps,currentZ);
453 
454  main = (3./4.*y*y - y + 1.) * ( (0.25*phi1-1./3.*lnZ-fCoulomb) + (0.25*psi1-2./3.*lnZ)/currentZ );
455  secondTerm = (1.-y)/8.*(phi1m2+psi1m2/currentZ);
456  }
457  G4double cross = main+secondTerm;
458  return cross;
459 }
460 
461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
462 
463 void G4eBremsstrahlungRelModel::SampleSecondaries(
464  std::vector<G4DynamicParticle*>* vdp,
465  const G4MaterialCutsCouple* couple,
466  const G4DynamicParticle* dp,
467  G4double cutEnergy,
468  G4double maxEnergy)
469 {
470  G4double kineticEnergy = dp->GetKineticEnergy();
471  if(kineticEnergy < lowKinEnergy) { return; }
472  G4double cut = std::min(cutEnergy, kineticEnergy);
473  G4double emax = std::min(maxEnergy, kineticEnergy);
474  if(cut >= emax) { return; }
475 
476  SetupForMaterial(particle, couple->GetMaterial(), kineticEnergy);
477 
478  const G4Element* elm =
479  SelectRandomAtom(couple,particle,kineticEnergy,cut,emax);
480  SetCurrentElement(elm->GetZ());
481 
482  kinEnergy = kineticEnergy;
483  totalEnergy = kineticEnergy + particleMass;
484  densityCorr = densityFactor*totalEnergy*totalEnergy;
485 
486  //G4double fmax= fMax;
487  G4bool highe = true;
488  if(totalEnergy < energyThresholdLPM) { highe = false; }
489 
490  G4double xmin = log(cut*cut + densityCorr);
491  G4double xmax = log(emax*emax + densityCorr);
492  G4double gammaEnergy, f, x;
493 
494  do {
495  x = exp(xmin + G4UniformRand()*(xmax - xmin)) - densityCorr;
496  if(x < 0.0) { x = 0.0; }
497  gammaEnergy = sqrt(x);
498  if(highe) { f = ComputeRelDXSectionPerAtom(gammaEnergy); }
499  else { f = ComputeDXSectionPerAtom(gammaEnergy); }
500 
501  if ( f > fMax ) {
502  G4cout << "### G4eBremsstrahlungRelModel Warning: Majoranta exceeded! "
503  << f << " > " << fMax
504  << " Egamma(MeV)= " << gammaEnergy
505  << " Ee(MeV)= " << kineticEnergy
506  << " " << GetName()
507  << G4endl;
508  }
509 
510  } while (f < fMax*G4UniformRand());
511 
512  //
513  // angles of the emitted gamma. ( Z - axis along the parent particle)
514  // use general interface
515  //
516 
517  G4ThreeVector gammaDirection =
518  GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
519  G4lrint(currentZ),
520  couple->GetMaterial());
521 
522  // create G4DynamicParticle object for the Gamma
523  G4DynamicParticle* gamma = new G4DynamicParticle(theGamma,gammaDirection,
524  gammaEnergy);
525  vdp->push_back(gamma);
526 
527  G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
528  G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
529  - gammaEnergy*gammaDirection).unit();
530 
531  // energy of primary
532  G4double finalE = kineticEnergy - gammaEnergy;
533 
534  // stop tracking and create new secondary instead of primary
535  if(gammaEnergy > SecondaryThreshold()) {
536  fParticleChange->ProposeTrackStatus(fStopAndKill);
537  fParticleChange->SetProposedKineticEnergy(0.0);
538  G4DynamicParticle* el =
539  new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
540  direction, finalE);
541  vdp->push_back(el);
542 
543  // continue tracking
544  } else {
545  fParticleChange->SetProposedMomentumDirection(direction);
546  fParticleChange->SetProposedKineticEnergy(finalE);
547  }
548 }
549 
550 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
551 
552 
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4