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