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