source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeOscillatorManager.cc @ 1348

Last change on this file since 1348 was 1347, checked in by garnier, 14 years ago

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