Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeBremsstrahlungModel.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 // Author: Luciano Pandola
29 //
30 // History:
31 // --------
32 // 23 Nov 2010 L Pandola First complete implementation, Penelope v2008
33 // 24 May 2011 L. Pandola Renamed (make default Penelope)
34 // 13 Mar 2012 L. Pandola Updated the interface for the angular generator
35 // 18 Jul 2012 L. Pandola Migrate to the new interface of the angular generator, which
36 // now provides the G4ThreeVector and takes care of rotation
37 //
38 
39 #include "G4PenelopeBremsstrahlungModel.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4PenelopeBremsstrahlungFS.hh"
43 #include "G4PenelopeBremsstrahlungAngular.hh"
44 #include "G4ParticleDefinition.hh"
45 #include "G4MaterialCutsCouple.hh"
46 #include "G4ProductionCutsTable.hh"
47 #include "G4DynamicParticle.hh"
48 #include "G4Gamma.hh"
49 #include "G4Electron.hh"
50 #include "G4Positron.hh"
51 #include "G4PenelopeOscillatorManager.hh"
52 #include "G4PenelopeCrossSection.hh"
53 #include "G4PhysicsFreeVector.hh"
54 #include "G4PhysicsLogVector.hh"
55 #include "G4PhysicsTable.hh"
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58 
59 G4PenelopeBremsstrahlungModel::G4PenelopeBremsstrahlungModel(const G4ParticleDefinition*,
60  const G4String& nam)
61  :G4VEmModel(nam),fParticleChange(0),isInitialised(false),energyGrid(0),
62  XSTableElectron(0),XSTablePositron(0)
63 
64 {
65  fIntrinsicLowEnergyLimit = 100.0*eV;
66  fIntrinsicHighEnergyLimit = 100.0*GeV;
67  nBins = 200;
68 
69  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
70  //
71  oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
72  //
73  verboseLevel= 0;
74 
75  // Verbosity scale:
76  // 0 = nothing
77  // 1 = warning for energy non-conservation
78  // 2 = details of energy budget
79  // 3 = calculation of cross sections, file openings, sampling of atoms
80  // 4 = entering in methods
81 
82  fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS();
83  fPenelopeAngular = new G4PenelopeBremsstrahlungAngular();
84 
85  // Atomic deexcitation model activated by default
86  SetDeexcitationFlag(true);
87 
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
92 G4PenelopeBremsstrahlungModel::~G4PenelopeBremsstrahlungModel()
93 {
94  ClearTables();
95  if (fPenelopeFSHelper)
96  delete fPenelopeFSHelper;
97  if (fPenelopeAngular)
98  delete fPenelopeAngular;
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 
103 void G4PenelopeBremsstrahlungModel::Initialise(const G4ParticleDefinition*,
104  const G4DataVector&)
105 {
106  if (verboseLevel > 3)
107  G4cout << "Calling G4PenelopeBremsstrahlungModel::Initialise()" << G4endl;
108 
109  //Clear and re-build the tables
110  ClearTables();
111 
112  //forces the cleaning of tables, in this specific case
113  if (fPenelopeAngular)
114  fPenelopeAngular->Initialize();
115 
116 
117  //Set the number of bins for the tables. 20 points per decade
118  nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
119  nBins = std::max(nBins,(size_t)100);
120  energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
121  HighEnergyLimit(),
122  nBins-1); //one hidden bin is added
123 
124 
125  XSTableElectron = new
126  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
127  XSTablePositron = new
128  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
129 
130  if (verboseLevel > 2) {
131  G4cout << "Penelope Bremsstrahlung model v2008 is initialized " << G4endl
132  << "Energy range: "
133  << LowEnergyLimit() / keV << " keV - "
134  << HighEnergyLimit() / GeV << " GeV."
135  << G4endl;
136  }
137 
138  if(isInitialised) return;
139  fParticleChange = GetParticleChangeForLoss();
140  isInitialised = true;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144 
145 G4double G4PenelopeBremsstrahlungModel::CrossSectionPerVolume(const G4Material* material,
146  const G4ParticleDefinition* theParticle,
147  G4double energy,
148  G4double cutEnergy,
149  G4double)
150 {
151  //
152  if (verboseLevel > 3)
153  G4cout << "Calling CrossSectionPerVolume() of G4PenelopeBremsstrahlungModel" << G4endl;
154 
155  SetupForMaterial(theParticle, material, energy);
156 
157  G4double crossPerMolecule = 0.;
158 
159  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
160  cutEnergy);
161 
162  if (theXS)
163  crossPerMolecule = theXS->GetHardCrossSection(energy);
164 
165  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
166  G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
167 
168  if (verboseLevel > 3)
169  G4cout << "Material " << material->GetName() << " has " << atPerMol <<
170  "atoms per molecule" << G4endl;
171 
172  G4double moleculeDensity = 0.;
173  if (atPerMol)
174  moleculeDensity = atomDensity/atPerMol;
175 
176  G4double crossPerVolume = crossPerMolecule*moleculeDensity;
177 
178  if (verboseLevel > 2)
179  {
180  G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
181  G4cout << "Mean free path for gamma emission > " << cutEnergy/keV << " keV at " <<
182  energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
183  }
184 
185  return crossPerVolume;
186 }
187 
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189 
190 //This is a dummy method. Never inkoved by the tracking, it just issues
191 //a warning if one tries to get Cross Sections per Atom via the
192 //G4EmCalculator.
193 G4double G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
194  G4double,
195  G4double,
196  G4double,
197  G4double,
198  G4double)
199 {
200  G4cout << "*** G4PenelopeBremsstrahlungModel -- WARNING ***" << G4endl;
201  G4cout << "Penelope Bremsstrahlung model v2008 does not calculate cross section _per atom_ " << G4endl;
202  G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
203  G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
204  return 0;
205 }
206 
207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
208 
209 G4double G4PenelopeBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* material,
210  const G4ParticleDefinition* theParticle,
211  G4double kineticEnergy,
212  G4double cutEnergy)
213 {
214  if (verboseLevel > 3)
215  G4cout << "Calling ComputeDEDX() of G4PenelopeBremsstrahlungModel" << G4endl;
216 
217  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
218  cutEnergy);
219  G4double sPowerPerMolecule = 0.0;
220  if (theXS)
221  sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
222 
223 
224  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
225  G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
226 
227  G4double moleculeDensity = 0.;
228  if (atPerMol)
229  moleculeDensity = atomDensity/atPerMol;
230 
231  G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
232 
233  if (verboseLevel > 2)
234  {
235  G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
236  G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
237  kineticEnergy/keV << " keV = " <<
238  sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
239  }
240  return sPowerPerVolume;
241 }
242 
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244 
245 void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>*fvect,
246  const G4MaterialCutsCouple* couple,
247  const G4DynamicParticle* aDynamicParticle,
248  G4double cutG,
249  G4double)
250 {
251  if (verboseLevel > 3)
252  G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
253 
254  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
255  const G4Material* material = couple->GetMaterial();
256 
257  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
258  {
259  fParticleChange->SetProposedKineticEnergy(0.);
260  fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
261  return ;
262  }
263 
264  G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
265  //This is the momentum
266  G4ThreeVector initialMomentum = aDynamicParticle->GetMomentum();
267 
268  //Not enough energy to produce a secondary! Return with nothing happened
269  if (kineticEnergy < cutG)
270  return;
271 
272  if (verboseLevel > 3)
273  G4cout << "Going to sample gamma energy for: " <<material->GetName() << " " <<
274  "energy = " << kineticEnergy/keV << ", cut = " << cutG/keV << G4endl;
275 
276  //Sample gamma's energy according to the spectrum
277  G4double gammaEnergy =
278  fPenelopeFSHelper->SampleGammaEnergy(kineticEnergy,material,cutG);
279 
280  if (verboseLevel > 3)
281  G4cout << "Sampled gamma energy: " << gammaEnergy/keV << " keV" << G4endl;
282 
283  //Now sample the direction for the Gamma. Notice that the rotation is already done
284  //Z is unused here, I plug 0. The information is in the material pointer
285  G4ThreeVector gammaDirection1 =
286  fPenelopeAngular->SampleDirection(aDynamicParticle,gammaEnergy,0,material);
287 
288  if (verboseLevel > 3)
289  G4cout << "Sampled cosTheta for e-: " << gammaDirection1.cosTheta() << G4endl;
290 
291  G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
292  if (residualPrimaryEnergy < 0)
293  {
294  //Ok we have a problem, all energy goes with the gamma
295  gammaEnergy += residualPrimaryEnergy;
296  residualPrimaryEnergy = 0.0;
297  }
298 
299  //Produce final state according to momentum conservation
300  G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
301  particleDirection1 = particleDirection1.unit(); //normalize
302 
303  //Update the primary particle
304  if (residualPrimaryEnergy > 0.)
305  {
306  fParticleChange->ProposeMomentumDirection(particleDirection1);
307  fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
308  }
309  else
310  {
311  fParticleChange->SetProposedKineticEnergy(0.);
312  }
313 
314  //Now produce the photon
315  G4DynamicParticle* theGamma = new G4DynamicParticle(G4Gamma::Gamma(),
316  gammaDirection1,
317  gammaEnergy);
318  fvect->push_back(theGamma);
319 
320  if (verboseLevel > 1)
321  {
322  G4cout << "-----------------------------------------------------------" << G4endl;
323  G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
324  G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
325  G4cout << "-----------------------------------------------------------" << G4endl;
326  G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
327  G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
328  G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV
329  << " keV" << G4endl;
330  G4cout << "-----------------------------------------------------------" << G4endl;
331  }
332 
333  if (verboseLevel > 0)
334  {
335  G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
336  if (energyDiff > 0.05*keV)
337  G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: "
338  <<
339  (residualPrimaryEnergy+gammaEnergy)/keV <<
340  " keV (final) vs. " <<
341  kineticEnergy/keV << " keV (initial)" << G4endl;
342  }
343  return;
344 }
345 
346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
347 
348 void G4PenelopeBremsstrahlungModel::ClearTables()
349 {
350  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>::iterator i;
351  if (XSTableElectron)
352  {
353  for (i=XSTableElectron->begin(); i != XSTableElectron->end(); i++)
354  {
355  G4PenelopeCrossSection* tab = i->second;
356  delete tab;
357  }
358  delete XSTableElectron;
359  XSTableElectron = 0;
360  }
361  if (XSTablePositron)
362  {
363  for (i=XSTablePositron->begin(); i != XSTablePositron->end(); i++)
364  {
365  G4PenelopeCrossSection* tab = i->second;
366  delete tab;
367  }
368  delete XSTablePositron;
369  XSTablePositron = 0;
370  }
371 
372  if (energyGrid)
373  delete energyGrid;
374 
375  if (fPenelopeFSHelper)
376  fPenelopeFSHelper->ClearTables();
377 
378  if (verboseLevel > 2)
379  G4cout << "G4PenelopeBremsstrahlungModel: cleared tables" << G4endl;
380 
381  return;
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
385 
386 G4double G4PenelopeBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
387  const G4MaterialCutsCouple*)
388 {
389  return fIntrinsicLowEnergyLimit;
390 }
391 
392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
393 
394 void G4PenelopeBremsstrahlungModel::BuildXSTable(const G4Material* mat,G4double cut)
395 {
396 //
397  //This method fills the G4PenelopeCrossSection containers for electrons or positrons
398  //and for the given material/cut couple.
399  //Equivalent of subroutines EBRaT and PINaT of Penelope
400  //
401  if (verboseLevel > 2)
402  {
403  G4cout << "G4PenelopeBremsstrahlungModel: going to build cross section table " << G4endl;
404  G4cout << "for e+/e- in " << mat->GetName() << G4endl;
405  }
406 
407  //Tables have been already created (checked by GetCrossSectionTableForCouple)
408  if (energyGrid->GetVectorLength() != nBins)
409  {
410  G4ExceptionDescription ed;
411  ed << "Energy Grid looks not initialized" << G4endl;
412  ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
413  G4Exception("G4PenelopeBremsstrahlungModel::BuildXSTable()",
414  "em2016",FatalException,ed);
415  }
416 
417  G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins);
418  G4PenelopeCrossSection* XSEntry2 = new G4PenelopeCrossSection(nBins);
419 
420  G4PhysicsTable* table = fPenelopeFSHelper->GetScaledXSTable(mat,cut);
421 
422 
423  //loop on the energy grid
424  for (size_t bin=0;bin<nBins;bin++)
425  {
426  G4double energy = energyGrid->GetLowEdgeEnergy(bin);
427  G4double XH0=0, XH1=0, XH2=0;
428  G4double XS0=0, XS1=0, XS2=0;
429 
430  //Global xs factor
431  G4double fact = fPenelopeFSHelper->GetEffectiveZSquared(mat)*
432  ((energy+electron_mass_c2)*(energy+electron_mass_c2)/
433  (energy*(energy+2.0*electron_mass_c2)));
434 
435  G4double restrictedCut = cut/energy;
436 
437  //Now I need the dSigma/dX profile - interpolated on energy - for
438  //the 32-point x grid. Interpolation is log-log
439  size_t nBinsX = fPenelopeFSHelper->GetNBinsX();
440  G4double* tempData = new G4double[nBinsX];
441  G4double logene = std::log(energy);
442  for (size_t ix=0;ix<nBinsX;ix++)
443  {
444  //find dSigma/dx for the given E. X belongs to the 32-point grid.
445  G4double val = (*table)[ix]->Value(logene);
446  tempData[ix] = std::exp(val); //back to the real value!
447  }
448 
449  G4double XH0A = 0.;
450  if (restrictedCut <= 1) //calculate only if we are above threshold!
451  XH0A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,-1) -
452  fPenelopeFSHelper->GetMomentumIntegral(tempData,restrictedCut,-1);
453  G4double XS1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
454  restrictedCut,0);
455  G4double XS2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
456  restrictedCut,1);
457  G4double XH1A=0, XH2A=0;
458  if (restrictedCut <=1)
459  {
460  XH1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,0) -
461  XS1A;
462  XH2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,1) -
463  XS2A;
464  }
465  delete[] tempData;
466 
467  XH0 = XH0A*fact;
468  XS1 = XS1A*fact*energy;
469  XH1 = XH1A*fact*energy;
470  XS2 = XS2A*fact*energy*energy;
471  XH2 = XH2A*fact*energy*energy;
472 
473  XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
474 
475  //take care of positrons
476  G4double posCorrection = GetPositronXSCorrection(mat,energy);
477  XSEntry2->AddCrossSectionPoint(bin,energy,XH0*posCorrection,
478  XH1*posCorrection,
479  XH2*posCorrection,
480  XS0,
481  XS1*posCorrection,
482  XS2*posCorrection);
483  }
484 
485  //Insert in the appropriate table
486  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
487  XSTableElectron->insert(std::make_pair(theKey,XSEntry));
488  XSTablePositron->insert(std::make_pair(theKey,XSEntry2));
489 
490  return;
491 }
492 
493 
494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
495 
496 G4PenelopeCrossSection*
497 G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
498  const G4Material* mat,
499  G4double cut)
500 {
501  if (part != G4Electron::Electron() && part != G4Positron::Positron())
502  {
503  G4ExceptionDescription ed;
504  ed << "Invalid particle: " << part->GetParticleName() << G4endl;
505  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
506  "em0001",FatalException,ed);
507  return NULL;
508  }
509 
510  if (part == G4Electron::Electron())
511  {
512  if (!XSTableElectron)
513  {
514  G4String excep = "The Cross Section Table for e- was not initialized correctly!";
515  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
516  "em2013",FatalException,excep);
517  return NULL;
518  }
519  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
520  if (XSTableElectron->count(theKey)) //table already built
521  return XSTableElectron->find(theKey)->second;
522  else
523  {
524  BuildXSTable(mat,cut);
525  if (XSTableElectron->count(theKey)) //now it should be ok!
526  return XSTableElectron->find(theKey)->second;
527  else
528  {
529  G4ExceptionDescription ed;
530  ed << "Unable to build e- table for " << mat->GetName() << G4endl;
531  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
532  "em2009",FatalException,ed);
533  }
534  }
535  }
536  if (part == G4Positron::Positron())
537  {
538  if (!XSTablePositron)
539  {
540  G4String excep = "The Cross Section Table for e+ was not initialized correctly!";
541  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
542  "em2013",FatalException,excep);
543  return NULL;
544  }
545  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
546  if (XSTablePositron->count(theKey)) //table already built
547  return XSTablePositron->find(theKey)->second;
548  else
549  {
550  BuildXSTable(mat,cut);
551  if (XSTablePositron->count(theKey)) //now it should be ok!
552  return XSTablePositron->find(theKey)->second;
553  else
554  {
555  G4ExceptionDescription ed;
556  ed << "Unable to build e+ table for " << mat->GetName() << G4endl;
557  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
558  "em2009",FatalException,ed);
559  }
560  }
561  }
562  return NULL;
563 }
564 
565 
566 
567 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
568 
569 G4double G4PenelopeBremsstrahlungModel::GetPositronXSCorrection(const G4Material* mat,
570  G4double energy)
571 {
572  //The electron-to-positron correction factor is set equal to the ratio of the
573  //radiative stopping powers for positrons and electrons, which has been calculated
574  //by Kim et al. (1986) (cf. Berger and Seltzer, 1982). Here, it is used an
575  //analytical approximation which reproduces the tabulated values with 0.5%
576  //accuracy
577  G4double t=std::log(1.0+1e6*energy/
578  (electron_mass_c2*fPenelopeFSHelper->GetEffectiveZSquared(mat)));
579  G4double corr = 1.0-std::exp(-t*(1.2359e-1-t*(6.1274e-2-t*
580  (3.1516e-2-t*(7.7446e-3-t*(1.0595e-3-t*
581  (7.0568e-5-t*
582  1.8080e-6)))))));
583  return corr;
584 }
Generated on Sat May 25 2013 14:33:34 for Geant4 by   doxygen 1.8.4