Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60
61G4PenelopeOscillatorManager::G4PenelopeOscillatorManager() :
62 oscillatorStoreIonisation(0),oscillatorStoreCompton(0),atomicNumber(0),
63 atomicMass(0),excitationEnergy(0),plasmaSquared(0),atomsPerMolecule(0),
64 atomTablePerMolecule(0)
65{
66 fReadElementData = false;
67 for (G4int i=0;i<5;i++)
68 {
69 for (G4int j=0;j<2000;j++)
70 elementData[i][j] = 0.;
71 }
72 verbosityLevel = 0;
73}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77G4PenelopeOscillatorManager::~G4PenelopeOscillatorManager()
78{
79 Clear();
80 delete instance;
81}
82
83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84
85G4PenelopeOscillatorManager* G4PenelopeOscillatorManager::instance = 0;
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88
89G4PenelopeOscillatorManager* G4PenelopeOscillatorManager::GetOscillatorManager()
90{
91 if (!instance)
92 instance = new G4PenelopeOscillatorManager();
93 return instance;
94}
95
96//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97
98void G4PenelopeOscillatorManager::Clear()
99{
100 if (verbosityLevel > 1)
101 G4cout << " G4PenelopeOscillatorManager::Clear() - Clean Oscillator Tables" << G4endl;
102
103 //Clean up OscillatorStoreIonisation
104 std::map<const G4Material*,G4PenelopeOscillatorTable*>::iterator i;
105 for (i=oscillatorStoreIonisation->begin();i != oscillatorStoreIonisation->end();i++)
106 {
107 G4PenelopeOscillatorTable* table = i->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 (i=oscillatorStoreCompton->begin();i != oscillatorStoreCompton->end();i++)
122 {
123 G4PenelopeOscillatorTable* table = i->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 NULL;
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 NULL;
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 NULL;
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 NULL;
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]->GetZ() << ") --> " <<
490 (*StechiometricFactors)[i] << " atoms/molecule " << G4endl;
491 }
492 }
493
494
495 for (G4int i=0;i<nElements;i++)
496 {
497 G4int iZ = (G4int) (*elementVector)[i]->GetZ();
498 totalZ += iZ * (*StechiometricFactors)[i];
499 totalMolecularWeight += (*elementVector)[i]->GetAtomicMassAmu() * (*StechiometricFactors)[i];
500 meanExcitationEnergy += iZ*std::log(meanAtomExcitationEnergy[iZ-1])*(*StechiometricFactors)[i];
501 /*
502 G4cout << iZ << " " << (*StechiometricFactors)[i] << " " << totalZ << " " <<
503 totalMolecularWeight/(g/mole) << " " << meanExcitationEnergy << " " <<
504 meanAtomExcitationEnergy[iZ-1]/eV <<
505 G4endl;
506 */
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 = std::exp(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 if (std::fabs(occup) > 0)
559 {
560 G4PenelopeOscillator newOscLocal;
561 newOscLocal.SetOscillatorStrength(std::fabs(occup)*(*StechiometricFactors)[k]);
562 newOscLocal.SetIonisationEnergy(elementData[3][i]);
563 newOscLocal.SetHartreeFactor(elementData[4][i]/fine_structure_const);
564 newOscLocal.SetParentZ(elementData[0][i]);
565 //keep track of the origianl shell level
566 newOscLocal.SetParentShellID((G4int)elementData[1][i]);
567 //register only K, L and M shells. Outer shells all grouped with
568 //shellIndex = 30
569 if (elementData[0][i] > 6 && elementData[1][i] < 10)
570 newOscLocal.SetShellFlag(((G4int)elementData[1][i]));
571 else
572 newOscLocal.SetShellFlag(30);
573 helper->push_back(newOscLocal);
574 if (occup < 0)
575 {
576 G4double ff = (*helper)[0].GetOscillatorStrength();
577 ff += std::fabs(occup)*(*StechiometricFactors)[k];
578 (*helper)[0].SetOscillatorStrength(ff);
579 }
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 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 strenght and energy of plasmon: " <<
610 conductionStrength << " and " << plasmaEnergy/eV << " eV" << G4endl;
611 }
612
613 if (conductionStrength < 0.5 || 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 strenght.. 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
699 //Sternheimer's adjustment factor
700 G4double adjustmentFactor = 0;
701 if (helper->size() > 1)
702 {
703 G4double TST = totalZ*std::log(meanExcitationEnergy/eV);
704 G4double AALow = 0.5;
705 G4double AAHigh = 10.;
706 do
707 {
708 adjustmentFactor = (AALow+AAHigh)*0.5;
709 G4double sumLocal = 0;
710 for (size_t i=0;i<helper->size();i++)
711 {
712 if (i == 0 && isAConductor)
713 {
714 G4double resEne = (*helper)[i].GetResonanceEnergy();
715 sumLocal += (*helper)[i].GetOscillatorStrength()*std::log(resEne/eV);
716 }
717 else
718 {
719 G4double ionEne = (*helper)[i].GetIonisationEnergy();
720 G4double oscStre = (*helper)[i].GetOscillatorStrength();
721 G4double WI2 = (adjustmentFactor*adjustmentFactor*ionEne*ionEne) +
722 2./3.*(oscStre/totalZ)*Omega*Omega;
723 G4double resEne = std::sqrt(WI2);
724 (*helper)[i].SetResonanceEnergy(resEne);
725 sumLocal += (*helper)[i].GetOscillatorStrength()*std::log(resEne/eV);
726 }
727 }
728 if (sumLocal < TST)
729 AALow = adjustmentFactor;
730 else
731 AAHigh = adjustmentFactor;
732 }while((AAHigh-AALow)>(1e-14*adjustmentFactor));
733 }
734 else
735 {
736 G4double ionEne = (*helper)[0].GetIonisationEnergy();
737 (*helper)[0].SetIonisationEnergy(std::fabs(ionEne));
738 (*helper)[0].SetResonanceEnergy(meanExcitationEnergy);
739 }
740 if (verbosityLevel > 1)
741 {
742 G4cout << "Sternheimer's adjustment factor: " << adjustmentFactor << G4endl;
743 }
744
745 //Check again for data consistency
746 G4double xcheck = (*helper)[0].GetOscillatorStrength()*std::log((*helper)[0].GetResonanceEnergy());
747 G4double TST = (*helper)[0].GetOscillatorStrength();
748 for (size_t i=1;i<helper->size();i++)
749 {
750 xcheck += (*helper)[i].GetOscillatorStrength()*std::log((*helper)[i].GetResonanceEnergy());
751 TST += (*helper)[i].GetOscillatorStrength();
752 }
753 if (std::fabs(TST-totalZ)>1e-8*totalZ)
754 {
755 G4ExceptionDescription ed;
756 ed << "Inconsistent oscillator data " << G4endl;
757 ed << TST << " " << totalZ << G4endl;
758 G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
759 "em2036",FatalException,ed);
760 }
761 xcheck = std::exp(xcheck/totalZ);
762 if (std::fabs(xcheck-meanExcitationEnergy) > 1e-8*meanExcitationEnergy)
763 {
764 G4ExceptionDescription ed;
765 ed << "Error in Sterheimer factor calculation " << G4endl;
766 ed << xcheck/eV << " " << meanExcitationEnergy/eV << G4endl;
767 G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
768 "em2037",FatalException,ed);
769 }
770
771 //Selection of the lowest ionisation energy for inner shells. Only the K, L and M shells with
772 //ionisation energy less than the N1 shell of the heaviest element in the material are considered as
773 //inner shells. As a results, the inner/outer shell character of an atomic shell depends on the
774 //composition of the material.
775 G4double Zmax = 0;
776 for (G4int k=0;k<nElements;k++)
777 {
778 G4double Z = (*elementVector)[k]->GetZ();
779 if (Z>Zmax) Zmax = Z;
780 }
781 //Find N1 level of the heaviest element (if any).
782 G4bool found = false;
783 G4double cutEnergy = 50*eV;
784 for (size_t i=0;i<helper->size() && !found;i++)
785 {
786 G4double Z = (*helper)[i].GetParentZ();
787 G4int shID = (*helper)[i].GetParentShellID(); //look for the N1 level
788 if (shID == 10 && Z == Zmax)
789 {
790 found = true;
791 if ((*helper)[i].GetIonisationEnergy() > cutEnergy)
792 cutEnergy = (*helper)[i].GetIonisationEnergy();
793 }
794 }
795 //Make that cutEnergy cannot be higher than 250 eV, namely the fluorescence level by
796 //Geant4
797 G4double lowEnergyLimitForFluorescence = 250*eV;
798 cutEnergy = std::min(cutEnergy,lowEnergyLimitForFluorescence);
799
800 if (verbosityLevel > 1)
801 G4cout << "Cutoff energy: " << cutEnergy/eV << " eV" << G4endl;
802
803 //
804 //Copy helper in the oscillatorTable for Ionisation
805 //
806 //Oscillator table Ionisation for the material
807 G4PenelopeOscillatorTable* theTable = new G4PenelopeOscillatorTable(); //vector of oscillator
808 G4PenelopeOscillatorResEnergyComparator comparator;
809 std::sort(helper->begin(),helper->end(),comparator);
810
811 //COPY THE HELPER (vector of object) to theTable (vector of Pointers).
812 for (size_t i=0;i<helper->size();i++)
813 {
814 //copy content --> one may need it later (e.g. to fill an other table, with variations)
815 G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
816 theTable->push_back(theOsc);
817 }
818
819 //Oscillators of outer shells with resonance energies differing by a factor less than
820 //Rgroup are grouped as a single oscillator
821 G4double Rgroup = 1.05;
822 size_t Nost = theTable->size();
823
824 size_t firstIndex = (isAConductor) ? 1 : 0; //for conductors, skip conduction oscillator
825 G4bool loopAgain = false;
826 G4int removedLevels = 0;
827 do
828 {
829 loopAgain = false;
830 if (Nost>firstIndex+1)
831 {
832 removedLevels = 0;
833 for (size_t i=firstIndex;i<theTable->size()-1;i++)
834 {
835 G4bool skipLoop = false;
836 G4int shellFlag = (*theTable)[i]->GetShellFlag();
837 G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
838 G4double resEne = (*theTable)[i]->GetResonanceEnergy();
839 G4double resEnePlus1 = (*theTable)[i+1]->GetResonanceEnergy();
840 G4double oscStre = (*theTable)[i]->GetOscillatorStrength();
841 G4double oscStrePlus1 = (*theTable)[i+1]->GetOscillatorStrength();
842 //if (shellFlag < 10 && ionEne>cutEnergy) in Penelope
843 if (ionEne>cutEnergy) //remove condition that shellFlag < 10!
844 skipLoop = true;
845 if (resEne<1.0*eV || resEnePlus1<1.0*eV)
846 skipLoop = true;
847 if (resEnePlus1 > Rgroup*resEne)
848 skipLoop = true;
849 if (!skipLoop)
850 {
851 G4double newRes = std::exp((oscStre*std::log(resEne)+
852 oscStrePlus1*std::log(resEnePlus1))
853 /(oscStre+oscStrePlus1));
854 (*theTable)[i]->SetResonanceEnergy(newRes);
855 G4double newIon = (oscStre*ionEne+
856 oscStrePlus1*(*theTable)[i+1]->GetIonisationEnergy())/
857 (oscStre+oscStrePlus1);
858 (*theTable)[i]->SetIonisationEnergy(newIon);
859 G4double newStre = oscStre+oscStrePlus1;
860 (*theTable)[i]->SetOscillatorStrength(newStre);
861 G4double newHartree = (oscStre*(*theTable)[i]->GetHartreeFactor()+
862 oscStrePlus1*(*theTable)[i+1]->GetHartreeFactor())/
863 (oscStre+oscStrePlus1);
864 (*theTable)[i]->SetHartreeFactor(newHartree);
865 if ((*theTable)[i]->GetParentZ() != (*theTable)[i+1]->GetParentZ())
866 (*theTable)[i]->SetParentZ(0.);
867 if (shellFlag < 10 || (*theTable)[i+1]->GetShellFlag() < 10)
868 {
869 G4int newFlag = std::min(shellFlag,(*theTable)[i+1]->GetShellFlag());
870 (*theTable)[i]->SetShellFlag(newFlag);
871 }
872 else
873 (*theTable)[i]->SetShellFlag(30);
874 //We've lost anyway the track of the original level
875 (*theTable)[i]->SetParentShellID((*theTable)[i]->GetShellFlag());
876
877
878 if (i<theTable->size()-2)
879 {
880 for (size_t ii=i+1;ii<theTable->size()-1;ii++)
881 (*theTable)[ii] = (*theTable)[ii+1];
882 }
883 //G4cout << theTable->size() << G4endl;
884 theTable->erase(theTable->begin()+theTable->size()-1); //delete last element
885 removedLevels++;
886 }
887 }
888 }
889 if (removedLevels)
890 {
891 Nost -= removedLevels;
892 loopAgain = true;
893 }
894 if (Rgroup < 1.414213 || Nost > 64)
895 {
896 Rgroup = Rgroup*Rgroup;
897 loopAgain = true;
898 }
899 }while(loopAgain);
900
901 if (verbosityLevel > 1)
902 {
903 G4cout << "Final grouping factor for Ionisation: " << Rgroup << G4endl;
904 }
905
906 //Final Electron/Positron model parameters
907 for (size_t i=0;i<theTable->size();i++)
908 {
909 //Set cutoff recoil energy for the resonant mode
910 G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
911 if (ionEne < 1e-3*eV)
912 {
913 G4double resEne = (*theTable)[i]->GetResonanceEnergy();
914 (*theTable)[i]->SetIonisationEnergy(0.*eV);
915 (*theTable)[i]->SetCutoffRecoilResonantEnergy(resEne);
916 }
917 else
918 (*theTable)[i]->SetCutoffRecoilResonantEnergy(ionEne);
919 }
920
921 //Last step
922 oscillatorStoreIonisation->insert(std::make_pair(material,theTable));
923
924
925 /*
926 SAME FOR COMPTON
927 */
928 //
929 //Copy helper in the oscillatorTable for Compton
930 //
931 //Oscillator table Ionisation for the material
932 G4PenelopeOscillatorTable* theTableC = new G4PenelopeOscillatorTable(); //vector of oscillator
933 //order by ionisation energy
934 std::sort(helper->begin(),helper->end());
935 //COPY THE HELPER (vector of object) to theTable (vector of Pointers).
936 for (size_t i=0;i<helper->size();i++)
937 {
938 //copy content --> one may need it later (e.g. to fill an other table, with variations)
939 G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
940 theTableC->push_back(theOsc);
941 }
942 //Oscillators of outer shells with resonance energies differing by a factor less than
943 //Rgroup are grouped as a single oscillator
944 Rgroup = 1.5;
945 Nost = theTableC->size();
946
947 firstIndex = (isAConductor) ? 1 : 0; //for conductors, skip conduction oscillator
948 loopAgain = false;
949 removedLevels = 0;
950 do
951 {
952 loopAgain = false;
953 if (Nost>firstIndex+1)
954 {
955 removedLevels = 0;
956 for (size_t i=firstIndex;i<theTableC->size()-1;i++)
957 {
958 G4bool skipLoop = false;
959 //G4int shellFlag = (*theTableC)[i]->GetShellFlag();
960 G4double ionEne = (*theTableC)[i]->GetIonisationEnergy();
961 G4double ionEnePlus1 = (*theTableC)[i+1]->GetIonisationEnergy();
962 G4double oscStre = (*theTableC)[i]->GetOscillatorStrength();
963 G4double oscStrePlus1 = (*theTableC)[i+1]->GetOscillatorStrength();
964 //if (shellFlag < 10 && ionEne>cutEnergy) in Penelope
965 if (ionEne>cutEnergy)
966 skipLoop = true;
967 if (ionEne<1.0*eV || ionEnePlus1<1.0*eV)
968 skipLoop = true;
969 if (ionEnePlus1 > Rgroup*ionEne)
970 skipLoop = true;
971
972 if (!skipLoop)
973 {
974 G4double newIon = (oscStre*ionEne+
975 oscStrePlus1*ionEnePlus1)/
976 (oscStre+oscStrePlus1);
977 (*theTableC)[i]->SetIonisationEnergy(newIon);
978 G4double newStre = oscStre+oscStrePlus1;
979 (*theTableC)[i]->SetOscillatorStrength(newStre);
980 G4double newHartree = (oscStre*(*theTableC)[i]->GetHartreeFactor()+
981 oscStrePlus1*(*theTableC)[i+1]->GetHartreeFactor())/
982 (oscStre+oscStrePlus1);
983 (*theTableC)[i]->SetHartreeFactor(newHartree);
984 if ((*theTableC)[i]->GetParentZ() != (*theTableC)[i+1]->GetParentZ())
985 (*theTableC)[i]->SetParentZ(0.);
986 (*theTableC)[i]->SetShellFlag(30);
987 (*theTableC)[i]->SetParentShellID((*theTableC)[i]->GetShellFlag());
988
989 if (i<theTableC->size()-2)
990 {
991 for (size_t ii=i+1;ii<theTableC->size()-1;ii++)
992 (*theTableC)[ii] = (*theTableC)[ii+1];
993 }
994 theTableC->erase(theTableC->begin()+theTableC->size()-1); //delete last element
995 removedLevels++;
996 }
997 }
998 }
999 if (removedLevels)
1000 {
1001 Nost -= removedLevels;
1002 loopAgain = true;
1003 }
1004 if (Rgroup < 2.0 || Nost > 64)
1005 {
1006 Rgroup = Rgroup*Rgroup;
1007 loopAgain = true;
1008 }
1009 }while(loopAgain);
1010
1011
1012 if (verbosityLevel > 1)
1013 {
1014 G4cout << "Final grouping factor for Compton: " << Rgroup << G4endl;
1015 }
1016
1017 //Last step
1018 oscillatorStoreCompton->insert(std::make_pair(material,theTableC));
1019
1020 /* //TESTING PURPOSES
1021 if (verbosityLevel > 1)
1022 {
1023 G4cout << "The table contains " << helper->size() << " oscillators " << G4endl;
1024 for (size_t k=0;k<helper->size();k++)
1025 {
1026 G4cout << "Oscillator # " << k << G4endl;
1027 G4cout << "Z = " << (*helper)[k].GetParentZ() << G4endl;
1028 G4cout << "Shell Flag = " << (*helper)[k].GetShellFlag() << G4endl;
1029 G4cout << "Compton index = " << (*helper)[k].GetHartreeFactor() << G4endl;
1030 G4cout << "Ionisation energy = " << (*helper)[k].GetIonisationEnergy()/eV << " eV" << G4endl;
1031 G4cout << "Occupation number = " << (*helper)[k].GetOscillatorStrength() << G4endl;
1032 G4cout << "Resonance energy = " << (*helper)[k].GetResonanceEnergy()/eV << " eV" << G4endl;
1033 }
1034
1035 for (size_t k=0;k<helper->size();k++)
1036 {
1037 G4cout << k << " " << (*helper)[k].GetOscillatorStrength() << " " <<
1038 (*helper)[k].GetIonisationEnergy()/eV << " " << (*helper)[k].GetResonanceEnergy()/eV << " " <<
1039 (*helper)[k].GetParentZ() << " " << (*helper)[k].GetShellFlag() << " " <<
1040 (*helper)[k].GetHartreeFactor() << G4endl;
1041 }
1042 }
1043 */
1044
1045
1046 //CLEAN UP theHelper and its content
1047 delete helper;
1048 if (verbosityLevel > 1)
1049 Dump(material);
1050
1051 return;
1052}
1053
1054//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1055
1056void G4PenelopeOscillatorManager::ReadElementData()
1057{
1058 if (verbosityLevel > 0)
1059 {
1060 G4cout << "G4PenelopeOscillatorManager::ReadElementData()" << G4endl;
1061 G4cout << "Going to read Element Data" << G4endl;
1062 }
1063 char* path = getenv("G4LEDATA");
1064 if (!path)
1065 {
1066 G4String excep = "G4PenelopeOscillatorManager - G4LEDATA environment variable not set!";
1067 G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
1068 "em0006",FatalException,excep);
1069 return;
1070 }
1071 G4String pathString(path);
1072 G4String pathFile = pathString + "/penelope/pdatconf.p08";
1073 std::ifstream file(pathFile);
1074
1075 if (!file.is_open())
1076 {
1077 G4String excep = "G4PenelopeOscillatorManager - data file " + pathFile + " not found!";
1078 G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
1079 "em0003",FatalException,excep);
1080 }
1081
1082 G4AtomicTransitionManager* theTransitionManager =
1083 G4AtomicTransitionManager::Instance();
1084
1085 //Read header (22 lines)
1086 G4String theHeader;
1087 for (G4int iline=0;iline<22;iline++)
1088 getline(file,theHeader);
1089 //Done
1090 G4int Z=0;
1091 G4int shellCode = 0;
1092 G4String shellId = "NULL";
1093 G4int occupationNumber = 0;
1094 G4double ionisationEnergy = 0.0*eV;
1095 G4double hartreeProfile = 0.;
1096 G4int shellCounter = 0;
1097 G4int oldZ = -1;
1098 G4int numberOfShells = 0;
1099 //Start reading data
1100 for (G4int i=0;!file.eof();i++)
1101 {
1102 file >> Z >> shellCode >> shellId >> occupationNumber >> ionisationEnergy >> hartreeProfile;
1103 if (Z>0 && i<2000)
1104 {
1105 elementData[0][i] = Z;
1106 elementData[1][i] = shellCode;
1107 elementData[2][i] = occupationNumber;
1108 //reset things
1109 if (Z != oldZ)
1110 {
1111 shellCounter = 0;
1112 oldZ = Z;
1113 numberOfShells = theTransitionManager->NumberOfShells(Z);
1114 }
1115 G4double bindingEnergy = -1*eV;
1116 if (shellCounter<numberOfShells)
1117 {
1118 G4AtomicShell* shell = theTransitionManager->Shell(Z,shellCounter);
1119 bindingEnergy = shell->BindingEnergy();
1120 }
1121 //Valid level found in the G4AtomicTransition database: keep it, otherwise use
1122 //the ionisation energy found in the Penelope database
1123 elementData[3][i] = (bindingEnergy>100*eV) ? bindingEnergy : ionisationEnergy*eV;
1124 //elementData[3][i] = ionisationEnergy*eV;
1125 elementData[4][i] = hartreeProfile;
1126 shellCounter++;
1127 }
1128 }
1129 file.close();
1130
1131 if (verbosityLevel > 1)
1132 {
1133 G4cout << "G4PenelopeOscillatorManager::ReadElementData(): Data file read" << G4endl;
1134 }
1135 fReadElementData = true;
1136 return;
1137
1138}
1139
1140//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1141G4double G4PenelopeOscillatorManager::GetMeanExcitationEnergy(const G4Material* mat)
1142{
1143 // (1) First time, create oscillatorStores and read data
1144 CheckForTablesCreated();
1145
1146 // (2) Check if the material has been already included
1147 if (excitationEnergy->count(mat))
1148 return excitationEnergy->find(mat)->second;
1149
1150 // (3) If we are here, it means that we have to create the table for the material
1151 BuildOscillatorTable(mat);
1152
1153 // (4) now, the oscillator store should be ok
1154 if (excitationEnergy->count(mat))
1155 return excitationEnergy->find(mat)->second;
1156 else
1157 {
1158 G4cout << "G4PenelopeOscillatorManager::GetMolecularExcitationEnergy() " << G4endl;
1159 G4cout << "Impossible to retrieve the excitation energy for " << mat->GetName() << G4endl;
1160 return 0;
1161 }
1162}
1163
1164//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1165G4double G4PenelopeOscillatorManager::GetPlasmaEnergySquared(const G4Material* mat)
1166{
1167 // (1) First time, create oscillatorStores and read data
1168 CheckForTablesCreated();
1169
1170 // (2) Check if the material has been already included
1171 if (plasmaSquared->count(mat))
1172 return plasmaSquared->find(mat)->second;
1173
1174 // (3) If we are here, it means that we have to create the table for the material
1175 BuildOscillatorTable(mat);
1176
1177 // (4) now, the oscillator store should be ok
1178 if (plasmaSquared->count(mat))
1179 return plasmaSquared->find(mat)->second;
1180 else
1181 {
1182 G4cout << "G4PenelopeOscillatorManager::GetPlasmaEnergySquared() " << G4endl;
1183 G4cout << "Impossible to retrieve the plasma energy for " << mat->GetName() << G4endl;
1184 return 0;
1185 }
1186}
1187
1188//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1189
1190G4double G4PenelopeOscillatorManager::GetAtomsPerMolecule(const G4Material* mat)
1191{
1192 // (1) First time, create oscillatorStores and read data
1193 CheckForTablesCreated();
1194
1195 // (2) Check if the material has been already included
1196 if (atomsPerMolecule->count(mat))
1197 return atomsPerMolecule->find(mat)->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 (atomsPerMolecule->count(mat))
1204 return atomsPerMolecule->find(mat)->second;
1205 else
1206 {
1207 G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
1208 G4cout << "Impossible to retrieve the number of atoms per molecule for "
1209 << mat->GetName() << G4endl;
1210 return 0;
1211 }
1212}
1213
1214//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1215
1216G4double G4PenelopeOscillatorManager::GetNumberOfZAtomsPerMolecule(const G4Material* mat,G4int Z)
1217{
1218 // (1) First time, create oscillatorStores and read data
1219 CheckForTablesCreated();
1220
1221 // (2) Check if the material/Z couple has been already included
1222 std::pair<const G4Material*,G4int> theKey = std::make_pair(mat,Z);
1223 if (atomTablePerMolecule->count(theKey))
1224 return atomTablePerMolecule->find(theKey)->second;
1225
1226 // (3) If we are here, it means that we have to create the table for the material
1227 BuildOscillatorTable(mat);
1228
1229 // (4) now, the oscillator store should be ok
1230 if (atomTablePerMolecule->count(theKey))
1231 return atomTablePerMolecule->find(theKey)->second;
1232 else
1233 {
1234 G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
1235 G4cout << "Impossible to retrieve the number of atoms per molecule for Z = "
1236 << Z << " in material " << mat->GetName() << G4endl;
1237 return 0;
1238 }
1239}
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4PenelopeOscillatorTable
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
Definition: G4PenelopeOscillatorManager.hh:65
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT 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:174
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:163
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:189
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:208
G4Material::GetFractionVector
const G4double * GetFractionVector() const
Definition: G4Material.hh:193
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:177
G4PenelopeOscillatorManager::GetNumberOfZAtomsPerMolecule
G4double GetNumberOfZAtomsPerMolecule(const G4Material *, G4int Z)
Definition: G4PenelopeOscillatorManager.cc:1216
G4PenelopeOscillatorManager::GetAtomsPerMolecule
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
Definition: G4PenelopeOscillatorManager.cc:1190
G4PenelopeOscillatorManager::GetOscillatorTableCompton
G4PenelopeOscillatorTable * GetOscillatorTableCompton(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:363
G4PenelopeOscillatorManager::GetOscillatorManager
static G4PenelopeOscillatorManager * GetOscillatorManager()
Definition: G4PenelopeOscillatorManager.cc:89
G4PenelopeOscillatorManager::GetOscillatorTableIonisation
G4PenelopeOscillatorTable * GetOscillatorTableIonisation(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:316
G4PenelopeOscillatorManager::GetPlasmaEnergySquared
G4double GetPlasmaEnergySquared(const G4Material *)
Returns the squared plasma energy.
Definition: G4PenelopeOscillatorManager.cc:1165
G4PenelopeOscillatorManager::GetTotalZ
G4double GetTotalZ(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:266
G4PenelopeOscillatorManager::GetMeanExcitationEnergy
G4double GetMeanExcitationEnergy(const G4Material *)
Returns the mean excitation energy.
Definition: G4PenelopeOscillatorManager.cc:1141
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
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:37