Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4PenelopeOscillatorManager.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 // Authors: Luciano Pandola (luciano.pandola at lngs.infn.it)
27 //
28 // History:
29 // -----------
30 //
31 // 03 Dec 2009 First working version, Luciano Pandola
32 // 16 Feb 2010 Added methods to store also total Z and A for the
33 // molecule, Luciano Pandola
34 // 19 Feb 2010 Scale the Hartree factors in the Compton Oscillator
35 // table by (1/fine_structure_const), since the models use
36 // always the ratio (hartreeFactor/fine_structure_const)
37 // 16 Mar 2010 Added methods to calculate and store mean exc energy
38 // and plasma energy (used for Ionisation). L Pandola
39 // 18 Mar 2010 Added method to retrieve number of atoms per
40 // molecule. L. Pandola
41 // 06 Sep 2011 Override the local Penelope database and use the main
42 // G4AtomicDeexcitation database to retrieve the shell
43 // binding energies. L. Pandola
44 // 15 Mar 2012 Added method to retrieve number of atom of given Z per
45 // molecule. Restore the original Penelope database for levels
46 // below 100 eV. L. Pandola
47 //
48 // -------------------------------------------------------------------
49 
51 
52 #include "globals.hh"
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
56 #include "G4AtomicShell.hh"
57 #include "G4Material.hh"
58 #include "G4Exp.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 
63  oscillatorStoreIonisation(0),oscillatorStoreCompton(0),atomicNumber(0),
64  atomicMass(0),excitationEnergy(0),plasmaSquared(0),atomsPerMolecule(0),
65  atomTablePerMolecule(0)
66 {
67  fReadElementData = false;
68  for (G4int i=0;i<5;i++)
69  {
70  for (G4int j=0;j<2000;j++)
71  elementData[i][j] = 0.;
72  }
73  verbosityLevel = 0;
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79 {
80  Clear();
81  delete instance;
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
91 {
92  if (!instance)
94  return instance;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
100 {
101  if (verbosityLevel > 1)
102  G4cout << " G4PenelopeOscillatorManager::Clear() - Clean Oscillator Tables" << G4endl;
103 
104  //Clean up OscillatorStoreIonisation
105  for (auto& item : (*oscillatorStoreIonisation))
106  {
107  G4PenelopeOscillatorTable* table = item.second;
108  if (table)
109  {
110  for (size_t k=0;k<table->size();k++) //clean individual oscillators
111  {
112  if ((*table)[k])
113  delete ((*table)[k]);
114  }
115  delete table;
116  }
117  }
119 
120  //Clean up OscillatorStoreCompton
121  for (auto& item : (*oscillatorStoreCompton))
122  {
123  G4PenelopeOscillatorTable* table = item.second;
124  if (table)
125  {
126  for (size_t k=0;k<table->size();k++) //clean individual oscillators
127  {
128  if ((*table)[k])
129  delete ((*table)[k]);
130  }
131  delete table;
132  }
133  }
134  delete oscillatorStoreCompton;
135 
136  if (atomicMass) delete atomicMass;
137  if (atomicNumber) delete atomicNumber;
139  if (plasmaSquared) delete plasmaSquared;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 
147 {
149  if (!theTable)
150  {
151  G4cout << " G4PenelopeOscillatorManager::Dump " << G4endl;
152  G4cout << "Problem in retrieving the Ionisation Oscillator Table for " << material->GetName() << G4endl;
153  return;
154  }
155  G4cout << "*********************************************************************" << G4endl;
156  G4cout << " Penelope Oscillator Table Ionisation for " << material->GetName() << G4endl;
157  G4cout << "*********************************************************************" << G4endl;
158  G4cout << "The table contains " << theTable->size() << " oscillators " << G4endl;
159  G4cout << "*********************************************************************" << G4endl;
160  if (theTable->size() < 10)
161  for (size_t k=0;k<theTable->size();k++)
162  {
163  G4cout << "Oscillator # " << k << " Z = " << (*theTable)[k]->GetParentZ() <<
164  " Shell Flag = " << (*theTable)[k]->GetShellFlag() <<
165  " Parent shell ID = " << (*theTable)[k]->GetParentShellID() << G4endl;
166  G4cout << "Ionisation energy = " << (*theTable)[k]->GetIonisationEnergy()/eV << " eV" << G4endl;
167  G4cout << "Occupation number = " << (*theTable)[k]->GetOscillatorStrength() << G4endl;
168  G4cout << "Resonance energy = " << (*theTable)[k]->GetResonanceEnergy()/eV << " eV" << G4endl;
169  G4cout << "Cufoff resonance energy = " <<
170  (*theTable)[k]->GetCutoffRecoilResonantEnergy()/eV << " eV" << G4endl;
171  G4cout << "*********************************************************************" << G4endl;
172  }
173  for (size_t k=0;k<theTable->size();k++)
174  {
175  G4cout << k << " " << (*theTable)[k]->GetOscillatorStrength() << " " <<
176  (*theTable)[k]->GetIonisationEnergy()/eV << " " << (*theTable)[k]->GetResonanceEnergy()/eV << " " <<
177  (*theTable)[k]->GetParentZ() << " " << (*theTable)[k]->GetShellFlag() << " " <<
178  (*theTable)[k]->GetParentShellID() << G4endl;
179  }
180  G4cout << "*********************************************************************" << G4endl;
181 
182 
183  //Compton table
184  theTable = GetOscillatorTableCompton(material);
185  if (!theTable)
186  {
187  G4cout << " G4PenelopeOscillatorManager::Dump " << G4endl;
188  G4cout << "Problem in retrieving the Compton Oscillator Table for " << material->GetName() << G4endl;
189  return;
190  }
191  G4cout << "*********************************************************************" << G4endl;
192  G4cout << " Penelope Oscillator Table Compton for " << material->GetName() << G4endl;
193  G4cout << "*********************************************************************" << G4endl;
194  G4cout << "The table contains " << theTable->size() << " oscillators " << G4endl;
195  G4cout << "*********************************************************************" << G4endl;
196  if (theTable->size() < 10)
197  for (size_t k=0;k<theTable->size();k++)
198  {
199  G4cout << "Oscillator # " << k << " Z = " << (*theTable)[k]->GetParentZ() <<
200  " Shell Flag = " << (*theTable)[k]->GetShellFlag() <<
201  " Parent shell ID = " << (*theTable)[k]->GetParentShellID() << G4endl;
202  G4cout << "Compton index = " << (*theTable)[k]->GetHartreeFactor() << G4endl;
203  G4cout << "Ionisation energy = " << (*theTable)[k]->GetIonisationEnergy()/eV << " eV" << G4endl;
204  G4cout << "Occupation number = " << (*theTable)[k]->GetOscillatorStrength() << G4endl;
205  G4cout << "*********************************************************************" << G4endl;
206  }
207  for (size_t k=0;k<theTable->size();k++)
208  {
209  G4cout << k << " " << (*theTable)[k]->GetOscillatorStrength() << " " <<
210  (*theTable)[k]->GetIonisationEnergy()/eV << " " << (*theTable)[k]->GetHartreeFactor() << " " <<
211  (*theTable)[k]->GetParentZ() << " " << (*theTable)[k]->GetShellFlag() << " " <<
212  (*theTable)[k]->GetParentShellID() << G4endl;
213  }
214  G4cout << "*********************************************************************" << G4endl;
215 
216  return;
217 }
218 
219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
220 
222 {
223  //Tables should be created at the same time, since they are both filled
224  //simultaneously
226  {
227  oscillatorStoreIonisation = new std::map<const G4Material*,G4PenelopeOscillatorTable*>;
228  if (!fReadElementData)
229  ReadElementData();
231  //It should be ok now
232  G4Exception("G4PenelopeOscillatorManager::GetOscillatorTableIonisation()",
233  "em2034",FatalException,
234  "Problem in allocating the Oscillator Store for Ionisation");
235  }
236 
238  {
239  oscillatorStoreCompton = new std::map<const G4Material*,G4PenelopeOscillatorTable*>;
240  if (!fReadElementData)
241  ReadElementData();
243  //It should be ok now
244  G4Exception("G4PenelopeOscillatorManager::GetOscillatorTableIonisation()",
245  "em2034",FatalException,
246  "Problem in allocating the Oscillator Store for Compton");
247  }
248 
249  if (!atomicNumber)
250  atomicNumber = new std::map<const G4Material*,G4double>;
251  if (!atomicMass)
252  atomicMass = new std::map<const G4Material*,G4double>;
253  if (!excitationEnergy)
254  excitationEnergy = new std::map<const G4Material*,G4double>;
255  if (!plasmaSquared)
256  plasmaSquared = new std::map<const G4Material*,G4double>;
257  if (!atomsPerMolecule)
258  atomsPerMolecule = new std::map<const G4Material*,G4double>;
260  atomTablePerMolecule = new std::map< std::pair<const G4Material*,G4int>, G4double>;
261 }
262 
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
265 
267 {
268  // (1) First time, create oscillatorStores and read data
270 
271  // (2) Check if the material has been already included
272  if (atomicNumber->count(mat))
273  return atomicNumber->find(mat)->second;
274 
275  // (3) If we are here, it means that we have to create the table for the material
277 
278  // (4) now, the oscillator store should be ok
279  if (atomicNumber->count(mat))
280  return atomicNumber->find(mat)->second;
281  else
282  {
283  G4cout << "G4PenelopeOscillatorManager::GetTotalZ() " << G4endl;
284  G4cout << "Impossible to retrieve the total Z for " << mat->GetName() << G4endl;
285  return 0;
286  }
287 }
288 
289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
290 
292 {
293  // (1) First time, create oscillatorStores and read data
295 
296  // (2) Check if the material has been already included
297  if (atomicMass->count(mat))
298  return atomicMass->find(mat)->second;
299 
300  // (3) If we are here, it means that we have to create the table for the material
302 
303  // (4) now, the oscillator store should be ok
304  if (atomicMass->count(mat))
305  return atomicMass->find(mat)->second;
306  else
307  {
308  G4cout << "G4PenelopeOscillatorManager::GetTotalA() " << G4endl;
309  G4cout << "Impossible to retrieve the total A for " << mat->GetName() << G4endl;
310  return 0;
311  }
312 }
313 
314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315 
317 {
318  // (1) First time, create oscillatorStores and read data
320 
321  // (2) Check if the material has been already included
322  if (oscillatorStoreIonisation->count(mat))
323  {
324  //Ok, it exists
325  return oscillatorStoreIonisation->find(mat)->second;
326  }
327 
328  // (3) If we are here, it means that we have to create the table for the material
330 
331  // (4) now, the oscillator store should be ok
332  if (oscillatorStoreIonisation->count(mat))
333  return oscillatorStoreIonisation->find(mat)->second;
334  else
335  {
336  G4cout << "G4PenelopeOscillatorManager::GetOscillatorTableIonisation() " << G4endl;
337  G4cout << "Impossible to create ionisation oscillator table for " << mat->GetName() << G4endl;
338  return nullptr;
339  }
340 }
341 
342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
343 
345  G4int index)
346 {
348  if (((size_t)index) < theTable->size())
349  return (*theTable)[index];
350  else
351  {
352  G4cout << "WARNING: Ionisation table for material " << material->GetName() << " has " <<
353  theTable->size() << " oscillators" << G4endl;
354  G4cout << "Oscillator #" << index << " cannot be retrieved" << G4endl;
355  G4cout << "Returning null pointer" << G4endl;
356  return nullptr;
357  }
358 }
359 
360 
361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
362 
364 {
365  // (1) First time, create oscillatorStore and read data
367 
368  // (2) Check if the material has been already included
369  if (oscillatorStoreCompton->count(mat))
370  {
371  //Ok, it exists
372  return oscillatorStoreCompton->find(mat)->second;
373  }
374 
375  // (3) If we are here, it means that we have to create the table for the material
377 
378  // (4) now, the oscillator store should be ok
379  if (oscillatorStoreCompton->count(mat))
380  return oscillatorStoreCompton->find(mat)->second;
381  else
382  {
383  G4cout << "G4PenelopeOscillatorManager::GetOscillatorTableCompton() " << G4endl;
384  G4cout << "Impossible to create Compton oscillator table for " << mat->GetName() << G4endl;
385  return nullptr;
386  }
387 }
388 
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390 
392  G4int index)
393 {
395  if (((size_t)index) < theTable->size())
396  return (*theTable)[index];
397  else
398  {
399  G4cout << "WARNING: Compton table for material " << material->GetName() << " has " <<
400  theTable->size() << " oscillators" << G4endl;
401  G4cout << "Oscillator #" << index << " cannot be retrieved" << G4endl;
402  G4cout << "Returning null pointer" << G4endl;
403  return nullptr;
404  }
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408 
410 {
411  //THIS CORRESPONDS TO THE ROUTINE PEMATW of PENELOPE
412 
413  G4double meanAtomExcitationEnergy[99] = {19.2*eV, 41.8*eV, 40.0*eV, 63.7*eV, 76.0*eV, 81.0*eV,
414  82.0*eV, 95.0*eV,115.0*eV,137.0*eV,149.0*eV,156.0*eV,
415  166.0*eV,
416  173.0*eV,173.0*eV,180.0*eV,174.0*eV,188.0*eV,190.0*eV,191.0*eV,
417  216.0*eV,233.0*eV,245.0*eV,257.0*eV,272.0*eV,286.0*eV,297.0*eV,
418  311.0*eV,322.0*eV,330.0*eV,334.0*eV,350.0*eV,347.0*eV,348.0*eV,
419  343.0*eV,352.0*eV,363.0*eV,366.0*eV,379.0*eV,393.0*eV,417.0*eV,
420  424.0*eV,428.0*eV,441.0*eV,449.0*eV,470.0*eV,470.0*eV,469.0*eV,
421  488.0*eV,488.0*eV,487.0*eV,485.0*eV,491.0*eV,482.0*eV,488.0*eV,
422  491.0*eV,501.0*eV,523.0*eV,535.0*eV,546.0*eV,560.0*eV,574.0*eV,
423  580.0*eV,591.0*eV,614.0*eV,628.0*eV,650.0*eV,658.0*eV,674.0*eV,
424  684.0*eV,694.0*eV,705.0*eV,718.0*eV,727.0*eV,736.0*eV,746.0*eV,
425  757.0*eV,790.0*eV,790.0*eV,800.0*eV,810.0*eV,823.0*eV,823.0*eV,
426  830.0*eV,825.0*eV,794.0*eV,827.0*eV,826.0*eV,841.0*eV,847.0*eV,
427  878.0*eV,890.0*eV,902.0*eV,921.0*eV,934.0*eV,939.0*eV,952.0*eV,
428  966.0*eV,980.0*eV};
429 
430  if (verbosityLevel > 0)
431  G4cout << "Going to build Oscillator Table for " << material->GetName() << G4endl;
432 
433  G4int nElements = material->GetNumberOfElements();
434  const G4ElementVector* elementVector = material->GetElementVector();
435 
436 
437  //At the moment, there's no way in Geant4 to know if a material
438  //is defined with atom numbers or fraction of weigth
439  const G4double* fractionVector = material->GetFractionVector();
440 
441 
442  //Take always the composition by fraction of mass. For the composition by
443  //atoms: it is calculated by Geant4 but with some rounding to integers
444  G4double totalZ = 0;
445  G4double totalMolecularWeight = 0;
446  G4double meanExcitationEnergy = 0;
447 
448  std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
449 
450  for (G4int i=0;i<nElements;i++)
451  {
452  //G4int iZ = (G4int) (*elementVector)[i]->GetZ();
453  G4double fraction = fractionVector[i];
454  G4double atomicWeigth = (*elementVector)[i]->GetAtomicMassAmu();
455  StechiometricFactors->push_back(fraction/atomicWeigth);
456  }
457  //Find max
458  G4double MaxStechiometricFactor = 0.;
459  for (G4int i=0;i<nElements;i++)
460  {
461  if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
462  MaxStechiometricFactor = (*StechiometricFactors)[i];
463  }
464  if (MaxStechiometricFactor<1e-16)
465  {
467  ed << "Problem with the mass composition of " << material->GetName() << G4endl;
468  ed << "MaxStechiometricFactor = " << MaxStechiometricFactor << G4endl;
469  G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
470  "em2035",FatalException,ed);
471  }
472  //Normalize
473  for (G4int i=0;i<nElements;i++)
474  (*StechiometricFactors)[i] /= MaxStechiometricFactor;
475 
476  // Equivalent atoms per molecule
477  G4double theatomsPerMolecule = 0;
478  for (G4int i=0;i<nElements;i++)
479  theatomsPerMolecule += (*StechiometricFactors)[i];
480  G4double moleculeDensity =
481  material->GetTotNbOfAtomsPerVolume()/theatomsPerMolecule; //molecules per unit volume
482 
483 
484  if (verbosityLevel > 1)
485  {
486  for (size_t i=0;i<StechiometricFactors->size();i++)
487  {
488  G4cout << "Element " << (*elementVector)[i]->GetSymbol() << " (Z = " <<
489  (*elementVector)[i]->GetZ() << ") --> " <<
490  (*StechiometricFactors)[i] << " atoms/molecule " << G4endl;
491  }
492  }
493 
494 
495  for (G4int i=0;i<nElements;i++)
496  {
497  G4int iZ = (G4int) (*elementVector)[i]->GetZ();
498  totalZ += iZ * (*StechiometricFactors)[i];
499  totalMolecularWeight += (*elementVector)[i]->GetAtomicMassAmu() * (*StechiometricFactors)[i];
500  meanExcitationEnergy += iZ*std::log(meanAtomExcitationEnergy[iZ-1])*(*StechiometricFactors)[i];
501  /*
502  G4cout << iZ << " " << (*StechiometricFactors)[i] << " " << totalZ << " " <<
503  totalMolecularWeight/(g/mole) << " " << meanExcitationEnergy << " " <<
504  meanAtomExcitationEnergy[iZ-1]/eV <<
505  G4endl;
506  */
507  std::pair<const G4Material*,G4int> theKey = std::make_pair(material,iZ);
508  if (!atomTablePerMolecule->count(theKey))
509  atomTablePerMolecule->insert(std::make_pair(theKey,(*StechiometricFactors)[i]));
510  }
511  meanExcitationEnergy = G4Exp(meanExcitationEnergy/totalZ);
512 
513  atomicNumber->insert(std::make_pair(material,totalZ));
514  atomicMass->insert(std::make_pair(material,totalMolecularWeight));
515  excitationEnergy->insert(std::make_pair(material,meanExcitationEnergy));
516  atomsPerMolecule->insert(std::make_pair(material,theatomsPerMolecule));
517 
518 
519  if (verbosityLevel > 1)
520  {
521  G4cout << "Calculated mean excitation energy for " << material->GetName() <<
522  " = " << meanExcitationEnergy/eV << " eV" << G4endl;
523  }
524 
525  std::vector<G4PenelopeOscillator> *helper = new std::vector<G4PenelopeOscillator>;
526 
527  //First Oscillator: conduction band. Tentativaly assumed to consist of valence electrons (each
528  //atom contributes a number of electrons equal to its lowest chemical valence)
529  G4PenelopeOscillator newOsc;
530  newOsc.SetOscillatorStrength(0.);
531  newOsc.SetIonisationEnergy(0*eV);
532  newOsc.SetHartreeFactor(0);
533  newOsc.SetParentZ(0);
534  newOsc.SetShellFlag(30);
535  newOsc.SetParentShellID(30); //does not correspond to any "real" level
536  helper->push_back(newOsc);
537 
538  //Load elements and oscillators
539  for (G4int k=0;k<nElements;k++)
540  {
541  G4double Z = (*elementVector)[k]->GetZ();
542  G4bool finished = false;
543  for (G4int i=0;i<2000 && !finished;i++)
544  {
545  /*
546  elementData[0][i] = Z;
547  elementData[1][i] = shellCode;
548  elementData[2][i] = occupationNumber;
549  elementData[3][i] = ionisationEnergy;
550  elementData[4][i] = hartreeProfile;
551  */
552  if (elementData[0][i] == Z)
553  {
554  G4int shellID = (G4int) elementData[1][i];
555  G4double occup = elementData[2][i];
556  if (shellID > 0)
557  {
558 
559  if (std::fabs(occup) > 0)
560  {
561  G4PenelopeOscillator newOscLocal;
562  newOscLocal.SetOscillatorStrength(std::fabs(occup)*(*StechiometricFactors)[k]);
563  newOscLocal.SetIonisationEnergy(elementData[3][i]);
565  newOscLocal.SetParentZ(elementData[0][i]);
566  //keep track of the origianl shell level
567  newOscLocal.SetParentShellID((G4int)elementData[1][i]);
568  //register only K, L and M shells. Outer shells all grouped with
569  //shellIndex = 30
570  if (elementData[0][i] > 6 && elementData[1][i] < 10)
571  newOscLocal.SetShellFlag(((G4int)elementData[1][i]));
572  else
573  newOscLocal.SetShellFlag(30);
574  helper->push_back(newOscLocal);
575  if (occup < 0)
576  {
577  G4double ff = (*helper)[0].GetOscillatorStrength();
578  ff += std::fabs(occup)*(*StechiometricFactors)[k];
579  (*helper)[0].SetOscillatorStrength(ff);
580  }
581  }
582  }
583  }
584  if (elementData[0][i] > Z)
585  finished = true;
586  }
587  }
588 
589  delete StechiometricFactors;
590 
591  //NOW: sort oscillators according to increasing ionisation energy
592  //Notice: it works because helper is a vector of _object_, not a
593  //vector to _pointers_
594  std::sort(helper->begin(),helper->end());
595 
596  // Plasma energy and conduction band excitation
597  static const G4double RydbergEnergy = 13.60569*eV;
598  G4double Omega = std::sqrt(4*pi*moleculeDensity*totalZ*Bohr_radius)*Bohr_radius*2.0*RydbergEnergy;
599  G4double conductionStrength = (*helper)[0].GetOscillatorStrength();
600  G4double plasmaEnergy = Omega*std::sqrt(conductionStrength/totalZ);
601 
602  plasmaSquared->insert(std::make_pair(material,Omega*Omega));
603 
604  G4bool isAConductor = false;
605  G4int nullOsc = 0;
606 
607  if (verbosityLevel > 1)
608  {
609  G4cout << "Estimated oscillator strenght and energy of plasmon: " <<
610  conductionStrength << " and " << plasmaEnergy/eV << " eV" << G4endl;
611  }
612 
613  if (conductionStrength < 0.01 || plasmaEnergy<1.0*eV) //this is an insulator
614  {
615  if (verbosityLevel >1 )
616  G4cout << material->GetName() << " is an insulator " << G4endl;
617  //remove conduction band oscillator
618  helper->erase(helper->begin());
619  }
620  else //this is a conductor, Outer shells moved to conduction band
621  {
622  if (verbosityLevel >1 )
623  G4cout << material->GetName() << " is a conductor " << G4endl;
624  isAConductor = true;
625  //copy the conduction strenght.. The number is going to change.
626  G4double conductionStrengthCopy = conductionStrength;
627  G4bool quit = false;
628  for (size_t i = 1; i<helper->size() && !quit ;i++)
629  {
630  G4double oscStre = (*helper)[i].GetOscillatorStrength();
631  //loop is repeated over here
632  if (oscStre < conductionStrengthCopy)
633  {
634  conductionStrengthCopy = conductionStrengthCopy-oscStre;
635  (*helper)[i].SetOscillatorStrength(0.);
636  nullOsc++;
637  }
638  else //this is passed only once - no goto -
639  {
640  quit = true;
641  (*helper)[i].SetOscillatorStrength(oscStre-conductionStrengthCopy);
642  if (std::fabs((*helper)[i].GetOscillatorStrength()) < 1e-12)
643  {
644  conductionStrength += (*helper)[i].GetOscillatorStrength();
645  (*helper)[i].SetOscillatorStrength(0.);
646  nullOsc++;
647  }
648  }
649  }
650 
651  //Update conduction band
652  (*helper)[0].SetOscillatorStrength(conductionStrength);
653  (*helper)[0].SetIonisationEnergy(0.);
654  (*helper)[0].SetResonanceEnergy(plasmaEnergy);
655  G4double hartree = 0.75/std::sqrt(3.0*pi*pi*moleculeDensity*
656  Bohr_radius*Bohr_radius*Bohr_radius*conductionStrength);
657  (*helper)[0].SetHartreeFactor(hartree/fine_structure_const);
658  }
659 
660  //Check f-sum rule
661  G4double sum = 0;
662  for (size_t i=0;i<helper->size();i++)
663  {
664  sum += (*helper)[i].GetOscillatorStrength();
665  }
666  if (std::fabs(sum-totalZ) > (1e-6*totalZ))
667  {
669  ed << "Inconsistent oscillator data for " << material->GetName() << G4endl;
670  ed << sum << " " << totalZ << G4endl;
671  G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
672  "em2036",FatalException,ed);
673  }
674  if (std::fabs(sum-totalZ) > (1e-12*totalZ))
675  {
676  G4double fact = totalZ/sum;
677  for (size_t i=0;i<helper->size();i++)
678  {
679  G4double ff = (*helper)[i].GetOscillatorStrength()*fact;
680  (*helper)[i].SetOscillatorStrength(ff);
681  }
682  }
683 
684  //Remove null items
685  for (G4int k=0;k<nullOsc;k++)
686  {
687  G4bool exit=false;
688  for (size_t i=0;i<helper->size() && !exit;i++)
689  {
690  if (std::fabs((*helper)[i].GetOscillatorStrength()) < 1e-12)
691  {
692  helper->erase(helper->begin()+i);
693  exit = true;
694  }
695  }
696  }
697 
698  //Sternheimer's adjustment factor
699  G4double adjustmentFactor = 0;
700  if (helper->size() > 1)
701  {
702  G4double TST = totalZ*std::log(meanExcitationEnergy/eV);
703  G4double AALow = 0.1;
704  G4double AAHigh = 10.;
705  do
706  {
707  adjustmentFactor = (AALow+AAHigh)*0.5;
708  G4double sumLocal = 0;
709  for (size_t i=0;i<helper->size();i++)
710  {
711  if (i == 0 && isAConductor)
712  {
713  G4double resEne = (*helper)[i].GetResonanceEnergy();
714  sumLocal += (*helper)[i].GetOscillatorStrength()*std::log(resEne/eV);
715  }
716  else
717  {
718  G4double ionEne = (*helper)[i].GetIonisationEnergy();
719  G4double oscStre = (*helper)[i].GetOscillatorStrength();
720  G4double WI2 = (adjustmentFactor*adjustmentFactor*ionEne*ionEne) +
721  2./3.*(oscStre/totalZ)*Omega*Omega;
722  G4double resEne = std::sqrt(WI2);
723  (*helper)[i].SetResonanceEnergy(resEne);
724  sumLocal += (*helper)[i].GetOscillatorStrength()*std::log(resEne/eV);
725  }
726  }
727  if (sumLocal < TST)
728  AALow = adjustmentFactor;
729  else
730  AAHigh = adjustmentFactor;
731  if (verbosityLevel > 3)
732  G4cout << "Sternheimer's adjustment factor loops: " << AALow << " " << AAHigh << " " <<
733  adjustmentFactor << " " << TST << " " <<
734  sumLocal << G4endl;
735  }while((AAHigh-AALow)>(1e-14*adjustmentFactor));
736  }
737  else
738  {
739  G4double ionEne = (*helper)[0].GetIonisationEnergy();
740  (*helper)[0].SetIonisationEnergy(std::fabs(ionEne));
741  (*helper)[0].SetResonanceEnergy(meanExcitationEnergy);
742  }
743  if (verbosityLevel > 1)
744  {
745  G4cout << "Sternheimer's adjustment factor: " << adjustmentFactor << G4endl;
746  }
747 
748  //Check again for data consistency
749  G4double xcheck = (*helper)[0].GetOscillatorStrength()*std::log((*helper)[0].GetResonanceEnergy());
750  G4double TST = (*helper)[0].GetOscillatorStrength();
751  for (size_t i=1;i<helper->size();i++)
752  {
753  xcheck += (*helper)[i].GetOscillatorStrength()*std::log((*helper)[i].GetResonanceEnergy());
754  TST += (*helper)[i].GetOscillatorStrength();
755  }
756  if (std::fabs(TST-totalZ)>1e-8*totalZ)
757  {
759  ed << "Inconsistent oscillator data " << G4endl;
760  ed << TST << " " << totalZ << G4endl;
761  G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
762  "em2036",FatalException,ed);
763  }
764  xcheck = G4Exp(xcheck/totalZ);
765  if (std::fabs(xcheck-meanExcitationEnergy) > 1e-8*meanExcitationEnergy)
766  {
768  ed << "Error in Sterheimer factor calculation " << G4endl;
769  ed << xcheck/eV << " " << meanExcitationEnergy/eV << G4endl;
770  G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
771  "em2037",FatalException,ed);
772  }
773 
774  //Selection of the lowest ionisation energy for inner shells. Only the K, L and M shells with
775  //ionisation energy less than the N1 shell of the heaviest element in the material are considered as
776  //inner shells. As a results, the inner/outer shell character of an atomic shell depends on the
777  //composition of the material.
778  G4double Zmax = 0;
779  for (G4int k=0;k<nElements;k++)
780  {
781  G4double Z = (*elementVector)[k]->GetZ();
782  if (Z>Zmax) Zmax = Z;
783  }
784  //Find N1 level of the heaviest element (if any).
785  G4bool found = false;
786  G4double cutEnergy = 50*eV;
787  for (size_t i=0;i<helper->size() && !found;i++)
788  {
789  G4double Z = (*helper)[i].GetParentZ();
790  G4int shID = (*helper)[i].GetParentShellID(); //look for the N1 level
791  if (shID == 10 && Z == Zmax)
792  {
793  found = true;
794  if ((*helper)[i].GetIonisationEnergy() > cutEnergy)
795  cutEnergy = (*helper)[i].GetIonisationEnergy();
796  }
797  }
798  //Make that cutEnergy cannot be higher than 250 eV, namely the fluorescence level by
799  //Geant4
800  G4double lowEnergyLimitForFluorescence = 250*eV;
801  cutEnergy = std::min(cutEnergy,lowEnergyLimitForFluorescence);
802 
803  if (verbosityLevel > 1)
804  G4cout << "Cutoff energy: " << cutEnergy/eV << " eV" << G4endl;
805  //
806  //Copy helper in the oscillatorTable for Ionisation
807  //
808  //Oscillator table Ionisation for the material
809  G4PenelopeOscillatorTable* theTable = new G4PenelopeOscillatorTable(); //vector of oscillator
811  std::sort(helper->begin(),helper->end(),comparator);
812 
813  //COPY THE HELPER (vector of object) to theTable (vector of Pointers).
814  for (size_t i=0;i<helper->size();i++)
815  {
816  //copy content --> one may need it later (e.g. to fill an other table, with variations)
817  G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
818  theTable->push_back(theOsc);
819  }
820 
821  //Oscillators of outer shells with resonance energies differing by a factor less than
822  //Rgroup are grouped as a single oscillator
823  G4double Rgroup = 1.05;
824  size_t Nost = theTable->size();
825 
826  size_t firstIndex = (isAConductor) ? 1 : 0; //for conductors, skip conduction oscillator
827  G4bool loopAgain = false;
828  G4int nLoops = 0;
829  G4int removedLevels = 0;
830  do
831  {
832  loopAgain = false;
833  nLoops++;
834  if (Nost>firstIndex+1)
835  {
836  removedLevels = 0;
837  for (size_t i=firstIndex;i<theTable->size()-1;i++)
838  {
839  G4bool skipLoop = false;
840  G4int shellFlag = (*theTable)[i]->GetShellFlag();
841  G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
842  G4double resEne = (*theTable)[i]->GetResonanceEnergy();
843  G4double resEnePlus1 = (*theTable)[i+1]->GetResonanceEnergy();
844  G4double oscStre = (*theTable)[i]->GetOscillatorStrength();
845  G4double oscStrePlus1 = (*theTable)[i+1]->GetOscillatorStrength();
846  //if (shellFlag < 10 && ionEne>cutEnergy) in Penelope
847  if (ionEne>cutEnergy) //remove condition that shellFlag < 10!
848  skipLoop = true;
849  if (resEne<1.0*eV || resEnePlus1<1.0*eV)
850  skipLoop = true;
851  if (resEnePlus1 > Rgroup*resEne)
852  skipLoop = true;
853  if (!skipLoop)
854  {
855  G4double newRes = G4Exp((oscStre*std::log(resEne)+
856  oscStrePlus1*std::log(resEnePlus1))
857  /(oscStre+oscStrePlus1));
858  (*theTable)[i]->SetResonanceEnergy(newRes);
859  G4double newIon = (oscStre*ionEne+
860  oscStrePlus1*(*theTable)[i+1]->GetIonisationEnergy())/
861  (oscStre+oscStrePlus1);
862  (*theTable)[i]->SetIonisationEnergy(newIon);
863  G4double newStre = oscStre+oscStrePlus1;
864  (*theTable)[i]->SetOscillatorStrength(newStre);
865  G4double newHartree = (oscStre*(*theTable)[i]->GetHartreeFactor()+
866  oscStrePlus1*(*theTable)[i+1]->GetHartreeFactor())/
867  (oscStre+oscStrePlus1);
868  (*theTable)[i]->SetHartreeFactor(newHartree);
869  if ((*theTable)[i]->GetParentZ() != (*theTable)[i+1]->GetParentZ())
870  (*theTable)[i]->SetParentZ(0.);
871  if (shellFlag < 10 || (*theTable)[i+1]->GetShellFlag() < 10)
872  {
873  G4int newFlag = std::min(shellFlag,(*theTable)[i+1]->GetShellFlag());
874  (*theTable)[i]->SetShellFlag(newFlag);
875  }
876  else
877  (*theTable)[i]->SetShellFlag(30);
878  //We've lost anyway the track of the original level
879  (*theTable)[i]->SetParentShellID((*theTable)[i]->GetShellFlag());
880 
881 
882  if (i<theTable->size()-2)
883  {
884  for (size_t ii=i+1;ii<theTable->size()-1;ii++)
885  (*theTable)[ii] = (*theTable)[ii+1];
886  }
887  //G4cout << theTable->size() << G4endl;
888  theTable->erase(theTable->begin()+theTable->size()-1); //delete last element
889  removedLevels++;
890  }
891  }
892  }
893  if (removedLevels)
894  {
895  Nost -= removedLevels;
896  loopAgain = true;
897  }
898  if (Rgroup < 1.414213 || Nost > 64)
899  {
900  Rgroup = Rgroup*Rgroup;
901  loopAgain = true;
902  }
903  //Add protection against infinite loops here
904  if (nLoops > 100 && !removedLevels)
905  loopAgain = false;
906  }while(loopAgain);
907 
908  if (verbosityLevel > 1)
909  {
910  G4cout << "Final grouping factor for Ionisation: " << Rgroup << G4endl;
911  }
912 
913  //Final Electron/Positron model parameters
914  for (size_t i=0;i<theTable->size();i++)
915  {
916  //Set cutoff recoil energy for the resonant mode
917  G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
918  if (ionEne < 1e-3*eV)
919  {
920  G4double resEne = (*theTable)[i]->GetResonanceEnergy();
921  (*theTable)[i]->SetIonisationEnergy(0.*eV);
922  (*theTable)[i]->SetCutoffRecoilResonantEnergy(resEne);
923  }
924  else
925  (*theTable)[i]->SetCutoffRecoilResonantEnergy(ionEne);
926  }
927 
928  //Last step
929  oscillatorStoreIonisation->insert(std::make_pair(material,theTable));
930 
931 
932  /*
933  SAME FOR COMPTON
934  */
935  //
936  //Copy helper in the oscillatorTable for Compton
937  //
938  //Oscillator table Ionisation for the material
939  G4PenelopeOscillatorTable* theTableC = new G4PenelopeOscillatorTable(); //vector of oscillator
940  //order by ionisation energy
941  std::sort(helper->begin(),helper->end());
942  //COPY THE HELPER (vector of object) to theTable (vector of Pointers).
943  for (size_t i=0;i<helper->size();i++)
944  {
945  //copy content --> one may need it later (e.g. to fill an other table, with variations)
946  G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
947  theTableC->push_back(theOsc);
948  }
949  //Oscillators of outer shells with resonance energies differing by a factor less than
950  //Rgroup are grouped as a single oscillator
951  Rgroup = 1.5;
952  Nost = theTableC->size();
953 
954  firstIndex = (isAConductor) ? 1 : 0; //for conductors, skip conduction oscillator
955  loopAgain = false;
956  removedLevels = 0;
957  do
958  {
959  nLoops++;
960  loopAgain = false;
961  if (Nost>firstIndex+1)
962  {
963  removedLevels = 0;
964  for (size_t i=firstIndex;i<theTableC->size()-1;i++)
965  {
966  G4bool skipLoop = false;
967  //G4int shellFlag = (*theTableC)[i]->GetShellFlag();
968  G4double ionEne = (*theTableC)[i]->GetIonisationEnergy();
969  G4double ionEnePlus1 = (*theTableC)[i+1]->GetIonisationEnergy();
970  G4double oscStre = (*theTableC)[i]->GetOscillatorStrength();
971  G4double oscStrePlus1 = (*theTableC)[i+1]->GetOscillatorStrength();
972  //if (shellFlag < 10 && ionEne>cutEnergy) in Penelope
973  if (ionEne>cutEnergy)
974  skipLoop = true;
975  if (ionEne<1.0*eV || ionEnePlus1<1.0*eV)
976  skipLoop = true;
977  if (ionEnePlus1 > Rgroup*ionEne)
978  skipLoop = true;
979 
980  if (!skipLoop)
981  {
982  G4double newIon = (oscStre*ionEne+
983  oscStrePlus1*ionEnePlus1)/
984  (oscStre+oscStrePlus1);
985  (*theTableC)[i]->SetIonisationEnergy(newIon);
986  G4double newStre = oscStre+oscStrePlus1;
987  (*theTableC)[i]->SetOscillatorStrength(newStre);
988  G4double newHartree = (oscStre*(*theTableC)[i]->GetHartreeFactor()+
989  oscStrePlus1*(*theTableC)[i+1]->GetHartreeFactor())/
990  (oscStre+oscStrePlus1);
991  (*theTableC)[i]->SetHartreeFactor(newHartree);
992  if ((*theTableC)[i]->GetParentZ() != (*theTableC)[i+1]->GetParentZ())
993  (*theTableC)[i]->SetParentZ(0.);
994  (*theTableC)[i]->SetShellFlag(30);
995  (*theTableC)[i]->SetParentShellID((*theTableC)[i]->GetShellFlag());
996 
997  if (i<theTableC->size()-2)
998  {
999  for (size_t ii=i+1;ii<theTableC->size()-1;ii++)
1000  (*theTableC)[ii] = (*theTableC)[ii+1];
1001  }
1002  theTableC->erase(theTableC->begin()+theTableC->size()-1); //delete last element
1003  removedLevels++;
1004  }
1005  }
1006  }
1007  if (removedLevels)
1008  {
1009  Nost -= removedLevels;
1010  loopAgain = true;
1011  }
1012  if (Rgroup < 2.0 || Nost > 64)
1013  {
1014  Rgroup = Rgroup*Rgroup;
1015  loopAgain = true;
1016  }
1017  //Add protection against infinite loops here
1018  if (nLoops > 100 && !removedLevels)
1019  loopAgain = false;
1020  }while(loopAgain);
1021 
1022 
1023  if (verbosityLevel > 1)
1024  {
1025  G4cout << "Final grouping factor for Compton: " << Rgroup << G4endl;
1026  }
1027 
1028  //Last step
1029  oscillatorStoreCompton->insert(std::make_pair(material,theTableC));
1030 
1031  /* //TESTING PURPOSES
1032  if (verbosityLevel > 1)
1033  {
1034  G4cout << "The table contains " << helper->size() << " oscillators " << G4endl;
1035  for (size_t k=0;k<helper->size();k++)
1036  {
1037  G4cout << "Oscillator # " << k << G4endl;
1038  G4cout << "Z = " << (*helper)[k].GetParentZ() << G4endl;
1039  G4cout << "Shell Flag = " << (*helper)[k].GetShellFlag() << G4endl;
1040  G4cout << "Compton index = " << (*helper)[k].GetHartreeFactor() << G4endl;
1041  G4cout << "Ionisation energy = " << (*helper)[k].GetIonisationEnergy()/eV << " eV" << G4endl;
1042  G4cout << "Occupation number = " << (*helper)[k].GetOscillatorStrength() << G4endl;
1043  G4cout << "Resonance energy = " << (*helper)[k].GetResonanceEnergy()/eV << " eV" << G4endl;
1044  }
1045 
1046  for (size_t k=0;k<helper->size();k++)
1047  {
1048  G4cout << k << " " << (*helper)[k].GetOscillatorStrength() << " " <<
1049  (*helper)[k].GetIonisationEnergy()/eV << " " << (*helper)[k].GetResonanceEnergy()/eV << " " <<
1050  (*helper)[k].GetParentZ() << " " << (*helper)[k].GetShellFlag() << " " <<
1051  (*helper)[k].GetHartreeFactor() << G4endl;
1052  }
1053  }
1054  */
1055 
1056 
1057  //CLEAN UP theHelper and its content
1058  delete helper;
1059  if (verbosityLevel > 1)
1060  Dump(material);
1061 
1062  return;
1063 }
1064 
1065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1066 
1068 {
1069  if (verbosityLevel > 0)
1070  {
1071  G4cout << "G4PenelopeOscillatorManager::ReadElementData()" << G4endl;
1072  G4cout << "Going to read Element Data" << G4endl;
1073  }
1074  char* path = getenv("G4LEDATA");
1075  if (!path)
1076  {
1077  G4String excep = "G4PenelopeOscillatorManager - G4LEDATA environment variable not set!";
1078  G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
1079  "em0006",FatalException,excep);
1080  return;
1081  }
1082  G4String pathString(path);
1083  G4String pathFile = pathString + "/penelope/pdatconf.p08";
1084  std::ifstream file(pathFile);
1085 
1086  if (!file.is_open())
1087  {
1088  G4String excep = "G4PenelopeOscillatorManager - data file " + pathFile + " not found!";
1089  G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
1090  "em0003",FatalException,excep);
1091  }
1092 
1093  G4AtomicTransitionManager* theTransitionManager =
1095  theTransitionManager->Initialise();
1096 
1097 
1098  //Read header (22 lines)
1099  G4String theHeader;
1100  for (G4int iline=0;iline<22;iline++)
1101  getline(file,theHeader);
1102  //Done
1103  G4int Z=0;
1104  G4int shellCode = 0;
1105  G4String shellId = "NULL";
1106  G4int occupationNumber = 0;
1107  G4double ionisationEnergy = 0.0*eV;
1108  G4double hartreeProfile = 0.;
1109  G4int shellCounter = 0;
1110  G4int oldZ = -1;
1111  G4int numberOfShells = 0;
1112  //Start reading data
1113  for (G4int i=0;!file.eof();i++)
1114  {
1115  file >> Z >> shellCode >> shellId >> occupationNumber >> ionisationEnergy >> hartreeProfile;
1116  if (Z>0 && i<2000)
1117  {
1118  elementData[0][i] = Z;
1119  elementData[1][i] = shellCode;
1120  elementData[2][i] = occupationNumber;
1121  //reset things
1122  if (Z != oldZ)
1123  {
1124  shellCounter = 0;
1125  oldZ = Z;
1126  numberOfShells = theTransitionManager->NumberOfShells(Z);
1127  }
1128  G4double bindingEnergy = -1*eV;
1129  if (shellCounter<numberOfShells)
1130  {
1131  G4AtomicShell* shell = theTransitionManager->Shell(Z,shellCounter);
1132  bindingEnergy = shell->BindingEnergy();
1133  }
1134  //Valid level found in the G4AtomicTransition database: keep it, otherwise use
1135  //the ionisation energy found in the Penelope database
1136  elementData[3][i] = (bindingEnergy>100*eV) ? bindingEnergy : ionisationEnergy*eV;
1137  //elementData[3][i] = ionisationEnergy*eV;
1138  elementData[4][i] = hartreeProfile;
1139  shellCounter++;
1140  }
1141  }
1142  file.close();
1143 
1144  if (verbosityLevel > 1)
1145  {
1146  G4cout << "G4PenelopeOscillatorManager::ReadElementData(): Data file read" << G4endl;
1147  }
1148  fReadElementData = true;
1149  return;
1150 
1151 }
1152 
1153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1155 {
1156  // (1) First time, create oscillatorStores and read data
1158 
1159  // (2) Check if the material has been already included
1160  if (excitationEnergy->count(mat))
1161  return excitationEnergy->find(mat)->second;
1162 
1163  // (3) If we are here, it means that we have to create the table for the material
1164  BuildOscillatorTable(mat);
1165 
1166  // (4) now, the oscillator store should be ok
1167  if (excitationEnergy->count(mat))
1168  return excitationEnergy->find(mat)->second;
1169  else
1170  {
1171  G4cout << "G4PenelopeOscillatorManager::GetMolecularExcitationEnergy() " << G4endl;
1172  G4cout << "Impossible to retrieve the excitation energy for " << mat->GetName() << G4endl;
1173  return 0;
1174  }
1175 }
1176 
1177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1179 {
1180  // (1) First time, create oscillatorStores and read data
1182 
1183  // (2) Check if the material has been already included
1184  if (plasmaSquared->count(mat))
1185  return plasmaSquared->find(mat)->second;
1186 
1187  // (3) If we are here, it means that we have to create the table for the material
1188  BuildOscillatorTable(mat);
1189 
1190  // (4) now, the oscillator store should be ok
1191  if (plasmaSquared->count(mat))
1192  return plasmaSquared->find(mat)->second;
1193  else
1194  {
1195  G4cout << "G4PenelopeOscillatorManager::GetPlasmaEnergySquared() " << G4endl;
1196  G4cout << "Impossible to retrieve the plasma energy for " << mat->GetName() << G4endl;
1197  return 0;
1198  }
1199 }
1200 
1201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1202 
1204 {
1205  // (1) First time, create oscillatorStores and read data
1207 
1208  // (2) Check if the material has been already included
1209  if (atomsPerMolecule->count(mat))
1210  return atomsPerMolecule->find(mat)->second;
1211 
1212  // (3) If we are here, it means that we have to create the table for the material
1213  BuildOscillatorTable(mat);
1214 
1215  // (4) now, the oscillator store should be ok
1216  if (atomsPerMolecule->count(mat))
1217  return atomsPerMolecule->find(mat)->second;
1218  else
1219  {
1220  G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
1221  G4cout << "Impossible to retrieve the number of atoms per molecule for "
1222  << mat->GetName() << G4endl;
1223  return 0;
1224  }
1225 }
1226 
1227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1228 
1230 {
1231  // (1) First time, create oscillatorStores and read data
1233 
1234  // (2) Check if the material/Z couple has been already included
1235  std::pair<const G4Material*,G4int> theKey = std::make_pair(mat,Z);
1236  if (atomTablePerMolecule->count(theKey))
1237  return atomTablePerMolecule->find(theKey)->second;
1238 
1239  // (3) If we are here, it means that we have to create the table for the material
1240  BuildOscillatorTable(mat);
1241 
1242  // (4) now, the oscillator store should be ok
1243  if (atomTablePerMolecule->count(theKey))
1244  return atomTablePerMolecule->find(theKey)->second;
1245  else
1246  {
1247  G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
1248  G4cout << "Impossible to retrieve the number of atoms per molecule for Z = "
1249  << Z << " in material " << mat->GetName() << G4endl;
1250  return 0;
1251  }
1252 }
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
void SetHartreeFactor(G4double hf)
G4PenelopeOscillator * GetOscillatorIonisation(const G4Material *, G4int)
const G4double * GetFractionVector() const
Definition: G4Material.hh:195
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4double GetMeanExcitationEnergy(const G4Material *)
Returns the mean excitation energy.
std::map< const G4Material *, G4double > * atomicNumber
#define G4endl
Definition: G4ios.hh:61
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
G4PenelopeOscillatorTable * GetOscillatorTableIonisation(const G4Material *)
#define G4ThreadLocal
Definition: tls.hh:69
G4double GetPlasmaEnergySquared(const G4Material *)
Returns the squared plasma energy.
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
void SetShellFlag(G4int theflag)
void SetParentZ(G4double parZ)
int Zmax
void SetOscillatorStrength(G4double ostr)
const G4String & GetName() const
Definition: G4Material.hh:179
Float_t Z
G4double GetTotalA(const G4Material *)
Returns the total A for the molecule.
static G4PenelopeOscillatorManager * GetOscillatorManager()
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
std::map< const G4Material *, G4PenelopeOscillatorTable * > * oscillatorStoreCompton
static G4ThreadLocal G4PenelopeOscillatorManager * instance
std::map< std::pair< const G4Material *, G4int >, G4double > * atomTablePerMolecule
std::map< const G4Material *, G4double > * atomicMass
void SetParentShellID(G4int psID)
static constexpr double eV
Definition: G4SIunits.hh:215
Float_t mat
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
std::vector< G4Element * > G4ElementVector
int G4int
Definition: G4Types.hh:78
G4PenelopeOscillator * GetOscillatorCompton(const G4Material *, G4int)
std::map< const G4Material *, G4PenelopeOscillatorTable * > * oscillatorStoreIonisation
void SetIonisationEnergy(G4double ie)
TFile * file
void BuildOscillatorTable(const G4Material *)
G4GLOB_DLL std::ostream G4cout
std::map< const G4Material *, G4double > * plasmaSquared
static constexpr double pi
Definition: G4SIunits.hh:75
std::map< const G4Material *, G4double > * atomsPerMolecule
TFile ff[ntarg]
Definition: Style.C:26
std::map< const G4Material *, G4double > * excitationEnergy
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:210
G4double GetTotalZ(const G4Material *)
static G4AtomicTransitionManager * Instance()
static constexpr double fine_structure_const
G4double bindingEnergy(G4int A, G4int Z)
G4double GetNumberOfZAtomsPerMolecule(const G4Material *, G4int Z)
static constexpr double Bohr_radius
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double BindingEnergy() const