Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeGammaConversionModel.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 // 13 Jan 2010 L Pandola First implementation (updated to Penelope08)
33 // 24 May 2011 L Pandola Renamed (make v2008 as default Penelope)
34 //
35 
36 #include "G4PenelopeGammaConversionModel.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4ParticleDefinition.hh"
40 #include "G4MaterialCutsCouple.hh"
41 #include "G4ProductionCutsTable.hh"
42 #include "G4DynamicParticle.hh"
43 #include "G4Element.hh"
44 #include "G4Gamma.hh"
45 #include "G4Electron.hh"
46 #include "G4Positron.hh"
47 #include "G4PhysicsFreeVector.hh"
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 
51 
52 G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel(const G4ParticleDefinition*,
53  const G4String& nam)
54  :G4VEmModel(nam),fParticleChange(0),logAtomicCrossSection(0),
55  fEffectiveCharge(0),fMaterialInvScreeningRadius(0),
56  fScreeningFunction(0),isInitialised(false)
57 {
58  fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
59  fIntrinsicHighEnergyLimit = 100.0*GeV;
60  fSmallEnergy = 1.1*MeV;
61  InitializeScreeningRadii();
62 
63  // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
64  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
65  //
66  verboseLevel= 0;
67  // Verbosity scale:
68  // 0 = nothing
69  // 1 = warning for energy non-conservation
70  // 2 = details of energy budget
71  // 3 = calculation of cross sections, file openings, sampling of atoms
72  // 4 = entering in methods
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
77 G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel()
78 {
79  std::map <G4int,G4PhysicsFreeVector*>::iterator i;
80  if (logAtomicCrossSection)
81  {
82  for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++)
83  if (i->second) delete i->second;
84  delete logAtomicCrossSection;
85  }
86  if (fEffectiveCharge)
87  delete fEffectiveCharge;
88  if (fMaterialInvScreeningRadius)
89  delete fMaterialInvScreeningRadius;
90  if (fScreeningFunction)
91  delete fScreeningFunction;
92 }
93 
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 
97 void G4PenelopeGammaConversionModel::Initialise(const G4ParticleDefinition*,
98  const G4DataVector&)
99 {
100  if (verboseLevel > 3)
101  G4cout << "Calling G4PenelopeGammaConversionModel::Initialise()" << G4endl;
102 
103  // logAtomicCrossSection is created only once, since it is never cleared
104  if (!logAtomicCrossSection)
105  logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
106 
107  //delete old material data...
108  if (fEffectiveCharge)
109  {
110  delete fEffectiveCharge;
111  fEffectiveCharge = 0;
112  }
113  if (fMaterialInvScreeningRadius)
114  {
115  delete fMaterialInvScreeningRadius;
116  fMaterialInvScreeningRadius = 0;
117  }
118  if (fScreeningFunction)
119  {
120  delete fScreeningFunction;
121  fScreeningFunction = 0;
122  }
123  //and create new ones
124  fEffectiveCharge = new std::map<const G4Material*,G4double>;
125  fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
126  fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
127 
128  if (verboseLevel > 0) {
129  G4cout << "Penelope Gamma Conversion model v2008 is initialized " << G4endl
130  << "Energy range: "
131  << LowEnergyLimit() / MeV << " MeV - "
132  << HighEnergyLimit() / GeV << " GeV"
133  << G4endl;
134  }
135 
136  if(isInitialised) return;
137  fParticleChange = GetParticleChangeForGamma();
138  isInitialised = true;
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 
143 G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom(
144  const G4ParticleDefinition*,
145  G4double energy,
146  G4double Z, G4double,
147  G4double, G4double)
148 {
149  //
150  // Penelope model v2008.
151  // Cross section (including triplet production) read from database and managed
152  // through the G4CrossSectionHandler utility. Cross section data are from
153  // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
154  //
155 
156  if (energy < fIntrinsicLowEnergyLimit)
157  return 0;
158 
159  G4int iZ = (G4int) Z;
160 
161  //read data files
162  if (!logAtomicCrossSection->count(iZ))
163  ReadDataFile(iZ);
164  //now it should be ok
165  if (!logAtomicCrossSection->count(iZ))
166  {
167  G4ExceptionDescription ed;
168  ed << "Unable to retrieve cross section table for Z=" << iZ << G4endl;
169  G4Exception("G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom()",
170  "em2018",FatalException,ed);
171  }
172 
173  G4double cs = 0;
174  G4double logene = std::log(energy);
175  G4PhysicsFreeVector* theVec = logAtomicCrossSection->find(iZ)->second;
176 
177  G4double logXS = theVec->Value(logene);
178  cs = std::exp(logXS);
179 
180  if (verboseLevel > 2)
181  G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z <<
182  " = " << cs/barn << " barn" << G4endl;
183  return cs;
184 }
185 
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187 
188 void
189 G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
190  const G4MaterialCutsCouple* couple,
191  const G4DynamicParticle* aDynamicGamma,
192  G4double,
193  G4double)
194 {
195  //
196  // Penelope model v2008.
197  // Final state is sampled according to the Bethe-Heitler model with Coulomb
198  // corrections, according to the semi-empirical model of
199  // J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531.
200  //
201  // The model uses the high energy Coulomb correction from
202  // H. Davies et al., Phys. Rev. 93 (1954) 788
203  // and atomic screening radii tabulated from
204  // J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023
205  // for Z= 1 to 92.
206  //
207  if (verboseLevel > 3)
208  G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl;
209 
210  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
211 
212  // Always kill primary
213  fParticleChange->ProposeTrackStatus(fStopAndKill);
214  fParticleChange->SetProposedKineticEnergy(0.);
215 
216  if (photonEnergy <= fIntrinsicLowEnergyLimit)
217  {
218  fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
219  return ;
220  }
221 
222  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
223  const G4Material* mat = couple->GetMaterial();
224 
225  //check if material data are available
226  if (!fEffectiveCharge->count(mat))
227  InitializeScreeningFunctions(mat);
228  if (!fEffectiveCharge->count(mat))
229  {
230  G4ExceptionDescription ed;
231  ed << "Unable to allocate the EffectiveCharge data for " <<
232  mat->GetName() << G4endl;
233  G4Exception("G4PenelopeGammaConversion::SampleSecondaries()",
234  "em2019",FatalException,ed);
235  }
236 
237  // eps is the fraction of the photon energy assigned to e- (including rest mass)
238  G4double eps = 0;
239  G4double eki = electron_mass_c2/photonEnergy;
240 
241  //Do it fast for photon energy < 1.1 MeV (close to threshold)
242  if (photonEnergy < fSmallEnergy)
243  eps = eki + (1.0-2.0*eki)*G4UniformRand();
244  else
245  {
246  //Complete calculation
247  G4double effC = fEffectiveCharge->find(mat)->second;
248  G4double alz = effC*fine_structure_const;
249  G4double T = std::sqrt(2.0*eki);
250  G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T
251  +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T
252  -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T
253  +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T;
254 
255  G4double F0b = fScreeningFunction->find(mat)->second.second;
256  G4double g0 = F0b + F00;
257  G4double invRad = fMaterialInvScreeningRadius->find(mat)->second;
258  G4double bmin = 4.0*eki/invRad;
259  std::pair<G4double,G4double> scree = GetScreeningFunctions(bmin);
260  G4double g1 = scree.first;
261  G4double g2 = scree.second;
262  G4double g1min = g1+g0;
263  G4double g2min = g2+g0;
264  G4double xr = 0.5-eki;
265  G4double a1 = 2.*g1min*xr*xr/3.;
266  G4double p1 = a1/(a1+g2min);
267 
268  G4bool loopAgain = false;
269  //Random sampling of eps
270  do{
271  loopAgain = false;
272  if (G4UniformRand() <= p1)
273  {
274  G4double ru2m1 = 2.0*G4UniformRand()-1.0;
275  if (ru2m1 < 0)
276  eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
277  else
278  eps = 0.5+xr*std::pow(ru2m1,1./3.);
279  G4double B = eki/(invRad*eps*(1.0-eps));
280  scree = GetScreeningFunctions(B);
281  g1 = scree.first;
282  g1 = std::max(g1+g0,0.);
283  if (G4UniformRand()*g1min > g1)
284  loopAgain = true;
285  }
286  else
287  {
288  eps = eki+2.0*xr*G4UniformRand();
289  G4double B = eki/(invRad*eps*(1.0-eps));
290  scree = GetScreeningFunctions(B);
291  g2 = scree.second;
292  g2 = std::max(g2+g0,0.);
293  if (G4UniformRand()*g2min > g2)
294  loopAgain = true;
295  }
296  }while(loopAgain);
297 
298  }
299  if (verboseLevel > 4)
300  G4cout << "Sampled eps = " << eps << G4endl;
301 
302  G4double electronTotEnergy = eps*photonEnergy;
303  G4double positronTotEnergy = (1.0-eps)*photonEnergy;
304 
305  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
306 
307  //electron kinematics
308  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
309  G4double costheta_el = G4UniformRand()*2.0-1.0;
310  G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
311  costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
312  G4double phi_el = twopi * G4UniformRand() ;
313  G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
314  G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
315  G4double dirZ_el = costheta_el;
316 
317  //positron kinematics
318  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
319  G4double costheta_po = G4UniformRand()*2.0-1.0;
320  kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
321  costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
322  G4double phi_po = twopi * G4UniformRand() ;
323  G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
324  G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
325  G4double dirZ_po = costheta_po;
326 
327  // Kinematics of the created pair:
328  // the electron and positron are assumed to have a symetric angular
329  // distribution with respect to the Z axis along the parent photon
330  G4double localEnergyDeposit = 0. ;
331 
332  if (electronKineEnergy > 0.0)
333  {
334  G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
335  electronDirection.rotateUz(photonDirection);
336  G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
337  electronDirection,
338  electronKineEnergy);
339  fvect->push_back(electron);
340  }
341  else
342  {
343  localEnergyDeposit += electronKineEnergy;
344  electronKineEnergy = 0;
345  }
346 
347  //Generate the positron. Real particle in any case, because it will annihilate. If below
348  //threshold, produce it at rest
349  if (positronKineEnergy < 0.0)
350  {
351  localEnergyDeposit += positronKineEnergy;
352  positronKineEnergy = 0; //produce it at rest
353  }
354  G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
355  positronDirection.rotateUz(photonDirection);
356  G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(),
357  positronDirection, positronKineEnergy);
358  fvect->push_back(positron);
359 
360  //Add rest of energy to the local energy deposit
361  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
362 
363  if (verboseLevel > 1)
364  {
365  G4cout << "-----------------------------------------------------------" << G4endl;
366  G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl;
367  G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
368  G4cout << "-----------------------------------------------------------" << G4endl;
369  if (electronKineEnergy)
370  G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV"
371  << G4endl;
372  if (positronKineEnergy)
373  G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
374  G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl;
375  if (localEnergyDeposit)
376  G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
377  G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+
378  localEnergyDeposit+2.0*electron_mass_c2)/keV <<
379  " keV" << G4endl;
380  G4cout << "-----------------------------------------------------------" << G4endl;
381  }
382  if (verboseLevel > 0)
383  {
384  G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
385  localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
386  if (energyDiff > 0.05*keV)
387  G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: "
388  << (electronKineEnergy+positronKineEnergy+
389  localEnergyDeposit+2.0*electron_mass_c2)/keV
390  << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl;
391  }
392 }
393 
394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
395 
396 void G4PenelopeGammaConversionModel::ReadDataFile(const G4int Z)
397 {
398  if (verboseLevel > 2)
399  {
400  G4cout << "G4PenelopeGammaConversionModel::ReadDataFile()" << G4endl;
401  G4cout << "Going to read Gamma Conversion data files for Z=" << Z << G4endl;
402  }
403 
404  char* path = getenv("G4LEDATA");
405  if (!path)
406  {
407  G4String excep =
408  "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
409  G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
410  "em0006",FatalException,excep);
411  return;
412  }
413 
414  /*
415  Read the cross section file
416  */
417  std::ostringstream ost;
418  if (Z>9)
419  ost << path << "/penelope/pairproduction/pdgpp" << Z << ".p08";
420  else
421  ost << path << "/penelope/pairproduction/pdgpp0" << Z << ".p08";
422  std::ifstream file(ost.str().c_str());
423  if (!file.is_open())
424  {
425  G4String excep = "G4PenelopeGammaConversionModel - data file " +
426  G4String(ost.str()) + " not found!";
427  G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
428  "em0003",FatalException,excep);
429  }
430 
431  //I have to know in advance how many points are in the data list
432  //to initialize the G4PhysicsFreeVector()
433  size_t ndata=0;
434  G4String line;
435  while( getline(file, line) )
436  ndata++;
437  ndata -= 1; //remove one header line
438  //G4cout << "Found: " << ndata << " lines" << G4endl;
439 
440  file.clear();
441  file.close();
442  file.open(ost.str().c_str());
443  G4int readZ =0;
444  file >> readZ;
445 
446  if (verboseLevel > 3)
447  G4cout << "Element Z=" << Z << G4endl;
448 
449  //check the right file is opened.
450  if (readZ != Z)
451  {
452  G4ExceptionDescription ed;
453  ed << "Corrupted data file for Z=" << Z << G4endl;
454  G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
455  "em0005",FatalException,ed);
456  }
457 
458  G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(ndata);
459  G4double ene=0,xs=0;
460  for (size_t i=0;i<ndata;i++)
461  {
462  file >> ene >> xs;
463  //dimensional quantities
464  ene *= eV;
465  xs *= barn;
466  if (xs < 1e-40*cm2) //protection against log(0)
467  xs = 1e-40*cm2;
468  theVec->PutValue(i,std::log(ene),std::log(xs));
469  }
470  file.close();
471 
472  if (!logAtomicCrossSection)
473  {
474  G4ExceptionDescription ed;
475  ed << "Problem with allocation of logAtomicCrossSection data table " << G4endl;
476  G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
477  "em2020",FatalException,ed);
478  delete theVec;
479  return;
480  }
481  logAtomicCrossSection->insert(std::make_pair(Z,theVec));
482 
483  return;
484 
485 }
486 
487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
488 
489 void G4PenelopeGammaConversionModel::InitializeScreeningRadii()
490 {
491  G4double temp[99] = {1.2281e+02,7.3167e+01,6.9228e+01,6.7301e+01,6.4696e+01,
492  6.1228e+01,5.7524e+01,5.4033e+01,5.0787e+01,4.7851e+01,4.6373e+01,
493  4.5401e+01,4.4503e+01,4.3815e+01,4.3074e+01,4.2321e+01,4.1586e+01,
494  4.0953e+01,4.0524e+01,4.0256e+01,3.9756e+01,3.9144e+01,3.8462e+01,
495  3.7778e+01,3.7174e+01,3.6663e+01,3.5986e+01,3.5317e+01,3.4688e+01,
496  3.4197e+01,3.3786e+01,3.3422e+01,3.3068e+01,3.2740e+01,3.2438e+01,
497  3.2143e+01,3.1884e+01,3.1622e+01,3.1438e+01,3.1142e+01,3.0950e+01,
498  3.0758e+01,3.0561e+01,3.0285e+01,3.0097e+01,2.9832e+01,2.9581e+01,
499  2.9411e+01,2.9247e+01,2.9085e+01,2.8930e+01,2.8721e+01,2.8580e+01,
500  2.8442e+01,2.8312e+01,2.8139e+01,2.7973e+01,2.7819e+01,2.7675e+01,
501  2.7496e+01,2.7285e+01,2.7093e+01,2.6911e+01,2.6705e+01,2.6516e+01,
502  2.6304e+01,2.6108e+01,2.5929e+01,2.5730e+01,2.5577e+01,2.5403e+01,
503  2.5245e+01,2.5100e+01,2.4941e+01,2.4790e+01,2.4655e+01,2.4506e+01,
504  2.4391e+01,2.4262e+01,2.4145e+01,2.4039e+01,2.3922e+01,2.3813e+01,
505  2.3712e+01,2.3621e+01,2.3523e+01,2.3430e+01,2.3331e+01,2.3238e+01,
506  2.3139e+01,2.3048e+01,2.2967e+01,2.2833e+01,2.2694e+01,2.2624e+01,
507  2.2545e+01,2.2446e+01,2.2358e+01,2.2264e+01};
508 
509  //copy temporary vector in class data member
510  for (G4int i=0;i<99;i++)
511  fAtomicScreeningRadius[i] = temp[i];
512 }
513 
514 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
515 
516 void G4PenelopeGammaConversionModel::InitializeScreeningFunctions(const G4Material* material)
517 {
518  // This is subroutine GPPa0 of Penelope
519  //
520  // 1) calculate the effective Z for the purpose
521  //
522  G4double zeff = 0;
523  G4int intZ = 0;
524  G4int nElements = material->GetNumberOfElements();
525  const G4ElementVector* elementVector = material->GetElementVector();
526 
527  //avoid calculations if only one building element!
528  if (nElements == 1)
529  {
530  zeff = (*elementVector)[0]->GetZ();
531  intZ = (G4int) zeff;
532  }
533  else // many elements...let's do the calculation
534  {
535  const G4double* fractionVector = material->GetVecNbOfAtomsPerVolume();
536 
537  G4double atot = 0;
538  for (G4int i=0;i<nElements;i++)
539  {
540  G4double Zelement = (*elementVector)[i]->GetZ();
541  G4double Aelement = (*elementVector)[i]->GetAtomicMassAmu();
542  atot += Aelement*fractionVector[i];
543  zeff += Zelement*Aelement*fractionVector[i]; //average with the number of nuclei
544  }
545  atot /= material->GetTotNbOfAtomsPerVolume();
546  zeff /= (material->GetTotNbOfAtomsPerVolume()*atot);
547 
548  intZ = (G4int) (zeff+0.25);
549  if (intZ <= 0)
550  intZ = 1;
551  if (intZ > 99)
552  intZ = 99;
553  }
554 
555  if (fEffectiveCharge)
556  fEffectiveCharge->insert(std::make_pair(material,zeff));
557 
558  //
559  // 2) Calculate Coulomb Correction
560  //
561  G4double alz = fine_structure_const*zeff;
562  G4double alzSquared = alz*alz;
563  G4double fc = alzSquared*(0.202059-alzSquared*
564  (0.03693-alzSquared*
565  (0.00835-alzSquared*(0.00201-alzSquared*
566  (0.00049-alzSquared*
567  (0.00012-alzSquared*0.00003)))))
568  +1.0/(alzSquared+1.0));
569  //
570  // 3) Screening functions and low-energy corrections
571  //
572  G4double matRadius = 2.0/ fAtomicScreeningRadius[intZ-1];
573  if (fMaterialInvScreeningRadius)
574  fMaterialInvScreeningRadius->insert(std::make_pair(material,matRadius));
575 
576  std::pair<G4double,G4double> myPair(0,0);
577  G4double f0a = 4.0*std::log(fAtomicScreeningRadius[intZ-1]);
578  G4double f0b = f0a - 4.0*fc;
579  myPair.first = f0a;
580  myPair.second = f0b;
581 
582  if (fScreeningFunction)
583  fScreeningFunction->insert(std::make_pair(material,myPair));
584 
585  if (verboseLevel > 2)
586  {
587  G4cout << "Average Z for material " << material->GetName() << " = " <<
588  zeff << G4endl;
589  G4cout << "Effective radius for material " << material->GetName() << " = " <<
590  fAtomicScreeningRadius[intZ-1] << " m_e*c/hbar --> BCB = " <<
591  matRadius << G4endl;
592  G4cout << "Screening parameters F0 for material " << material->GetName() << " = " <<
593  f0a << "," << f0b << G4endl;
594  }
595  return;
596 }
597 
598 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
599 
600 std::pair<G4double,G4double>
601 G4PenelopeGammaConversionModel::GetScreeningFunctions(G4double B)
602 {
603  // This is subroutine SCHIFF of Penelope
604  //
605  // Screening Functions F1(B) and F2(B) in the Bethe-Heitler differential cross
606  // section for pair production
607  //
608  std::pair<G4double,G4double> result(0.,0.);
609  G4double BSquared = B*B;
610  G4double f1 = 2.0-2.0*std::log(1.0+BSquared);
611  G4double f2 = f1 - 6.66666666e-1; // (-2/3)
612  if (B < 1.0e-10)
613  f1 = f1-twopi*B;
614  else
615  {
616  G4double a0 = 4.0*B*std::atan(1./B);
617  f1 = f1 - a0;
618  f2 += 2.0*BSquared*(4.0-a0-3.0*std::log((1.0+BSquared)/BSquared));
619  }
620  G4double g1 = 0.5*(3.0*f1-f2);
621  G4double g2 = 0.25*(3.0*f1+f2);
622 
623  result.first = g1;
624  result.second = g2;
625 
626  return result;
627 }
Generated on Sat May 25 2013 14:33:34 for Geant4 by   doxygen 1.8.4