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