source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeOscillatorManager.cc

Last change on this file was 1347, checked in by garnier, 15 years ago

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