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