source: trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08GammaConversionModel.cc @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

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