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