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