Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60 
61 G4PenelopeOscillatorManager::G4PenelopeOscillatorManager() :
62  oscillatorStoreIonisation(0),oscillatorStoreCompton(0),atomicNumber(0),
63  atomicMass(0),excitationEnergy(0),plasmaSquared(0),atomsPerMolecule(0),
64  atomTablePerMolecule(0)
65 {
66  fReadElementData = false;
67  for (G4int i=0;i<5;i++)
68  {
69  for (G4int j=0;j<2000;j++)
70  elementData[i][j] = 0.;
71  }
72  verbosityLevel = 0;
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
77 G4PenelopeOscillatorManager::~G4PenelopeOscillatorManager()
78 {
79  Clear();
80  delete instance;
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
85 G4PenelopeOscillatorManager* G4PenelopeOscillatorManager::instance = 0;
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
89 G4PenelopeOscillatorManager* G4PenelopeOscillatorManager::GetOscillatorManager()
90 {
91  if (!instance)
92  instance = new G4PenelopeOscillatorManager();
93  return instance;
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
98 void G4PenelopeOscillatorManager::Clear()
99 {
100  if (verbosityLevel > 1)
101  G4cout << " G4PenelopeOscillatorManager::Clear() - Clean Oscillator Tables" << G4endl;
102 
103  //Clean up OscillatorStoreIonisation
104  std::map<const G4Material*,G4PenelopeOscillatorTable*>::iterator i;
105  for (i=oscillatorStoreIonisation->begin();i != oscillatorStoreIonisation->end();i++)
106  {
107  G4PenelopeOscillatorTable* table = i->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 (i=oscillatorStoreCompton->begin();i != oscillatorStoreCompton->end();i++)
122  {
123  G4PenelopeOscillatorTable* table = i->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 NULL;
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 NULL;
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 NULL;
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 NULL;
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 = std::exp(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  if (std::fabs(occup) > 0)
559  {
560  G4PenelopeOscillator newOscLocal;
561  newOscLocal.SetOscillatorStrength(std::fabs(occup)*(*StechiometricFactors)[k]);
562  newOscLocal.SetIonisationEnergy(elementData[3][i]);
563  newOscLocal.SetHartreeFactor(elementData[4][i]/fine_structure_const);
564  newOscLocal.SetParentZ(elementData[0][i]);
565  //keep track of the origianl shell level
566  newOscLocal.SetParentShellID((G4int)elementData[1][i]);
567  //register only K, L and M shells. Outer shells all grouped with
568  //shellIndex = 30
569  if (elementData[0][i] > 6 && elementData[1][i] < 10)
570  newOscLocal.SetShellFlag(((G4int)elementData[1][i]));
571  else
572  newOscLocal.SetShellFlag(30);
573  helper->push_back(newOscLocal);
574  if (occup < 0)
575  {
576  G4double ff = (*helper)[0].GetOscillatorStrength();
577  ff += std::fabs(occup)*(*StechiometricFactors)[k];
578  (*helper)[0].SetOscillatorStrength(ff);
579  }
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  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.5 || 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 
699  //Sternheimer's adjustment factor
700  G4double adjustmentFactor = 0;
701  if (helper->size() > 1)
702  {
703  G4double TST = totalZ*std::log(meanExcitationEnergy/eV);
704  G4double AALow = 0.5;
705  G4double AAHigh = 10.;
706  do
707  {
708  adjustmentFactor = (AALow+AAHigh)*0.5;
709  G4double sumLocal = 0;
710  for (size_t i=0;i<helper->size();i++)
711  {
712  if (i == 0 && isAConductor)
713  {
714  G4double resEne = (*helper)[i].GetResonanceEnergy();
715  sumLocal += (*helper)[i].GetOscillatorStrength()*std::log(resEne/eV);
716  }
717  else
718  {
719  G4double ionEne = (*helper)[i].GetIonisationEnergy();
720  G4double oscStre = (*helper)[i].GetOscillatorStrength();
721  G4double WI2 = (adjustmentFactor*adjustmentFactor*ionEne*ionEne) +
722  2./3.*(oscStre/totalZ)*Omega*Omega;
723  G4double resEne = std::sqrt(WI2);
724  (*helper)[i].SetResonanceEnergy(resEne);
725  sumLocal += (*helper)[i].GetOscillatorStrength()*std::log(resEne/eV);
726  }
727  }
728  if (sumLocal < TST)
729  AALow = adjustmentFactor;
730  else
731  AAHigh = adjustmentFactor;
732  }while((AAHigh-AALow)>(1e-14*adjustmentFactor));
733  }
734  else
735  {
736  G4double ionEne = (*helper)[0].GetIonisationEnergy();
737  (*helper)[0].SetIonisationEnergy(std::fabs(ionEne));
738  (*helper)[0].SetResonanceEnergy(meanExcitationEnergy);
739  }
740  if (verbosityLevel > 1)
741  {
742  G4cout << "Sternheimer's adjustment factor: " << adjustmentFactor << G4endl;
743  }
744 
745  //Check again for data consistency
746  G4double xcheck = (*helper)[0].GetOscillatorStrength()*std::log((*helper)[0].GetResonanceEnergy());
747  G4double TST = (*helper)[0].GetOscillatorStrength();
748  for (size_t i=1;i<helper->size();i++)
749  {
750  xcheck += (*helper)[i].GetOscillatorStrength()*std::log((*helper)[i].GetResonanceEnergy());
751  TST += (*helper)[i].GetOscillatorStrength();
752  }
753  if (std::fabs(TST-totalZ)>1e-8*totalZ)
754  {
755  G4ExceptionDescription ed;
756  ed << "Inconsistent oscillator data " << G4endl;
757  ed << TST << " " << totalZ << G4endl;
758  G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
759  "em2036",FatalException,ed);
760  }
761  xcheck = std::exp(xcheck/totalZ);
762  if (std::fabs(xcheck-meanExcitationEnergy) > 1e-8*meanExcitationEnergy)
763  {
764  G4ExceptionDescription ed;
765  ed << "Error in Sterheimer factor calculation " << G4endl;
766  ed << xcheck/eV << " " << meanExcitationEnergy/eV << G4endl;
767  G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
768  "em2037",FatalException,ed);
769  }
770 
771  //Selection of the lowest ionisation energy for inner shells. Only the K, L and M shells with
772  //ionisation energy less than the N1 shell of the heaviest element in the material are considered as
773  //inner shells. As a results, the inner/outer shell character of an atomic shell depends on the
774  //composition of the material.
775  G4double Zmax = 0;
776  for (G4int k=0;k<nElements;k++)
777  {
778  G4double Z = (*elementVector)[k]->GetZ();
779  if (Z>Zmax) Zmax = Z;
780  }
781  //Find N1 level of the heaviest element (if any).
782  G4bool found = false;
783  G4double cutEnergy = 50*eV;
784  for (size_t i=0;i<helper->size() && !found;i++)
785  {
786  G4double Z = (*helper)[i].GetParentZ();
787  G4int shID = (*helper)[i].GetParentShellID(); //look for the N1 level
788  if (shID == 10 && Z == Zmax)
789  {
790  found = true;
791  if ((*helper)[i].GetIonisationEnergy() > cutEnergy)
792  cutEnergy = (*helper)[i].GetIonisationEnergy();
793  }
794  }
795  //Make that cutEnergy cannot be higher than 250 eV, namely the fluorescence level by
796  //Geant4
797  G4double lowEnergyLimitForFluorescence = 250*eV;
798  cutEnergy = std::min(cutEnergy,lowEnergyLimitForFluorescence);
799 
800  if (verbosityLevel > 1)
801  G4cout << "Cutoff energy: " << cutEnergy/eV << " eV" << G4endl;
802 
803  //
804  //Copy helper in the oscillatorTable for Ionisation
805  //
806  //Oscillator table Ionisation for the material
807  G4PenelopeOscillatorTable* theTable = new G4PenelopeOscillatorTable(); //vector of oscillator
808  G4PenelopeOscillatorResEnergyComparator comparator;
809  std::sort(helper->begin(),helper->end(),comparator);
810 
811  //COPY THE HELPER (vector of object) to theTable (vector of Pointers).
812  for (size_t i=0;i<helper->size();i++)
813  {
814  //copy content --> one may need it later (e.g. to fill an other table, with variations)
815  G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
816  theTable->push_back(theOsc);
817  }
818 
819  //Oscillators of outer shells with resonance energies differing by a factor less than
820  //Rgroup are grouped as a single oscillator
821  G4double Rgroup = 1.05;
822  size_t Nost = theTable->size();
823 
824  size_t firstIndex = (isAConductor) ? 1 : 0; //for conductors, skip conduction oscillator
825  G4bool loopAgain = false;
826  G4int removedLevels = 0;
827  do
828  {
829  loopAgain = false;
830  if (Nost>firstIndex+1)
831  {
832  removedLevels = 0;
833  for (size_t i=firstIndex;i<theTable->size()-1;i++)
834  {
835  G4bool skipLoop = false;
836  G4int shellFlag = (*theTable)[i]->GetShellFlag();
837  G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
838  G4double resEne = (*theTable)[i]->GetResonanceEnergy();
839  G4double resEnePlus1 = (*theTable)[i+1]->GetResonanceEnergy();
840  G4double oscStre = (*theTable)[i]->GetOscillatorStrength();
841  G4double oscStrePlus1 = (*theTable)[i+1]->GetOscillatorStrength();
842  //if (shellFlag < 10 && ionEne>cutEnergy) in Penelope
843  if (ionEne>cutEnergy) //remove condition that shellFlag < 10!
844  skipLoop = true;
845  if (resEne<1.0*eV || resEnePlus1<1.0*eV)
846  skipLoop = true;
847  if (resEnePlus1 > Rgroup*resEne)
848  skipLoop = true;
849  if (!skipLoop)
850  {
851  G4double newRes = std::exp((oscStre*std::log(resEne)+
852  oscStrePlus1*std::log(resEnePlus1))
853  /(oscStre+oscStrePlus1));
854  (*theTable)[i]->SetResonanceEnergy(newRes);
855  G4double newIon = (oscStre*ionEne+
856  oscStrePlus1*(*theTable)[i+1]->GetIonisationEnergy())/
857  (oscStre+oscStrePlus1);
858  (*theTable)[i]->SetIonisationEnergy(newIon);
859  G4double newStre = oscStre+oscStrePlus1;
860  (*theTable)[i]->SetOscillatorStrength(newStre);
861  G4double newHartree = (oscStre*(*theTable)[i]->GetHartreeFactor()+
862  oscStrePlus1*(*theTable)[i+1]->GetHartreeFactor())/
863  (oscStre+oscStrePlus1);
864  (*theTable)[i]->SetHartreeFactor(newHartree);
865  if ((*theTable)[i]->GetParentZ() != (*theTable)[i+1]->GetParentZ())
866  (*theTable)[i]->SetParentZ(0.);
867  if (shellFlag < 10 || (*theTable)[i+1]->GetShellFlag() < 10)
868  {
869  G4int newFlag = std::min(shellFlag,(*theTable)[i+1]->GetShellFlag());
870  (*theTable)[i]->SetShellFlag(newFlag);
871  }
872  else
873  (*theTable)[i]->SetShellFlag(30);
874  //We've lost anyway the track of the original level
875  (*theTable)[i]->SetParentShellID((*theTable)[i]->GetShellFlag());
876 
877 
878  if (i<theTable->size()-2)
879  {
880  for (size_t ii=i+1;ii<theTable->size()-1;ii++)
881  (*theTable)[ii] = (*theTable)[ii+1];
882  }
883  //G4cout << theTable->size() << G4endl;
884  theTable->erase(theTable->begin()+theTable->size()-1); //delete last element
885  removedLevels++;
886  }
887  }
888  }
889  if (removedLevels)
890  {
891  Nost -= removedLevels;
892  loopAgain = true;
893  }
894  if (Rgroup < 1.414213 || Nost > 64)
895  {
896  Rgroup = Rgroup*Rgroup;
897  loopAgain = true;
898  }
899  }while(loopAgain);
900 
901  if (verbosityLevel > 1)
902  {
903  G4cout << "Final grouping factor for Ionisation: " << Rgroup << G4endl;
904  }
905 
906  //Final Electron/Positron model parameters
907  for (size_t i=0;i<theTable->size();i++)
908  {
909  //Set cutoff recoil energy for the resonant mode
910  G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
911  if (ionEne < 1e-3*eV)
912  {
913  G4double resEne = (*theTable)[i]->GetResonanceEnergy();
914  (*theTable)[i]->SetIonisationEnergy(0.*eV);
915  (*theTable)[i]->SetCutoffRecoilResonantEnergy(resEne);
916  }
917  else
918  (*theTable)[i]->SetCutoffRecoilResonantEnergy(ionEne);
919  }
920 
921  //Last step
922  oscillatorStoreIonisation->insert(std::make_pair(material,theTable));
923 
924 
925  /*
926  SAME FOR COMPTON
927  */
928  //
929  //Copy helper in the oscillatorTable for Compton
930  //
931  //Oscillator table Ionisation for the material
932  G4PenelopeOscillatorTable* theTableC = new G4PenelopeOscillatorTable(); //vector of oscillator
933  //order by ionisation energy
934  std::sort(helper->begin(),helper->end());
935  //COPY THE HELPER (vector of object) to theTable (vector of Pointers).
936  for (size_t i=0;i<helper->size();i++)
937  {
938  //copy content --> one may need it later (e.g. to fill an other table, with variations)
939  G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
940  theTableC->push_back(theOsc);
941  }
942  //Oscillators of outer shells with resonance energies differing by a factor less than
943  //Rgroup are grouped as a single oscillator
944  Rgroup = 1.5;
945  Nost = theTableC->size();
946 
947  firstIndex = (isAConductor) ? 1 : 0; //for conductors, skip conduction oscillator
948  loopAgain = false;
949  removedLevels = 0;
950  do
951  {
952  loopAgain = false;
953  if (Nost>firstIndex+1)
954  {
955  removedLevels = 0;
956  for (size_t i=firstIndex;i<theTableC->size()-1;i++)
957  {
958  G4bool skipLoop = false;
959  //G4int shellFlag = (*theTableC)[i]->GetShellFlag();
960  G4double ionEne = (*theTableC)[i]->GetIonisationEnergy();
961  G4double ionEnePlus1 = (*theTableC)[i+1]->GetIonisationEnergy();
962  G4double oscStre = (*theTableC)[i]->GetOscillatorStrength();
963  G4double oscStrePlus1 = (*theTableC)[i+1]->GetOscillatorStrength();
964  //if (shellFlag < 10 && ionEne>cutEnergy) in Penelope
965  if (ionEne>cutEnergy)
966  skipLoop = true;
967  if (ionEne<1.0*eV || ionEnePlus1<1.0*eV)
968  skipLoop = true;
969  if (ionEnePlus1 > Rgroup*ionEne)
970  skipLoop = true;
971 
972  if (!skipLoop)
973  {
974  G4double newIon = (oscStre*ionEne+
975  oscStrePlus1*ionEnePlus1)/
976  (oscStre+oscStrePlus1);
977  (*theTableC)[i]->SetIonisationEnergy(newIon);
978  G4double newStre = oscStre+oscStrePlus1;
979  (*theTableC)[i]->SetOscillatorStrength(newStre);
980  G4double newHartree = (oscStre*(*theTableC)[i]->GetHartreeFactor()+
981  oscStrePlus1*(*theTableC)[i+1]->GetHartreeFactor())/
982  (oscStre+oscStrePlus1);
983  (*theTableC)[i]->SetHartreeFactor(newHartree);
984  if ((*theTableC)[i]->GetParentZ() != (*theTableC)[i+1]->GetParentZ())
985  (*theTableC)[i]->SetParentZ(0.);
986  (*theTableC)[i]->SetShellFlag(30);
987  (*theTableC)[i]->SetParentShellID((*theTableC)[i]->GetShellFlag());
988 
989  if (i<theTableC->size()-2)
990  {
991  for (size_t ii=i+1;ii<theTableC->size()-1;ii++)
992  (*theTableC)[ii] = (*theTableC)[ii+1];
993  }
994  theTableC->erase(theTableC->begin()+theTableC->size()-1); //delete last element
995  removedLevels++;
996  }
997  }
998  }
999  if (removedLevels)
1000  {
1001  Nost -= removedLevels;
1002  loopAgain = true;
1003  }
1004  if (Rgroup < 2.0 || Nost > 64)
1005  {
1006  Rgroup = Rgroup*Rgroup;
1007  loopAgain = true;
1008  }
1009  }while(loopAgain);
1010 
1011 
1012  if (verbosityLevel > 1)
1013  {
1014  G4cout << "Final grouping factor for Compton: " << Rgroup << G4endl;
1015  }
1016 
1017  //Last step
1018  oscillatorStoreCompton->insert(std::make_pair(material,theTableC));
1019 
1020  /* //TESTING PURPOSES
1021  if (verbosityLevel > 1)
1022  {
1023  G4cout << "The table contains " << helper->size() << " oscillators " << G4endl;
1024  for (size_t k=0;k<helper->size();k++)
1025  {
1026  G4cout << "Oscillator # " << k << G4endl;
1027  G4cout << "Z = " << (*helper)[k].GetParentZ() << G4endl;
1028  G4cout << "Shell Flag = " << (*helper)[k].GetShellFlag() << G4endl;
1029  G4cout << "Compton index = " << (*helper)[k].GetHartreeFactor() << G4endl;
1030  G4cout << "Ionisation energy = " << (*helper)[k].GetIonisationEnergy()/eV << " eV" << G4endl;
1031  G4cout << "Occupation number = " << (*helper)[k].GetOscillatorStrength() << G4endl;
1032  G4cout << "Resonance energy = " << (*helper)[k].GetResonanceEnergy()/eV << " eV" << G4endl;
1033  }
1034 
1035  for (size_t k=0;k<helper->size();k++)
1036  {
1037  G4cout << k << " " << (*helper)[k].GetOscillatorStrength() << " " <<
1038  (*helper)[k].GetIonisationEnergy()/eV << " " << (*helper)[k].GetResonanceEnergy()/eV << " " <<
1039  (*helper)[k].GetParentZ() << " " << (*helper)[k].GetShellFlag() << " " <<
1040  (*helper)[k].GetHartreeFactor() << G4endl;
1041  }
1042  }
1043  */
1044 
1045 
1046  //CLEAN UP theHelper and its content
1047  delete helper;
1048  if (verbosityLevel > 1)
1049  Dump(material);
1050 
1051  return;
1052 }
1053 
1054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1055 
1056 void G4PenelopeOscillatorManager::ReadElementData()
1057 {
1058  if (verbosityLevel > 0)
1059  {
1060  G4cout << "G4PenelopeOscillatorManager::ReadElementData()" << G4endl;
1061  G4cout << "Going to read Element Data" << G4endl;
1062  }
1063  char* path = getenv("G4LEDATA");
1064  if (!path)
1065  {
1066  G4String excep = "G4PenelopeOscillatorManager - G4LEDATA environment variable not set!";
1067  G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
1068  "em0006",FatalException,excep);
1069  return;
1070  }
1071  G4String pathString(path);
1072  G4String pathFile = pathString + "/penelope/pdatconf.p08";
1073  std::ifstream file(pathFile);
1074 
1075  if (!file.is_open())
1076  {
1077  G4String excep = "G4PenelopeOscillatorManager - data file " + pathFile + " not found!";
1078  G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
1079  "em0003",FatalException,excep);
1080  }
1081 
1082  G4AtomicTransitionManager* theTransitionManager =
1083  G4AtomicTransitionManager::Instance();
1084 
1085  //Read header (22 lines)
1086  G4String theHeader;
1087  for (G4int iline=0;iline<22;iline++)
1088  getline(file,theHeader);
1089  //Done
1090  G4int Z=0;
1091  G4int shellCode = 0;
1092  G4String shellId = "NULL";
1093  G4int occupationNumber = 0;
1094  G4double ionisationEnergy = 0.0*eV;
1095  G4double hartreeProfile = 0.;
1096  G4int shellCounter = 0;
1097  G4int oldZ = -1;
1098  G4int numberOfShells = 0;
1099  //Start reading data
1100  for (G4int i=0;!file.eof();i++)
1101  {
1102  file >> Z >> shellCode >> shellId >> occupationNumber >> ionisationEnergy >> hartreeProfile;
1103  if (Z>0 && i<2000)
1104  {
1105  elementData[0][i] = Z;
1106  elementData[1][i] = shellCode;
1107  elementData[2][i] = occupationNumber;
1108  //reset things
1109  if (Z != oldZ)
1110  {
1111  shellCounter = 0;
1112  oldZ = Z;
1113  numberOfShells = theTransitionManager->NumberOfShells(Z);
1114  }
1115  G4double bindingEnergy = -1*eV;
1116  if (shellCounter<numberOfShells)
1117  {
1118  G4AtomicShell* shell = theTransitionManager->Shell(Z,shellCounter);
1119  bindingEnergy = shell->BindingEnergy();
1120  }
1121  //Valid level found in the G4AtomicTransition database: keep it, otherwise use
1122  //the ionisation energy found in the Penelope database
1123  elementData[3][i] = (bindingEnergy>100*eV) ? bindingEnergy : ionisationEnergy*eV;
1124  //elementData[3][i] = ionisationEnergy*eV;
1125  elementData[4][i] = hartreeProfile;
1126  shellCounter++;
1127  }
1128  }
1129  file.close();
1130 
1131  if (verbosityLevel > 1)
1132  {
1133  G4cout << "G4PenelopeOscillatorManager::ReadElementData(): Data file read" << G4endl;
1134  }
1135  fReadElementData = true;
1136  return;
1137 
1138 }
1139 
1140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1141 G4double G4PenelopeOscillatorManager::GetMeanExcitationEnergy(const G4Material* mat)
1142 {
1143  // (1) First time, create oscillatorStores and read data
1144  CheckForTablesCreated();
1145 
1146  // (2) Check if the material has been already included
1147  if (excitationEnergy->count(mat))
1148  return excitationEnergy->find(mat)->second;
1149 
1150  // (3) If we are here, it means that we have to create the table for the material
1151  BuildOscillatorTable(mat);
1152 
1153  // (4) now, the oscillator store should be ok
1154  if (excitationEnergy->count(mat))
1155  return excitationEnergy->find(mat)->second;
1156  else
1157  {
1158  G4cout << "G4PenelopeOscillatorManager::GetMolecularExcitationEnergy() " << G4endl;
1159  G4cout << "Impossible to retrieve the excitation energy for " << mat->GetName() << G4endl;
1160  return 0;
1161  }
1162 }
1163 
1164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1165 G4double G4PenelopeOscillatorManager::GetPlasmaEnergySquared(const G4Material* mat)
1166 {
1167  // (1) First time, create oscillatorStores and read data
1168  CheckForTablesCreated();
1169 
1170  // (2) Check if the material has been already included
1171  if (plasmaSquared->count(mat))
1172  return plasmaSquared->find(mat)->second;
1173 
1174  // (3) If we are here, it means that we have to create the table for the material
1175  BuildOscillatorTable(mat);
1176 
1177  // (4) now, the oscillator store should be ok
1178  if (plasmaSquared->count(mat))
1179  return plasmaSquared->find(mat)->second;
1180  else
1181  {
1182  G4cout << "G4PenelopeOscillatorManager::GetPlasmaEnergySquared() " << G4endl;
1183  G4cout << "Impossible to retrieve the plasma energy for " << mat->GetName() << G4endl;
1184  return 0;
1185  }
1186 }
1187 
1188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1189 
1190 G4double G4PenelopeOscillatorManager::GetAtomsPerMolecule(const G4Material* mat)
1191 {
1192  // (1) First time, create oscillatorStores and read data
1193  CheckForTablesCreated();
1194 
1195  // (2) Check if the material has been already included
1196  if (atomsPerMolecule->count(mat))
1197  return atomsPerMolecule->find(mat)->second;
1198 
1199  // (3) If we are here, it means that we have to create the table for the material
1200  BuildOscillatorTable(mat);
1201 
1202  // (4) now, the oscillator store should be ok
1203  if (atomsPerMolecule->count(mat))
1204  return atomsPerMolecule->find(mat)->second;
1205  else
1206  {
1207  G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
1208  G4cout << "Impossible to retrieve the number of atoms per molecule for "
1209  << mat->GetName() << G4endl;
1210  return 0;
1211  }
1212 }
1213 
1214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1215 
1216 G4double G4PenelopeOscillatorManager::GetNumberOfZAtomsPerMolecule(const G4Material* mat,G4int Z)
1217 {
1218  // (1) First time, create oscillatorStores and read data
1219  CheckForTablesCreated();
1220 
1221  // (2) Check if the material/Z couple has been already included
1222  std::pair<const G4Material*,G4int> theKey = std::make_pair(mat,Z);
1223  if (atomTablePerMolecule->count(theKey))
1224  return atomTablePerMolecule->find(theKey)->second;
1225 
1226  // (3) If we are here, it means that we have to create the table for the material
1227  BuildOscillatorTable(mat);
1228 
1229  // (4) now, the oscillator store should be ok
1230  if (atomTablePerMolecule->count(theKey))
1231  return atomTablePerMolecule->find(theKey)->second;
1232  else
1233  {
1234  G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
1235  G4cout << "Impossible to retrieve the number of atoms per molecule for Z = "
1236  << Z << " in material " << mat->GetName() << G4endl;
1237  return 0;
1238  }
1239 }
Generated on Sat May 25 2013 14:33:34 for Geant4 by   doxygen 1.8.4