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