Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eBremParametrizedModel.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 // GEANT4 tag $Name: geant4-09-04 $
28 //
29 // -------------------------------------------------------------------
30 //
31 // GEANT4 Class file
32 //
33 //
34 // File name: G4eBremParametrizedModel
35 //
36 // Author: Andreas Schaelicke
37 //
38 // Creation date: 06.04.2011
39 //
40 // Modifications:
41 //
42 // Main References:
43 // - based on G4eBremsstrahlungModel and G4eBremsstrahlungRelModel
44 // -------------------------------------------------------------------
45 //
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 #include "G4eBremParametrizedModel.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4Electron.hh"
53 #include "G4Positron.hh"
54 #include "G4Gamma.hh"
55 #include "Randomize.hh"
56 #include "G4Material.hh"
57 #include "G4Element.hh"
58 #include "G4ElementVector.hh"
59 #include "G4ProductionCutsTable.hh"
60 #include "G4ParticleChangeForLoss.hh"
61 #include "G4LossTableManager.hh"
62 #include "G4ModifiedTsai.hh"
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 
66 const G4double G4eBremParametrizedModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
67  0.5917, 0.7628, 0.8983, 0.9801 };
68 const G4double G4eBremParametrizedModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
69  0.1813, 0.1569, 0.1112, 0.0506 };
70 
71 
72 using namespace std;
73 
74 G4eBremParametrizedModel::G4eBremParametrizedModel(const G4ParticleDefinition* p,
75  const G4String& nam)
76  : G4VEmModel(nam),
77  particle(0),
78  isElectron(true),
79  fMigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi),
80  bremFactor(fine_structure_const*classic_electr_radius*classic_electr_radius*16./3.),
81  isInitialised(false)
82 {
83  theGamma = G4Gamma::Gamma();
84 
85  minThreshold = 0.1*keV;
86  lowKinEnergy = 10.*MeV;
87  SetLowEnergyLimit(lowKinEnergy);
88 
89  nist = G4NistManager::Instance();
90 
91  SetAngularDistribution(new G4ModifiedTsai());
92 
93  particleMass = kinEnergy = totalEnergy = currentZ = z13 = z23 = lnZ = Fel = Finel
94  = densityFactor = densityCorr = fMax = fCoulomb = 0.;
95 
96  InitialiseConstants();
97  if(p) { SetParticle(p); }
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
102 void G4eBremParametrizedModel::InitialiseConstants()
103 {
104  facFel = log(184.15);
105  facFinel = log(1194.);
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109 
110 G4eBremParametrizedModel::~G4eBremParametrizedModel()
111 {
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115 
116 void G4eBremParametrizedModel::SetParticle(const G4ParticleDefinition* p)
117 {
118  particle = p;
119  particleMass = p->GetPDGMass();
120  if(p == G4Electron::Electron()) { isElectron = true; }
121  else { isElectron = false;}
122 }
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125 
126 G4double G4eBremParametrizedModel::MinEnergyCut(const G4ParticleDefinition*,
127  const G4MaterialCutsCouple*)
128 {
129  return minThreshold;
130 }
131 
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133 
134 void G4eBremParametrizedModel::SetupForMaterial(const G4ParticleDefinition*,
135  const G4Material* mat,
136  G4double kineticEnergy)
137 {
138  densityFactor = mat->GetElectronDensity()*fMigdalConstant;
139 
140  // calculate threshold for density effect
141  kinEnergy = kineticEnergy;
142  totalEnergy = kineticEnergy + particleMass;
143  densityCorr = densityFactor*totalEnergy*totalEnergy;
144 }
145 
146 
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148 
149 void G4eBremParametrizedModel::Initialise(const G4ParticleDefinition* p,
150  const G4DataVector& cuts)
151 {
152  if(p) { SetParticle(p); }
153 
154  lowKinEnergy = LowEnergyLimit();
155 
156  currentZ = 0.;
157 
158  InitialiseElementSelectors(p, cuts);
159 
160  if(isInitialised) { return; }
161  fParticleChange = GetParticleChangeForLoss();
162  isInitialised = true;
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166 
167 G4double G4eBremParametrizedModel::ComputeDEDXPerVolume(
168  const G4Material* material,
169  const G4ParticleDefinition* p,
170  G4double kineticEnergy,
171  G4double cutEnergy)
172 {
173  if(!particle) { SetParticle(p); }
174  if(kineticEnergy < lowKinEnergy) { return 0.0; }
175  G4double cut = std::min(cutEnergy, kineticEnergy);
176  if(cut == 0.0) { return 0.0; }
177 
178  SetupForMaterial(particle, material,kineticEnergy);
179 
180  const G4ElementVector* theElementVector = material->GetElementVector();
181  const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
182 
183  G4double dedx = 0.0;
184 
185  // loop for elements in the material
186  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
187 
188  G4VEmModel::SetCurrentElement((*theElementVector)[i]);
189  SetCurrentElement((*theElementVector)[i]->GetZ());
190 
191  dedx += theAtomicNumDensityVector[i]*currentZ*currentZ*ComputeBremLoss(cut);
192  }
193  dedx *= bremFactor;
194 
195 
196  return dedx;
197 }
198 
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
200 
201 G4double G4eBremParametrizedModel::ComputeBremLoss(G4double cut)
202 {
203  G4double loss = 0.0;
204 
205  // number of intervals and integration step
206  G4double vcut = cut/totalEnergy;
207  G4int n = (G4int)(20*vcut) + 3;
208  G4double delta = vcut/G4double(n);
209 
210  G4double e0 = 0.0;
211  G4double xs;
212 
213  // integration
214  for(G4int l=0; l<n; l++) {
215 
216  for(G4int i=0; i<8; i++) {
217 
218  G4double eg = (e0 + xgi[i]*delta)*totalEnergy;
219 
220  xs = ComputeDXSectionPerAtom(eg);
221 
222  loss += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
223  }
224  e0 += delta;
225  }
226 
227  loss *= delta*totalEnergy;
228 
229  return loss;
230 }
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 
234 G4double G4eBremParametrizedModel::ComputeCrossSectionPerAtom(
235  const G4ParticleDefinition* p,
236  G4double kineticEnergy,
237  G4double Z, G4double,
238  G4double cutEnergy,
239  G4double maxEnergy)
240 {
241  if(!particle) { SetParticle(p); }
242  if(kineticEnergy < lowKinEnergy) { return 0.0; }
243  G4double cut = std::min(cutEnergy, kineticEnergy);
244  G4double tmax = std::min(maxEnergy, kineticEnergy);
245 
246  if(cut >= tmax) { return 0.0; }
247 
248  SetCurrentElement(Z);
249 
250  G4double cross = ComputeXSectionPerAtom(cut);
251 
252  // allow partial integration
253  if(tmax < kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
254 
255  cross *= Z*Z*bremFactor;
256 
257  return cross;
258 }
259 
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
261 
262 
263 G4double G4eBremParametrizedModel::ComputeXSectionPerAtom(G4double cut)
264 {
265  G4double cross = 0.0;
266 
267  // number of intervals and integration step
268  G4double vcut = log(cut/totalEnergy);
269  G4double vmax = log(kinEnergy/totalEnergy);
270  G4int n = (G4int)(0.45*(vmax - vcut)) + 4;
271  // n=1; // integration test
272  G4double delta = (vmax - vcut)/G4double(n);
273 
274  G4double e0 = vcut;
275  G4double xs;
276 
277  // integration
278  for(G4int l=0; l<n; l++) {
279 
280  for(G4int i=0; i<8; i++) {
281 
282  G4double eg = exp(e0 + xgi[i]*delta)*totalEnergy;
283 
284  xs = ComputeDXSectionPerAtom(eg);
285 
286  cross += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
287  }
288  e0 += delta;
289  }
290 
291  cross *= delta;
292 
293  return cross;
294 }
295 
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
297 
298 //
299 // GEANT4 internal units.
300 //
301  static const G4double
302  ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
303  ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
304  ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
305 
306  static const G4double
307  bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
308  bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
309  bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
310 
311  static const G4double
312  al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
313  al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
314  al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
315 
316  static const G4double
317  bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
318  bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
319  bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
320 
321  static const G4double tlow = 1.*MeV;
322 
323 G4double ScreenFunction1(G4double ScreenVariable)
324 
325 // compute the value of the screening function 3*PHI1 - PHI2
326 
327 {
328  G4double screenVal;
329 
330  if (ScreenVariable > 1.)
331  screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
332  else
333  screenVal = 42.392 - ScreenVariable* (7.796 - 1.961*ScreenVariable);
334 
335  return screenVal;
336 }
337 
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339 
340 G4double ScreenFunction2(G4double ScreenVariable)
341 
342 // compute the value of the screening function 1.5*PHI1 - 0.5*PHI2
343 
344 {
345  G4double screenVal;
346 
347  if (ScreenVariable > 1.)
348  screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
349  else
350  screenVal = 41.734 - ScreenVariable* (6.484 - 1.250*ScreenVariable);
351 
352  return screenVal;
353 }
354 
355 
356 // Parametrized cross section
357 G4double ComputeParametrizedDXSectionPerAtom(G4double kineticEnergy, G4double gammaEnergy, G4double Z) {
358  G4double lnZ = std::log(Z); // 3.*(anElement->GetIonisation()->GetlogZ3());
359  G4double FZ = lnZ* (4.- 0.55*lnZ);
360  G4double ZZ = std::pow (Z*(Z+1.),1./3.); // anElement->GetIonisation()->GetZZ3();
361  G4double Z3 = std::pow (Z,1./3.); // (anElement->GetIonisation()->GetZ3())
362 
363  G4double totalEnergy = kineticEnergy + electron_mass_c2;
364 
365  // G4double x, epsil, greject, migdal, grejmax, q;
366  G4double epsil, greject;
367  G4double U = log(kineticEnergy/electron_mass_c2);
368  G4double U2 = U*U;
369 
370  // precalculated parameters
371  G4double ah, bh;
372 
373 if (kineticEnergy > tlow) {
374 
375  G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12);
376  G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22);
377  G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
378 
379  G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12);
380  G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22);
381  G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
382 
383  ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
384  bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
385 
386  // limit of the screening variable
387  G4double screenfac =
388  136.*electron_mass_c2/(Z3*totalEnergy);
389 
390  epsil = gammaEnergy/totalEnergy; // epsil = x*kineticEnergy/totalEnergy;
391  G4double screenvar = screenfac*epsil/(1.0-epsil);
392  G4double F1 = max(ScreenFunction1(screenvar) - FZ ,0.);
393  G4double F2 = max(ScreenFunction2(screenvar) - FZ ,0.);
394 
395 
396  greject = (F1 - epsil* (ah*F1 - bh*epsil*F2))/8.; // 1./(42.392 - FZ);
397 
398  std::cout << " yy = "<<epsil<<std::endl;
399  std::cout << " F1/(...) "<<F1/(42.392 - FZ)<<std::endl;
400  std::cout << " F2/(...) "<<F2/(42.392 - FZ)<<std::endl;
401  std::cout << " (42.392 - FZ) " << (42.392 - FZ) <<std::endl;
402 
403  } else {
404 
405  G4double al0 = al00 + ZZ* (al01 + ZZ* al02);
406  G4double al1 = al10 + ZZ* (al11 + ZZ* al12);
407  G4double al2 = al20 + ZZ* (al21 + ZZ* al22);
408 
409  G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02);
410  G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12);
411  G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
412 
413  ah = al0 + al1*U + al2*U2;
414  bh = bl0 + bl1*U + bl2*U2;
415 
416  G4double x=gammaEnergy/kineticEnergy;
417  greject=(1. + x* (ah + bh*x));
418 
419  /*
420  // Compute the maximum of the rejection function
421  grejmax = max(1. + xmin* (ah + bh*xmin), 1.+ah+bh);
422  G4double xm = -ah/(2.*bh);
423  if ( xmin < xm && xm < xmax) grejmax = max(grejmax, 1.+ xm* (ah + bh*xm));
424  */
425  }
426 
427  return greject;
428 }
429 
430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
431 
432 G4double G4eBremParametrizedModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
433 {
434 
435  if(gammaEnergy < 0.0) { return 0.0; }
436 
437  G4double y = gammaEnergy/totalEnergy;
438 
439  G4double main=0.;
440  //secondTerm=0.;
441 
442  // ** form factors complete screening case **
443  // only valid for high energies (and if LPM suppression does not play a role)
444  main = (3./4.*y*y - y + 1.) * ( (Fel-fCoulomb) + Finel/currentZ );
445  // secondTerm = (1.-y)/12.*(1.+1./currentZ);
446 
447  std::cout<<" F1(0) "<<ScreenFunction1(0.) <<std::endl;
448  std::cout<<" F1(0) "<<ScreenFunction2(0.) <<std::endl;
449  std::cout<<"Ekin = "<<kinEnergy<<std::endl;
450  std::cout<<"Z = "<<currentZ<<std::endl;
451  std::cout<<"main = "<<main<<std::endl;
452  std::cout<<" y = "<<y<<std::endl;
453  std::cout<<" Fel-fCoulomb "<< (Fel-fCoulomb) <<std::endl;
454 
455  G4double main2 = ComputeParametrizedDXSectionPerAtom(kinEnergy,gammaEnergy,currentZ);
456  std::cout<<"main2 = "<<main2<<std::endl;
457  std::cout<<"main2tot = "<<main2 * ( (Fel-fCoulomb) + Finel/currentZ ) / (Fel-fCoulomb);
458 
459 
460  G4double cross = main2; //main+secondTerm;
461  return cross;
462 }
463 
464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 
466 void G4eBremParametrizedModel::SampleSecondaries(
467  std::vector<G4DynamicParticle*>* vdp,
468  const G4MaterialCutsCouple* couple,
469  const G4DynamicParticle* dp,
470  G4double cutEnergy,
471  G4double maxEnergy)
472 {
473  G4double kineticEnergy = dp->GetKineticEnergy();
474  if(kineticEnergy < lowKinEnergy) { return; }
475  G4double cut = std::min(cutEnergy, kineticEnergy);
476  G4double emax = std::min(maxEnergy, kineticEnergy);
477  if(cut >= emax) { return; }
478 
479  SetupForMaterial(particle, couple->GetMaterial(),kineticEnergy);
480 
481  const G4Element* elm =
482  SelectRandomAtom(couple,particle,kineticEnergy,cut,emax);
483  SetCurrentElement(elm->GetZ());
484 
485  kinEnergy = kineticEnergy;
486  totalEnergy = kineticEnergy + particleMass;
487  densityCorr = densityFactor*totalEnergy*totalEnergy;
488 
489  G4double xmin = log(cut*cut + densityCorr);
490  G4double xmax = log(emax*emax + densityCorr);
491  G4double gammaEnergy, f, x;
492 
493  do {
494  x = exp(xmin + G4UniformRand()*(xmax - xmin)) - densityCorr;
495  if(x < 0.0) x = 0.0;
496  gammaEnergy = sqrt(x);
497  f = ComputeDXSectionPerAtom(gammaEnergy);
498 
499  if ( f > fMax ) {
500  G4cout << "### G4eBremParametrizedModel Warning: Majoranta exceeded! "
501  << f << " > " << fMax
502  << " Egamma(MeV)= " << gammaEnergy
503  << " E(mEV)= " << kineticEnergy
504  << G4endl;
505  }
506 
507  } while (f < fMax*G4UniformRand());
508 
509  //
510  // angles of the emitted gamma. ( Z - axis along the parent particle)
511  // use general interface
512  //
513  G4ThreeVector gammaDirection =
514  GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
515  G4lrint(currentZ),
516  couple->GetMaterial());
517 
518  // create G4DynamicParticle object for the Gamma
519  G4DynamicParticle* gamma = new G4DynamicParticle(theGamma,gammaDirection,
520  gammaEnergy);
521  vdp->push_back(gamma);
522 
523  G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
524  G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
525  - gammaEnergy*gammaDirection).unit();
526 
527  // energy of primary
528  G4double finalE = kineticEnergy - gammaEnergy;
529 
530  // stop tracking and create new secondary instead of primary
531  if(gammaEnergy > SecondaryThreshold()) {
532  fParticleChange->ProposeTrackStatus(fStopAndKill);
533  fParticleChange->SetProposedKineticEnergy(0.0);
534  G4DynamicParticle* el =
535  new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
536  direction, finalE);
537  vdp->push_back(el);
538 
539  // continue tracking
540  } else {
541  fParticleChange->SetProposedMomentumDirection(direction);
542  fParticleChange->SetProposedKineticEnergy(finalE);
543  }
544 }
545 
546 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
547 
548 
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4