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

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

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