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

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

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