Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4PenelopeOscillatorManager.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 // 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 
50 #include "G4PenelopeOscillatorManager.hh"
51 
52 #include "globals.hh"
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4AtomicTransitionManager.hh"
56 #include "G4AtomicShell.hh"
57 #include "G4Material.hh"
58 #include "G4Exp.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 
62 G4PenelopeOscillatorManager::G4PenelopeOscillatorManager() :
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 
78 G4PenelopeOscillatorManager::~G4PenelopeOscillatorManager()
79 {
80  Clear();
81  delete instance;
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
86 G4ThreadLocal G4PenelopeOscillatorManager* G4PenelopeOscillatorManager::instance = 0;
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
90 G4PenelopeOscillatorManager* G4PenelopeOscillatorManager::GetOscillatorManager()
91 {
92  if (!instance)
93  instance = new G4PenelopeOscillatorManager();
94  return instance;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
99 void G4PenelopeOscillatorManager::Clear()
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  }
118  delete oscillatorStoreIonisation;
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;
138  if (excitationEnergy) delete excitationEnergy;
139  if (plasmaSquared) delete plasmaSquared;
140  if (atomsPerMolecule) delete atomsPerMolecule;
141  if (atomTablePerMolecule) delete atomTablePerMolecule;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 
146 void G4PenelopeOscillatorManager::Dump(const G4Material* material)
147 {
148  G4PenelopeOscillatorTable* theTable = GetOscillatorTableIonisation(material);
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 
221 void G4PenelopeOscillatorManager::CheckForTablesCreated()
222 {
223  //Tables should be created at the same time, since they are both filled
224  //simultaneously
225  if (!oscillatorStoreIonisation)
226  {
227  oscillatorStoreIonisation = new std::map<const G4Material*,G4PenelopeOscillatorTable*>;
228  if (!fReadElementData)
229  ReadElementData();
230  if (!oscillatorStoreIonisation)
231  //It should be ok now
232  G4Exception("G4PenelopeOscillatorManager::GetOscillatorTableIonisation()",
233  "em2034",FatalException,
234  "Problem in allocating the Oscillator Store for Ionisation");
235  }
236 
237  if (!oscillatorStoreCompton)
238  {
239  oscillatorStoreCompton = new std::map<const G4Material*,G4PenelopeOscillatorTable*>;
240  if (!fReadElementData)
241  ReadElementData();
242  if (!oscillatorStoreCompton)
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>;
259  if (!atomTablePerMolecule)
260  atomTablePerMolecule = new std::map< std::pair<const G4Material*,G4int>, G4double>;
261 }
262 
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
265 
266 G4double G4PenelopeOscillatorManager::GetTotalZ(const G4Material* mat)
267 {
268  // (1) First time, create oscillatorStores and read data
269  CheckForTablesCreated();
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
276  BuildOscillatorTable(mat);
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 
291 G4double G4PenelopeOscillatorManager::GetTotalA(const G4Material* mat)
292 {
293  // (1) First time, create oscillatorStores and read data
294  CheckForTablesCreated();
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
301  BuildOscillatorTable(mat);
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 
316 G4PenelopeOscillatorTable* G4PenelopeOscillatorManager::GetOscillatorTableIonisation(const G4Material* mat)
317 {
318  // (1) First time, create oscillatorStores and read data
319  CheckForTablesCreated();
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
329  BuildOscillatorTable(mat);
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 
344 G4PenelopeOscillator* G4PenelopeOscillatorManager::GetOscillatorIonisation(const G4Material* material,
345  G4int index)
346 {
347  G4PenelopeOscillatorTable* theTable = GetOscillatorTableIonisation(material);
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 
363 G4PenelopeOscillatorTable* G4PenelopeOscillatorManager::GetOscillatorTableCompton(const G4Material* mat)
364 {
365  // (1) First time, create oscillatorStore and read data
366  CheckForTablesCreated();
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
376  BuildOscillatorTable(mat);
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 
391 G4PenelopeOscillator* G4PenelopeOscillatorManager::GetOscillatorCompton(const G4Material* material,
392  G4int index)
393 {
394  G4PenelopeOscillatorTable* theTable = GetOscillatorTableCompton(material);
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 
409 void G4PenelopeOscillatorManager::BuildOscillatorTable(const G4Material* material)
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  {
466  G4ExceptionDescription ed;
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]);
564  newOscLocal.SetHartreeFactor(elementData[4][i]/fine_structure_const);
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  {
668  G4ExceptionDescription ed;
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  {
758  G4ExceptionDescription ed;
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  {
767  G4ExceptionDescription ed;
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  //
807  //Copy helper in the oscillatorTable for Ionisation
808  //
809  //Oscillator table Ionisation for the material
810  G4PenelopeOscillatorTable* theTable = new G4PenelopeOscillatorTable(); //vector of oscillator
811  G4PenelopeOscillatorResEnergyComparator comparator;
812  std::sort(helper->begin(),helper->end(),comparator);
813 
814  //COPY THE HELPER (vector of object) to theTable (vector of Pointers).
815  for (size_t i=0;i<helper->size();i++)
816  {
817  //copy content --> one may need it later (e.g. to fill an other table, with variations)
818  G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
819  theTable->push_back(theOsc);
820  }
821 
822  //Oscillators of outer shells with resonance energies differing by a factor less than
823  //Rgroup are grouped as a single oscillator
824  G4double Rgroup = 1.05;
825  size_t Nost = theTable->size();
826 
827  size_t firstIndex = (isAConductor) ? 1 : 0; //for conductors, skip conduction oscillator
828  G4bool loopAgain = false;
829  G4int removedLevels = 0;
830  do
831  {
832  loopAgain = false;
833  if (Nost>firstIndex+1)
834  {
835  removedLevels = 0;
836  for (size_t i=firstIndex;i<theTable->size()-1;i++)
837  {
838  G4bool skipLoop = false;
839  G4int shellFlag = (*theTable)[i]->GetShellFlag();
840  G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
841  G4double resEne = (*theTable)[i]->GetResonanceEnergy();
842  G4double resEnePlus1 = (*theTable)[i+1]->GetResonanceEnergy();
843  G4double oscStre = (*theTable)[i]->GetOscillatorStrength();
844  G4double oscStrePlus1 = (*theTable)[i+1]->GetOscillatorStrength();
845  //if (shellFlag < 10 && ionEne>cutEnergy) in Penelope
846  if (ionEne>cutEnergy) //remove condition that shellFlag < 10!
847  skipLoop = true;
848  if (resEne<1.0*eV || resEnePlus1<1.0*eV)
849  skipLoop = true;
850  if (resEnePlus1 > Rgroup*resEne)
851  skipLoop = true;
852  if (!skipLoop)
853  {
854  G4double newRes = G4Exp((oscStre*std::log(resEne)+
855  oscStrePlus1*std::log(resEnePlus1))
856  /(oscStre+oscStrePlus1));
857  (*theTable)[i]->SetResonanceEnergy(newRes);
858  G4double newIon = (oscStre*ionEne+
859  oscStrePlus1*(*theTable)[i+1]->GetIonisationEnergy())/
860  (oscStre+oscStrePlus1);
861  (*theTable)[i]->SetIonisationEnergy(newIon);
862  G4double newStre = oscStre+oscStrePlus1;
863  (*theTable)[i]->SetOscillatorStrength(newStre);
864  G4double newHartree = (oscStre*(*theTable)[i]->GetHartreeFactor()+
865  oscStrePlus1*(*theTable)[i+1]->GetHartreeFactor())/
866  (oscStre+oscStrePlus1);
867  (*theTable)[i]->SetHartreeFactor(newHartree);
868  if ((*theTable)[i]->GetParentZ() != (*theTable)[i+1]->GetParentZ())
869  (*theTable)[i]->SetParentZ(0.);
870  if (shellFlag < 10 || (*theTable)[i+1]->GetShellFlag() < 10)
871  {
872  G4int newFlag = std::min(shellFlag,(*theTable)[i+1]->GetShellFlag());
873  (*theTable)[i]->SetShellFlag(newFlag);
874  }
875  else
876  (*theTable)[i]->SetShellFlag(30);
877  //We've lost anyway the track of the original level
878  (*theTable)[i]->SetParentShellID((*theTable)[i]->GetShellFlag());
879 
880 
881  if (i<theTable->size()-2)
882  {
883  for (size_t ii=i+1;ii<theTable->size()-1;ii++)
884  (*theTable)[ii] = (*theTable)[ii+1];
885  }
886  //G4cout << theTable->size() << G4endl;
887  theTable->erase(theTable->begin()+theTable->size()-1); //delete last element
888  removedLevels++;
889  }
890  }
891  }
892  if (removedLevels)
893  {
894  Nost -= removedLevels;
895  loopAgain = true;
896  }
897  if (Rgroup < 1.414213 || Nost > 64)
898  {
899  Rgroup = Rgroup*Rgroup;
900  loopAgain = true;
901  }
902  }while(loopAgain);
903 
904  if (verbosityLevel > 1)
905  {
906  G4cout << "Final grouping factor for Ionisation: " << Rgroup << G4endl;
907  }
908 
909  //Final Electron/Positron model parameters
910  for (size_t i=0;i<theTable->size();i++)
911  {
912  //Set cutoff recoil energy for the resonant mode
913  G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
914  if (ionEne < 1e-3*eV)
915  {
916  G4double resEne = (*theTable)[i]->GetResonanceEnergy();
917  (*theTable)[i]->SetIonisationEnergy(0.*eV);
918  (*theTable)[i]->SetCutoffRecoilResonantEnergy(resEne);
919  }
920  else
921  (*theTable)[i]->SetCutoffRecoilResonantEnergy(ionEne);
922  }
923 
924  //Last step
925  oscillatorStoreIonisation->insert(std::make_pair(material,theTable));
926 
927 
928  /*
929  SAME FOR COMPTON
930  */
931  //
932  //Copy helper in the oscillatorTable for Compton
933  //
934  //Oscillator table Ionisation for the material
935  G4PenelopeOscillatorTable* theTableC = new G4PenelopeOscillatorTable(); //vector of oscillator
936  //order by ionisation energy
937  std::sort(helper->begin(),helper->end());
938  //COPY THE HELPER (vector of object) to theTable (vector of Pointers).
939  for (size_t i=0;i<helper->size();i++)
940  {
941  //copy content --> one may need it later (e.g. to fill an other table, with variations)
942  G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
943  theTableC->push_back(theOsc);
944  }
945  //Oscillators of outer shells with resonance energies differing by a factor less than
946  //Rgroup are grouped as a single oscillator
947  Rgroup = 1.5;
948  Nost = theTableC->size();
949 
950  firstIndex = (isAConductor) ? 1 : 0; //for conductors, skip conduction oscillator
951  loopAgain = false;
952  removedLevels = 0;
953  do
954  {
955  loopAgain = false;
956  if (Nost>firstIndex+1)
957  {
958  removedLevels = 0;
959  for (size_t i=firstIndex;i<theTableC->size()-1;i++)
960  {
961  G4bool skipLoop = false;
962  //G4int shellFlag = (*theTableC)[i]->GetShellFlag();
963  G4double ionEne = (*theTableC)[i]->GetIonisationEnergy();
964  G4double ionEnePlus1 = (*theTableC)[i+1]->GetIonisationEnergy();
965  G4double oscStre = (*theTableC)[i]->GetOscillatorStrength();
966  G4double oscStrePlus1 = (*theTableC)[i+1]->GetOscillatorStrength();
967  //if (shellFlag < 10 && ionEne>cutEnergy) in Penelope
968  if (ionEne>cutEnergy)
969  skipLoop = true;
970  if (ionEne<1.0*eV || ionEnePlus1<1.0*eV)
971  skipLoop = true;
972  if (ionEnePlus1 > Rgroup*ionEne)
973  skipLoop = true;
974 
975  if (!skipLoop)
976  {
977  G4double newIon = (oscStre*ionEne+
978  oscStrePlus1*ionEnePlus1)/
979  (oscStre+oscStrePlus1);
980  (*theTableC)[i]->SetIonisationEnergy(newIon);
981  G4double newStre = oscStre+oscStrePlus1;
982  (*theTableC)[i]->SetOscillatorStrength(newStre);
983  G4double newHartree = (oscStre*(*theTableC)[i]->GetHartreeFactor()+
984  oscStrePlus1*(*theTableC)[i+1]->GetHartreeFactor())/
985  (oscStre+oscStrePlus1);
986  (*theTableC)[i]->SetHartreeFactor(newHartree);
987  if ((*theTableC)[i]->GetParentZ() != (*theTableC)[i+1]->GetParentZ())
988  (*theTableC)[i]->SetParentZ(0.);
989  (*theTableC)[i]->SetShellFlag(30);
990  (*theTableC)[i]->SetParentShellID((*theTableC)[i]->GetShellFlag());
991 
992  if (i<theTableC->size()-2)
993  {
994  for (size_t ii=i+1;ii<theTableC->size()-1;ii++)
995  (*theTableC)[ii] = (*theTableC)[ii+1];
996  }
997  theTableC->erase(theTableC->begin()+theTableC->size()-1); //delete last element
998  removedLevels++;
999  }
1000  }
1001  }
1002  if (removedLevels)
1003  {
1004  Nost -= removedLevels;
1005  loopAgain = true;
1006  }
1007  if (Rgroup < 2.0 || Nost > 64)
1008  {
1009  Rgroup = Rgroup*Rgroup;
1010  loopAgain = true;
1011  }
1012  }while(loopAgain);
1013 
1014 
1015  if (verbosityLevel > 1)
1016  {
1017  G4cout << "Final grouping factor for Compton: " << Rgroup << G4endl;
1018  }
1019 
1020  //Last step
1021  oscillatorStoreCompton->insert(std::make_pair(material,theTableC));
1022 
1023  /* //TESTING PURPOSES
1024  if (verbosityLevel > 1)
1025  {
1026  G4cout << "The table contains " << helper->size() << " oscillators " << G4endl;
1027  for (size_t k=0;k<helper->size();k++)
1028  {
1029  G4cout << "Oscillator # " << k << G4endl;
1030  G4cout << "Z = " << (*helper)[k].GetParentZ() << G4endl;
1031  G4cout << "Shell Flag = " << (*helper)[k].GetShellFlag() << G4endl;
1032  G4cout << "Compton index = " << (*helper)[k].GetHartreeFactor() << G4endl;
1033  G4cout << "Ionisation energy = " << (*helper)[k].GetIonisationEnergy()/eV << " eV" << G4endl;
1034  G4cout << "Occupation number = " << (*helper)[k].GetOscillatorStrength() << G4endl;
1035  G4cout << "Resonance energy = " << (*helper)[k].GetResonanceEnergy()/eV << " eV" << G4endl;
1036  }
1037 
1038  for (size_t k=0;k<helper->size();k++)
1039  {
1040  G4cout << k << " " << (*helper)[k].GetOscillatorStrength() << " " <<
1041  (*helper)[k].GetIonisationEnergy()/eV << " " << (*helper)[k].GetResonanceEnergy()/eV << " " <<
1042  (*helper)[k].GetParentZ() << " " << (*helper)[k].GetShellFlag() << " " <<
1043  (*helper)[k].GetHartreeFactor() << G4endl;
1044  }
1045  }
1046  */
1047 
1048 
1049  //CLEAN UP theHelper and its content
1050  delete helper;
1051  if (verbosityLevel > 1)
1052  Dump(material);
1053 
1054  return;
1055 }
1056 
1057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1058 
1059 void G4PenelopeOscillatorManager::ReadElementData()
1060 {
1061  if (verbosityLevel > 0)
1062  {
1063  G4cout << "G4PenelopeOscillatorManager::ReadElementData()" << G4endl;
1064  G4cout << "Going to read Element Data" << G4endl;
1065  }
1066  char* path = getenv("G4LEDATA");
1067  if (!path)
1068  {
1069  G4String excep = "G4PenelopeOscillatorManager - G4LEDATA environment variable not set!";
1070  G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
1071  "em0006",FatalException,excep);
1072  return;
1073  }
1074  G4String pathString(path);
1075  G4String pathFile = pathString + "/penelope/pdatconf.p08";
1076  std::ifstream file(pathFile);
1077 
1078  if (!file.is_open())
1079  {
1080  G4String excep = "G4PenelopeOscillatorManager - data file " + pathFile + " not found!";
1081  G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
1082  "em0003",FatalException,excep);
1083  }
1084 
1085  G4AtomicTransitionManager* theTransitionManager =
1086  G4AtomicTransitionManager::Instance();
1087  theTransitionManager->Initialise();
1088 
1089 
1090  //Read header (22 lines)
1091  G4String theHeader;
1092  for (G4int iline=0;iline<22;iline++)
1093  getline(file,theHeader);
1094  //Done
1095  G4int Z=0;
1096  G4int shellCode = 0;
1097  G4String shellId = "NULL";
1098  G4int occupationNumber = 0;
1099  G4double ionisationEnergy = 0.0*eV;
1100  G4double hartreeProfile = 0.;
1101  G4int shellCounter = 0;
1102  G4int oldZ = -1;
1103  G4int numberOfShells = 0;
1104  //Start reading data
1105  for (G4int i=0;!file.eof();i++)
1106  {
1107  file >> Z >> shellCode >> shellId >> occupationNumber >> ionisationEnergy >> hartreeProfile;
1108  if (Z>0 && i<2000)
1109  {
1110  elementData[0][i] = Z;
1111  elementData[1][i] = shellCode;
1112  elementData[2][i] = occupationNumber;
1113  //reset things
1114  if (Z != oldZ)
1115  {
1116  shellCounter = 0;
1117  oldZ = Z;
1118  numberOfShells = theTransitionManager->NumberOfShells(Z);
1119  }
1120  G4double bindingEnergy = -1*eV;
1121  if (shellCounter<numberOfShells)
1122  {
1123  G4AtomicShell* shell = theTransitionManager->Shell(Z,shellCounter);
1124  bindingEnergy = shell->BindingEnergy();
1125  }
1126  //Valid level found in the G4AtomicTransition database: keep it, otherwise use
1127  //the ionisation energy found in the Penelope database
1128  elementData[3][i] = (bindingEnergy>100*eV) ? bindingEnergy : ionisationEnergy*eV;
1129  //elementData[3][i] = ionisationEnergy*eV;
1130  elementData[4][i] = hartreeProfile;
1131  shellCounter++;
1132  }
1133  }
1134  file.close();
1135 
1136  if (verbosityLevel > 1)
1137  {
1138  G4cout << "G4PenelopeOscillatorManager::ReadElementData(): Data file read" << G4endl;
1139  }
1140  fReadElementData = true;
1141  return;
1142 
1143 }
1144 
1145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1146 G4double G4PenelopeOscillatorManager::GetMeanExcitationEnergy(const G4Material* mat)
1147 {
1148  // (1) First time, create oscillatorStores and read data
1149  CheckForTablesCreated();
1150 
1151  // (2) Check if the material has been already included
1152  if (excitationEnergy->count(mat))
1153  return excitationEnergy->find(mat)->second;
1154 
1155  // (3) If we are here, it means that we have to create the table for the material
1156  BuildOscillatorTable(mat);
1157 
1158  // (4) now, the oscillator store should be ok
1159  if (excitationEnergy->count(mat))
1160  return excitationEnergy->find(mat)->second;
1161  else
1162  {
1163  G4cout << "G4PenelopeOscillatorManager::GetMolecularExcitationEnergy() " << G4endl;
1164  G4cout << "Impossible to retrieve the excitation energy for " << mat->GetName() << G4endl;
1165  return 0;
1166  }
1167 }
1168 
1169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1170 G4double G4PenelopeOscillatorManager::GetPlasmaEnergySquared(const G4Material* mat)
1171 {
1172  // (1) First time, create oscillatorStores and read data
1173  CheckForTablesCreated();
1174 
1175  // (2) Check if the material has been already included
1176  if (plasmaSquared->count(mat))
1177  return plasmaSquared->find(mat)->second;
1178 
1179  // (3) If we are here, it means that we have to create the table for the material
1180  BuildOscillatorTable(mat);
1181 
1182  // (4) now, the oscillator store should be ok
1183  if (plasmaSquared->count(mat))
1184  return plasmaSquared->find(mat)->second;
1185  else
1186  {
1187  G4cout << "G4PenelopeOscillatorManager::GetPlasmaEnergySquared() " << G4endl;
1188  G4cout << "Impossible to retrieve the plasma energy for " << mat->GetName() << G4endl;
1189  return 0;
1190  }
1191 }
1192 
1193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1194 
1195 G4double G4PenelopeOscillatorManager::GetAtomsPerMolecule(const G4Material* mat)
1196 {
1197  // (1) First time, create oscillatorStores and read data
1198  CheckForTablesCreated();
1199 
1200  // (2) Check if the material has been already included
1201  if (atomsPerMolecule->count(mat))
1202  return atomsPerMolecule->find(mat)->second;
1203 
1204  // (3) If we are here, it means that we have to create the table for the material
1205  BuildOscillatorTable(mat);
1206 
1207  // (4) now, the oscillator store should be ok
1208  if (atomsPerMolecule->count(mat))
1209  return atomsPerMolecule->find(mat)->second;
1210  else
1211  {
1212  G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
1213  G4cout << "Impossible to retrieve the number of atoms per molecule for "
1214  << mat->GetName() << G4endl;
1215  return 0;
1216  }
1217 }
1218 
1219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1220 
1221 G4double G4PenelopeOscillatorManager::GetNumberOfZAtomsPerMolecule(const G4Material* mat,G4int Z)
1222 {
1223  // (1) First time, create oscillatorStores and read data
1224  CheckForTablesCreated();
1225 
1226  // (2) Check if the material/Z couple has been already included
1227  std::pair<const G4Material*,G4int> theKey = std::make_pair(mat,Z);
1228  if (atomTablePerMolecule->count(theKey))
1229  return atomTablePerMolecule->find(theKey)->second;
1230 
1231  // (3) If we are here, it means that we have to create the table for the material
1232  BuildOscillatorTable(mat);
1233 
1234  // (4) now, the oscillator store should be ok
1235  if (atomTablePerMolecule->count(theKey))
1236  return atomTablePerMolecule->find(theKey)->second;
1237  else
1238  {
1239  G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
1240  G4cout << "Impossible to retrieve the number of atoms per molecule for Z = "
1241  << Z << " in material " << mat->GetName() << G4endl;
1242  return 0;
1243  }
1244 }
G4Exp.hh
G4PenelopeOscillatorManager::GetOscillatorTableIonisation
G4PenelopeOscillatorTable * GetOscillatorTableIonisation(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:316
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4PenelopeOscillatorManager::GetPlasmaEnergySquared
G4double GetPlasmaEnergySquared(const G4Material *)
Returns the squared plasma energy.
Definition: G4PenelopeOscillatorManager.cc:1170
CLHEP::Bohr_radius
static constexpr double Bohr_radius
Definition: PhysicalConstants.h:99
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:178
G4AtomicShell.hh
G4PenelopeOscillatorManager.hh
G4AtomicShell::BindingEnergy
G4double BindingEnergy() const
Definition: G4AtomicShell.hh:80
G4PenelopeOscillatorManager::GetOscillatorIonisation
G4PenelopeOscillator * GetOscillatorIonisation(const G4Material *, G4int)
Definition: G4PenelopeOscillatorManager.cc:344
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:89
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
G4int
int G4int
Definition: G4Types.hh:78
G4PenelopeOscillatorManager::GetOscillatorTableCompton
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:363
G4PenelopeOscillator::SetParentZ
void SetParentZ(G4double parZ)
Definition: G4PenelopeOscillator.hh:79
G4PenelopeOscillator::SetOscillatorStrength
void SetOscillatorStrength(G4double ostr)
Definition: G4PenelopeOscillator.hh:73
G4cout
G4GLOB_DLL std::ostream G4cout
G4PenelopeOscillator::SetHartreeFactor
void SetHartreeFactor(G4double hf)
Definition: G4PenelopeOscillator.hh:64
G4PenelopeOscillator::SetIonisationEnergy
void SetIonisationEnergy(G4double ie)
Definition: G4PenelopeOscillator.hh:67
G4bool
bool G4bool
Definition: G4Types.hh:79
G4Material.hh
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4PenelopeOscillatorManager::GetOscillatorManager
static G4PenelopeOscillatorManager * GetOscillatorManager()
Definition: G4PenelopeOscillatorManager.cc:90
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
G4PenelopeOscillator::SetParentShellID
void SetParentShellID(G4int psID)
Definition: G4PenelopeOscillator.hh:82
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
G4AtomicTransitionManager.hh
G4PenelopeOscillatorManager::GetTotalA
G4double GetTotalA(const G4Material *)
Returns the total A for the molecule.
Definition: G4PenelopeOscillatorManager.cc:291
G4PenelopeOscillatorTable
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
Definition: G4PenelopeOscillatorManager.hh:63
FatalException
Definition: G4ExceptionSeverity.hh:60
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4endl
#define G4endl
Definition: G4ios.hh:61
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4PenelopeOscillatorManager::GetOscillatorCompton
G4PenelopeOscillator * GetOscillatorCompton(const G4Material *, G4int)
Definition: G4PenelopeOscillatorManager.cc:391
CLHEP::fine_structure_const
static constexpr double fine_structure_const
Definition: PhysicalConstants.h:96
G4PenelopeOscillatorManager::GetNumberOfZAtomsPerMolecule
G4double GetNumberOfZAtomsPerMolecule(const G4Material *, G4int Z)
Definition: G4PenelopeOscillatorManager.cc:1221
G4PenelopeOscillatorManager::GetAtomsPerMolecule
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
Definition: G4PenelopeOscillatorManager.cc:1195
G4PenelopeOscillatorManager::GetMeanExcitationEnergy
G4double GetMeanExcitationEnergy(const G4Material *)
Returns the mean excitation energy.
Definition: G4PenelopeOscillatorManager.cc:1146
G4Material::GetFractionVector
const G4double * GetFractionVector() const
Definition: G4Material.hh:194
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:37
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:45
G4AtomicTransitionManager::Shell
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
Definition: G4AtomicTransitionManager.cc:95
G4PenelopeOscillatorManager::GetTotalZ
G4double GetTotalZ(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:266
G4PenelopeOscillator::SetShellFlag
void SetShellFlag(G4int theflag)
Definition: G4PenelopeOscillator.hh:76
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262