ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeOscillatorManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 
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  //Copy helper in the oscillatorTable for Ionisation
807  //
808  //Oscillator table Ionisation for the material
809  G4PenelopeOscillatorTable* theTable = new G4PenelopeOscillatorTable(); //vector of oscillator
810  G4PenelopeOscillatorResEnergyComparator comparator;
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 
1067 void G4PenelopeOscillatorManager::ReadElementData()
1068 {
1069  if (verbosityLevel > 0)
1070  {
1071  G4cout << "G4PenelopeOscillatorManager::ReadElementData()" << G4endl;
1072  G4cout << "Going to read Element Data" << G4endl;
1073  }
1074  char* path = std::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 =
1094  G4AtomicTransitionManager::Instance();
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....
1154 G4double G4PenelopeOscillatorManager::GetMeanExcitationEnergy(const G4Material* mat)
1155 {
1156  // (1) First time, create oscillatorStores and read data
1157  CheckForTablesCreated();
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....
1178 G4double G4PenelopeOscillatorManager::GetPlasmaEnergySquared(const G4Material* mat)
1179 {
1180  // (1) First time, create oscillatorStores and read data
1181  CheckForTablesCreated();
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 
1203 G4double G4PenelopeOscillatorManager::GetAtomsPerMolecule(const G4Material* mat)
1204 {
1205  // (1) First time, create oscillatorStores and read data
1206  CheckForTablesCreated();
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 
1229 G4double G4PenelopeOscillatorManager::GetNumberOfZAtomsPerMolecule(const G4Material* mat,G4int Z)
1230 {
1231  // (1) First time, create oscillatorStores and read data
1232  CheckForTablesCreated();
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 }