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