Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PenelopeOscillatorManager.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4PenelopeOscillatorManager.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4PenelopeOscillatorManager.cc (Version 9.5.p1)


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