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