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