| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Authors: Luciano Pandola (luciano.pandola a
27 //
28 // History:
29 // -----------
30 //
31 // 03 Dec 2009 First working version, Lucian
32 // 16 Feb 2010 Added methods to store also t
33 // molecule, Luciano Pandola
34 // 19 Feb 2010 Scale the Hartree factors in
35 // table by (1/fine_structure_co
36 // always the ratio (hartreeFact
37 // 16 Mar 2010 Added methods to calculate an
38 // and plasma energy (used for I
39 // 18 Mar 2010 Added method to retrieve numb
40 // molecule. L. Pandola
41 // 06 Sep 2011 Override the local Penelope d
42 // G4AtomicDeexcitation database
43 // binding energies. L. Pandola
44 // 15 Mar 2012 Added method to retrieve numb
45 // molecule. Restore the origina
46 // below 100 eV. L. Pandola
47 //
48 // -------------------------------------------
49
50 #include "G4PenelopeOscillatorManager.hh"
51
52 #include "globals.hh"
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4AtomicTransitionManager.hh"
56 #include "G4AtomicShell.hh"
57 #include "G4Material.hh"
58 #include "G4Exp.hh"
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61
62 G4ThreadLocal G4PenelopeOscillatorManager* G4P
63
64 //....oooOO0OOooo........oooOO0OOooo........oo
65
66 G4PenelopeOscillatorManager::G4PenelopeOscilla
67 fOscillatorStoreIonisation(nullptr),fOscilla
68 fAtomicNumber(nullptr),fAtomicMass(nullptr),
69 fPlasmaSquared(nullptr),fAtomsPerMolecule(nu
70 fAtomTablePerMolecule(nullptr)
71 {
72 fReadElementData = false;
73 for (G4int i=0;i<5;i++)
74 {
75 for (G4int j=0;j<2000;j++)
76 fElementData[i][j] = 0.;
77 }
78 fVerbosityLevel = 0;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 G4PenelopeOscillatorManager::~G4PenelopeOscill
84 {
85 Clear();
86 delete instance;
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90
91 G4PenelopeOscillatorManager* G4PenelopeOscilla
92 {
93 if (!instance)
94 instance = new G4PenelopeOscillatorManager
95 return instance;
96 }
97
98 //....oooOO0OOooo........oooOO0OOooo........oo
99
100 void G4PenelopeOscillatorManager::Clear()
101 {
102 if (fVerbosityLevel > 1)
103 G4cout << " G4PenelopeOscillatorManager::C
104
105 //Clean up OscillatorStoreIonisation
106 for (auto& item : (*fOscillatorStoreIonisati
107 {
108 G4PenelopeOscillatorTable* table = item.
109 if (table)
110 {
111 for (std::size_t k=0;k<table->size();++k)
112 {
113 if ((*table)[k])
114 delete ((*table)[k]);
115 }
116 delete table;
117 }
118 }
119 delete fOscillatorStoreIonisation;
120
121 //Clean up OscillatorStoreCompton
122 for (auto& item : (*fOscillatorStoreCompton)
123 {
124 G4PenelopeOscillatorTable* table = item.
125 if (table)
126 {
127 for (std::size_t k=0;k<table->size();++k)
128 {
129 if ((*table)[k])
130 delete ((*table)[k]);
131 }
132 delete table;
133 }
134 }
135 delete fOscillatorStoreCompton;
136
137 if (fAtomicMass) delete fAtomicMass;
138 if (fAtomicNumber) delete fAtomicNumber;
139 if (fExcitationEnergy) delete fExcitationEne
140 if (fPlasmaSquared) delete fPlasmaSquared;
141 if (fAtomsPerMolecule) delete fAtomsPerMolec
142 if (fAtomTablePerMolecule) delete fAtomTable
143 }
144
145 //....oooOO0OOooo........oooOO0OOooo........oo
146
147 void G4PenelopeOscillatorManager::Dump(const G
148 {
149 G4PenelopeOscillatorTable* theTable = GetOsc
150 if (!theTable)
151 {
152 G4cout << " G4PenelopeOscillatorManager:
153 G4cout << "Problem in retrieving the Ion
154 << material->GetName() << G4endl;
155 return;
156 }
157 G4cout << "*********************************
158 G4cout << " Penelope Oscillator Table Ionisa
159 G4cout << "*********************************
160 G4cout << "The table contains " << theTable-
161 G4cout << "*********************************
162 if (theTable->size() < 10)
163 for (std::size_t k=0;k<theTable->size();++
164 {
165 G4cout << "Oscillator # " << k << " Z = " <<
166 " Shell Flag = " << (*theTable)[k]->GetShe
167 " Parent shell ID = " << (*theTable)[k]->G
168 G4cout << "Ionisation energy = " << (*theTab
169 G4cout << "Occupation number = " << (*theTab
170 G4cout << "Resonance energy = " << (*theTabl
171 G4cout << "Cufoff resonance energy = " <<
172 (*theTable)[k]->GetCutoffRecoilResonantEne
173 G4cout << "*********************************
174 }
175 for (std::size_t k=0;k<theTable->size();++k)
176 {
177 G4cout << k << " " << (*theTable)[k]->Ge
178 (*theTable)[k]->GetIonisationEnergy()/eV <<
179 << (*theTable)[k]->GetResonanceEnergy()
180 (*theTable)[k]->GetParentZ() << " " << (*the
181 (*theTable)[k]->GetParentShellID() << G4endl
182 }
183 G4cout << "*********************************
184
185 //Compton table
186 theTable = GetOscillatorTableCompton(materia
187 if (!theTable)
188 {
189 G4cout << " G4PenelopeOscillatorManager:
190 G4cout << "Problem in retrieving the Com
191 material->GetName() << G4endl;
192 return;
193 }
194 G4cout << "*********************************
195 G4cout << " Penelope Oscillator Table Compto
196 G4cout << "*********************************
197 G4cout << "The table contains " << theTable-
198 G4cout << "*********************************
199 if (theTable->size() < 10)
200 for (std::size_t k=0;k<theTable->size();++
201 {
202 G4cout << "Oscillator # " << k << " Z = " <<
203 " Shell Flag = " << (*theTable)[k]->GetShe
204 " Parent shell ID = " << (*theTable)[k]->
205 G4cout << "Compton index = " << (*theTable)[
206 G4cout << "Ionisation energy = " << (*theTab
207 G4cout << "Occupation number = " << (*theTab
208 G4cout << "*********************************
209 }
210 for (std::size_t k=0;k<theTable->size();++k)
211 {
212 G4cout << k << " " << (*theTable)[k]->Ge
213 (*theTable)[k]->GetIonisationEnergy()/eV <<
214 (*theTable)[k]->GetParentZ() << " " << (*the
215 (*theTable)[k]->GetParentShellID() << G4endl
216 }
217 G4cout << "*********************************
218
219 return;
220 }
221
222 //....oooOO0OOooo........oooOO0OOooo........oo
223
224 void G4PenelopeOscillatorManager::CheckForTabl
225 {
226 //Tables should be created at the same time,
227 //simultaneously
228 if (!fOscillatorStoreIonisation)
229 {
230 fOscillatorStoreIonisation = new std::ma
231 if (!fReadElementData)
232 ReadElementData();
233 if (!fOscillatorStoreIonisation)
234 //It should be ok now
235 G4Exception("G4PenelopeOscillatorManager::Ge
236 "em2034",FatalException,
237 "Problem in allocating the Oscillator
238 }
239
240 if (!fOscillatorStoreCompton)
241 {
242 fOscillatorStoreCompton = new std::map<c
243 if (!fReadElementData)
244 ReadElementData();
245 if (!fOscillatorStoreCompton)
246 //It should be ok now
247 G4Exception("G4PenelopeOscillatorManager::Ge
248 "em2034",FatalException,
249 "Problem in allocating the Oscillator
250 }
251
252 if (!fAtomicNumber)
253 fAtomicNumber = new std::map<const G4Mater
254 if (!fAtomicMass)
255 fAtomicMass = new std::map<const G4Materia
256 if (!fExcitationEnergy)
257 fExcitationEnergy = new std::map<const G4M
258 if (!fPlasmaSquared)
259 fPlasmaSquared = new std::map<const G4Mate
260 if (!fAtomsPerMolecule)
261 fAtomsPerMolecule = new std::map<const G4M
262 if (!fAtomTablePerMolecule)
263 fAtomTablePerMolecule = new std::map< std:
264 }
265
266
267 //....oooOO0OOooo........oooOO0OOooo........oo
268
269 G4double G4PenelopeOscillatorManager::GetTotal
270 {
271 // (1) First time, create fOscillatorStores
272 CheckForTablesCreated();
273
274 // (2) Check if the material has been alread
275 if (fAtomicNumber->count(mat))
276 return fAtomicNumber->find(mat)->second;
277
278 // (3) If we are here, it means that we have
279 BuildOscillatorTable(mat);
280
281 // (4) now, the oscillator store should be o
282 if (fAtomicNumber->count(mat))
283 return fAtomicNumber->find(mat)->second;
284 else
285 {
286 G4cout << "G4PenelopeOscillatorManager::
287 G4cout << "Impossible to retrieve the to
288 return 0;
289 }
290 }
291
292 //....oooOO0OOooo........oooOO0OOooo........oo
293
294 G4double G4PenelopeOscillatorManager::GetTotal
295 {
296 // (1) First time, create fOscillatorStores
297 CheckForTablesCreated();
298
299 // (2) Check if the material has been alread
300 if (fAtomicMass->count(mat))
301 return fAtomicMass->find(mat)->second;
302
303 // (3) If we are here, it means that we have
304 BuildOscillatorTable(mat);
305
306 // (4) now, the oscillator store should be o
307 if (fAtomicMass->count(mat))
308 return fAtomicMass->find(mat)->second;
309 else
310 {
311 G4cout << "G4PenelopeOscillatorManager::
312 G4cout << "Impossible to retrieve the to
313 return 0;
314 }
315 }
316
317 //....oooOO0OOooo........oooOO0OOooo........oo
318
319 G4PenelopeOscillatorTable* G4PenelopeOscillato
320 {
321 // (1) First time, create fOscillatorStores
322 CheckForTablesCreated();
323
324 // (2) Check if the material has been alread
325 if (fOscillatorStoreIonisation->count(mat))
326 {
327 //Ok, it exists
328 return fOscillatorStoreIonisation->find(
329 }
330
331 // (3) If we are here, it means that we have
332 BuildOscillatorTable(mat);
333
334 // (4) now, the oscillator store should be o
335 if (fOscillatorStoreIonisation->count(mat))
336 return fOscillatorStoreIonisation->find(ma
337 else
338 {
339 G4cout << "G4PenelopeOscillatorManager::
340 G4cout << "Impossible to create ionisati
341 return nullptr;
342 }
343 }
344
345 //....oooOO0OOooo........oooOO0OOooo........oo
346
347 G4PenelopeOscillator* G4PenelopeOscillatorMana
348 G4int index)
349 {
350 G4PenelopeOscillatorTable* theTable = GetOsc
351 if (((std::size_t)index) < theTable->size())
352 return (*theTable)[index];
353 else
354 {
355 G4cout << "WARNING: Ionisation table for
356 theTable->size() << " oscillators" << G4endl
357 G4cout << "Oscillator #" << index << " c
358 G4cout << "Returning null pointer" << G4
359 return nullptr;
360 }
361 }
362
363
364 //....oooOO0OOooo........oooOO0OOooo........oo
365
366 G4PenelopeOscillatorTable* G4PenelopeOscillato
367 {
368 // (1) First time, create fOscillatorStore a
369 CheckForTablesCreated();
370
371 // (2) Check if the material has been alread
372 if (fOscillatorStoreCompton->count(mat))
373 {
374 //Ok, it exists
375 return fOscillatorStoreCompton->find(mat
376 }
377
378 // (3) If we are here, it means that we have
379 BuildOscillatorTable(mat);
380
381 // (4) now, the oscillator store should be o
382 if (fOscillatorStoreCompton->count(mat))
383 return fOscillatorStoreCompton->find(mat)-
384 else
385 {
386 G4cout << "G4PenelopeOscillatorManager::
387 G4cout << "Impossible to create Compton
388 return nullptr;
389 }
390 }
391
392 //....oooOO0OOooo........oooOO0OOooo........oo
393
394 G4PenelopeOscillator* G4PenelopeOscillatorMana
395 G4int index)
396 {
397 G4PenelopeOscillatorTable* theTable = GetOsc
398 if (((std::size_t)index) < theTable->size())
399 return (*theTable)[index];
400 else
401 {
402 G4cout << "WARNING: Compton table for ma
403 theTable->size() << " oscillators" << G4endl
404 G4cout << "Oscillator #" << index << " c
405 G4cout << "Returning null pointer" << G4
406 return nullptr;
407 }
408 }
409
410 //....oooOO0OOooo........oooOO0OOooo........oo
411
412 void G4PenelopeOscillatorManager::BuildOscilla
413 {
414 //THIS CORRESPONDS TO THE ROUTINE PEMATW of
415
416 G4double meanAtomExcitationEnergy[99] = {19.
417 82.0*eV, 95.0*eV,115.0*eV,137.0*e
418 166.0*eV,
419 173.0*eV,173.0*eV,180.0*eV,174.0*
420 216.0*eV,233.0*eV,245.0*eV,257.0*
421 311.0*eV,322.0*eV,330.0*eV,334.0*
422 343.0*eV,352.0*eV,363.0*eV,366.0*
423 424.0*eV,428.0*eV,441.0*eV,449.0*
424 488.0*eV,488.0*eV,487.0*eV,485.0*
425 491.0*eV,501.0*eV,523.0*eV,535.0*
426 580.0*eV,591.0*eV,614.0*eV,628.0*
427 684.0*eV,694.0*eV,705.0*eV,718.0*
428 757.0*eV,790.0*eV,790.0*eV,800.0*
429 830.0*eV,825.0*eV,794.0*eV,827.0*
430 878.0*eV,890.0*eV,902.0*eV,921.0*
431 966.0*eV,980.0*eV};
432
433 if (fVerbosityLevel > 0)
434 G4cout << "Going to build Oscillator Table
435
436 G4int nElements = (G4int)material->GetNumber
437 const G4ElementVector* elementVector = mater
438
439 //At the moment, there's no way in Geant4 to
440 //is defined with atom numbers or fraction o
441 const G4double* fractionVector = material->G
442
443 //Take always the composition by fraction of
444 //atoms: it is calculated by Geant4 but with
445 G4double totalZ = 0;
446 G4double totalMolecularWeight = 0;
447 G4double meanExcitationEnergy = 0;
448
449 std::vector<G4double> *StechiometricFactors
450
451 for (G4int i=0;i<nElements;i++)
452 {
453 //G4int iZ = (G4int) (*elementVector)[i]
454 G4double fraction = fractionVector[i];
455 G4double atomicWeigth = (*elementVector)
456 StechiometricFactors->push_back(fraction
457 }
458 //Find max
459 G4double MaxStechiometricFactor = 0.;
460 for (G4int i=0;i<nElements;i++)
461 {
462 if ((*StechiometricFactors)[i] > MaxStec
463 MaxStechiometricFactor = (*StechiometricFact
464 }
465 if (MaxStechiometricFactor<1e-16)
466 {
467 G4ExceptionDescription ed;
468 ed << "Problem with the mass composition
469 ed << "MaxStechiometricFactor = " << Max
470 G4Exception("G4PenelopeOscillatorManager
471 "em2035",FatalException,ed);
472 }
473 //Normalize
474 for (G4int i=0;i<nElements;++i)
475 (*StechiometricFactors)[i] /= MaxStechiom
476
477 // Equivalent atoms per molecule
478 G4double theatomsPerMolecule = 0;
479 for (G4int i=0;i<nElements;i++)
480 theatomsPerMolecule += (*StechiometricFact
481 G4double moleculeDensity =
482 material->GetTotNbOfAtomsPerVolume()/theat
483
484 if (fVerbosityLevel > 1)
485 {
486 for (std::size_t i=0;i<StechiometricFact
487 {
488 G4cout << "Element " << (*elementVector)[i
489 (*elementVector)[i]->GetZasInt() << ") -
490 (*StechiometricFactors)[i] << " atoms/mo
491 }
492 }
493
494 for (G4int i=0;i<nElements;++i)
495 {
496 G4int iZ = (*elementVector)[i]->GetZasIn
497 totalZ += iZ * (*StechiometricFactors)[i
498 totalMolecularWeight += (*elementVector)
499 meanExcitationEnergy += iZ*G4Log(meanAto
500 std::pair<const G4Material*,G4int> theKe
501 if (!fAtomTablePerMolecule->count(theKey
502 fAtomTablePerMolecule->insert(std::make_pair
503 }
504 meanExcitationEnergy = G4Exp(meanExcitationE
505
506 fAtomicNumber->insert(std::make_pair(materia
507 fAtomicMass->insert(std::make_pair(material,
508 fExcitationEnergy->insert(std::make_pair(mat
509 fAtomsPerMolecule->insert(std::make_pair(mat
510
511 if (fVerbosityLevel > 1)
512 {
513 G4cout << "Calculated mean excitation en
514 " = " << meanExcitationEnergy/eV << " eV" <<
515 }
516
517 std::vector<G4PenelopeOscillator> *helper =
518
519 //First Oscillator: conduction band. Tentati
520 //atom contributes a number of electrons equ
521 G4PenelopeOscillator newOsc;
522 newOsc.SetOscillatorStrength(0.);
523 newOsc.SetIonisationEnergy(0*eV);
524 newOsc.SetHartreeFactor(0);
525 newOsc.SetParentZ(0);
526 newOsc.SetShellFlag(30);
527 newOsc.SetParentShellID(30); //does not corr
528 helper->push_back(newOsc);
529
530 //Load elements and oscillators
531 for (G4int k=0;k<nElements;k++)
532 {
533 G4double Z = (*elementVector)[k]->GetZ()
534 G4bool finished = false;
535 for (G4int i=0;i<2000 && !finished;i++)
536 {
537 /*
538 fElementData[0][i] = Z;
539 fElementData[1][i] = shellCode;
540 fElementData[2][i] = occupationNumber;
541 fElementData[3][i] = ionisationEnergy;
542 fElementData[4][i] = hartreeProfile;
543 */
544 if (fElementData[0][i] == Z)
545 {
546 G4int shellID = (G4int) fElementData[1
547 G4double occup = fElementData[2][i];
548 if (shellID > 0)
549 {
550
551 if (std::fabs(occup) > 0)
552 {
553 G4PenelopeOscillator newOscLocal;
554 newOscLocal.SetOscillatorStrength(st
555 newOscLocal.SetIonisationEnergy(fEle
556 newOscLocal.SetHartreeFactor(fElemen
557 newOscLocal.SetParentZ(fElementData[
558 //keep track of the origianl shell l
559 newOscLocal.SetParentShellID((G4int)
560 //register only K, L and M shells. O
561 //shellIndex = 30
562 if (fElementData[0][i] > 6 && fEleme
563 newOscLocal.SetShellFlag(((G4int)fElemen
564 else
565 newOscLocal.SetShellFlag(30);
566 helper->push_back(newOscLocal);
567 if (occup < 0)
568 {
569 G4double ff = (*helper)[0].GetOscillat
570 ff += std::fabs(occup)*(*Stechiometric
571 (*helper)[0].SetOscillatorStrength(ff)
572 }
573 }
574 }
575 }
576 if (fElementData[0][i] > Z)
577 finished = true;
578 }
579 }
580
581 delete StechiometricFactors;
582
583 //NOW: sort oscillators according to increas
584 //Notice: it works because helper is a vecto
585 //vector to _pointers_
586 std::sort(helper->begin(),helper->end());
587
588 // Plasma energy and conduction band excitat
589 static const G4double RydbergEnergy = 13.605
590 G4double Omega = std::sqrt(4*pi*moleculeDens
591 G4double conductionStrength = (*helper)[0].G
592 G4double plasmaEnergy = Omega*std::sqrt(cond
593
594 fPlasmaSquared->insert(std::make_pair(materi
595
596 G4bool isAConductor = false;
597 G4int nullOsc = 0;
598
599 if (fVerbosityLevel > 1)
600 {
601 G4cout << "Estimated oscillator strength
602 conductionStrength << " and " << plasmaEnerg
603 }
604
605 if (conductionStrength < 0.01 || plasmaEnerg
606 {
607 if (fVerbosityLevel >1 )
608 G4cout << material->GetName() << " is an ins
609 //remove conduction band oscillator
610 helper->erase(helper->begin());
611 }
612 else //this is a conductor, Outer shells mov
613 {
614 if (fVerbosityLevel >1 )
615 G4cout << material->GetName() << " is a cond
616 isAConductor = true;
617 //copy the conduction strength.. The num
618 G4double conductionStrengthCopy = conduc
619 G4bool quit = false;
620 for (std::size_t i = 1; i<helper->size()
621 {
622 G4double oscStre = (*helper)[i].GetOscilla
623 //loop is repeated over here
624 if (oscStre < conductionStrengthCopy)
625 {
626 conductionStrengthCopy = conductionStr
627 (*helper)[i].SetOscillatorStrength(0.)
628 nullOsc++;
629 }
630 else //this is passed only once - no goto
631 {
632 quit = true;
633 (*helper)[i].SetOscillatorStrength(osc
634 if (std::fabs((*helper)[i].GetOscillat
635 {
636 conductionStrength += (*helper)[i].GetOs
637 (*helper)[i].SetOscillatorStrength(0.);
638 nullOsc++;
639 }
640 }
641 }
642 //Update conduction band
643 (*helper)[0].SetOscillatorStrength(condu
644 (*helper)[0].SetIonisationEnergy(0.);
645 (*helper)[0].SetResonanceEnergy(plasmaEn
646 G4double hartree = 0.75/std::sqrt(3.0*pi
647 Bohr_radius*Bohr_radius*Bohr_radius*
648 (*helper)[0].SetHartreeFactor(hartree/fi
649 }
650
651 //Check f-sum rule
652 G4double sum = 0;
653 for (std::size_t i=0;i<helper->size();++i)
654 {
655 sum += (*helper)[i].GetOscillatorStrengt
656 }
657 if (std::fabs(sum-totalZ) > (1e-6*totalZ))
658 {
659 G4ExceptionDescription ed;
660 ed << "Inconsistent oscillator data for
661 ed << sum << " " << totalZ << G4endl;
662 G4Exception("G4PenelopeOscillatorManager
663 "em2036",FatalException,ed);
664 }
665 if (std::fabs(sum-totalZ) > (1e-12*totalZ))
666 {
667 G4double fact = totalZ/sum;
668 for (std::size_t i=0;i<helper->size();++
669 {
670 G4double ff = (*helper)[i].GetOscillatorSt
671 (*helper)[i].SetOscillatorStrength(ff);
672 }
673 }
674
675 //Remove null items
676 for (G4int k=0;k<nullOsc;k++)
677 {
678 G4bool exit=false;
679 for (std::size_t i=0;i<helper->size() &&
680 {
681 if (std::fabs((*helper)[i].GetOscillatorSt
682 {
683 helper->erase(helper->begin()+i);
684 exit = true;
685 }
686 }
687 }
688
689 //Sternheimer's adjustment factor
690 G4double adjustmentFactor = 0;
691 if (helper->size() > 1)
692 {
693 G4double TST = totalZ*G4Log(meanExcitati
694 G4double AALow = 0.1;
695 G4double AAHigh = 10.;
696 do
697 {
698 adjustmentFactor = (AALow+AAHigh)*0.5;
699 G4double sumLocal = 0;
700 for (std::size_t i=0;i<helper->size();++i)
701 {
702 if (i == 0 && isAConductor)
703 {
704 G4double resEne = (*helper)[i].GetResona
705 sumLocal += (*helper)[i].GetOscillatorSt
706 }
707 else
708 {
709 G4double ionEne = (*helper)[i].GetIonisa
710 G4double oscStre = (*helper)[i].GetOscil
711 G4double WI2 = (adjustmentFactor*adjustm
712 2./3.*(oscStre/totalZ)*Omega*Omega;
713 G4double resEne = std::sqrt(WI2);
714 (*helper)[i].SetResonanceEnergy(resEne);
715 sumLocal += (*helper)[i].GetOscillatorS
716 }
717 }
718 if (sumLocal < TST)
719 AALow = adjustmentFactor;
720 else
721 AAHigh = adjustmentFactor;
722 if (fVerbosityLevel > 3)
723 G4cout << "Sternheimer's adjustment fact
724 adjustmentFactor << " " << TST << " "
725 sumLocal << G4endl;
726 }while((AAHigh-AALow)>(1e-14*adjustmentFacto
727 }
728 else
729 {
730 G4double ionEne = (*helper)[0].GetIonisa
731 (*helper)[0].SetIonisationEnergy(std::fa
732 (*helper)[0].SetResonanceEnergy(meanExci
733 }
734 if (fVerbosityLevel > 1)
735 {
736 G4cout << "Sternheimer's adjustment fact
737 }
738
739 //Check again for data consistency
740 G4double xcheck = (*helper)[0].GetOscillator
741 G4double TST = (*helper)[0].GetOscillatorStr
742 for (std::size_t i=1;i<helper->size();++i)
743 {
744 xcheck += (*helper)[i].GetOscillatorStre
745 TST += (*helper)[i].GetOscillatorStrengt
746 }
747 if (std::fabs(TST-totalZ)>1e-8*totalZ)
748 {
749 G4ExceptionDescription ed;
750 ed << "Inconsistent oscillator data " <<
751 ed << TST << " " << totalZ << G4endl;
752 G4Exception("G4PenelopeOscillatorManager
753 "em2036",FatalException,ed);
754 }
755 xcheck = G4Exp(xcheck/totalZ);
756 if (std::fabs(xcheck-meanExcitationEnergy) >
757 {
758 G4ExceptionDescription ed;
759 ed << "Error in Sterheimer factor calcul
760 ed << xcheck/eV << " " << meanExcitation
761 G4Exception("G4PenelopeOscillatorManager
762 "em2037",FatalException,ed);
763 }
764
765 //Selection of the lowest ionisation energy
766 //ionisation energy less than the N1 shell o
767 //inner shells. As a results, the inner/oute
768 //composition of the material.
769 G4double Zmax = 0;
770 for (G4int k=0;k<nElements;k++)
771 {
772 G4double Z = (*elementVector)[k]->GetZ()
773 if (Z>Zmax) Zmax = Z;
774 }
775 //Find N1 level of the heaviest element (if
776 G4bool found = false;
777 G4double cutEnergy = 50*eV;
778 for (std::size_t i=0;i<helper->size() && !fo
779 {
780 G4double Z = (*helper)[i].GetParentZ();
781 G4int shID = (*helper)[i].GetParentShell
782 if (shID == 10 && Z == Zmax)
783 {
784 found = true;
785 if ((*helper)[i].GetIonisationEnergy() > c
786 cutEnergy = (*helper)[i].GetIonisationEn
787 }
788 }
789 //Make that cutEnergy cannot be higher than
790 //Geant4
791 G4double lowEnergyLimitForFluorescence = 250
792 cutEnergy = std::min(cutEnergy,lowEnergyLimi
793
794 if (fVerbosityLevel > 1)
795 G4cout << "Cutoff energy: " << cutEnergy
796 //
797 //Copy helper in the oscillatorTable for Ion
798 //
799 //Oscillator table Ionisation for the materi
800 G4PenelopeOscillatorTable* theTable = new G4
801 G4PenelopeOscillatorResEnergyComparator comp
802 std::sort(helper->begin(),helper->end(),comp
803
804 //COPY THE HELPER (vector of object) to theT
805 for (std::size_t i=0;i<helper->size();++i)
806 {
807 //copy content --> one may need it later
808 G4PenelopeOscillator* theOsc = new G4Pen
809 theTable->push_back(theOsc);
810 }
811
812 //Oscillators of outer shells with resonance
813 //Rgroup are grouped as a single oscillator
814 G4double Rgroup = 1.05;
815 std::size_t Nost = theTable->size();
816
817 std::size_t firstIndex = (isAConductor) ? 1
818 G4bool loopAgain = false;
819 G4int nLoops = 0;
820 G4int removedLevels = 0;
821 do
822 {
823 loopAgain = false;
824 nLoops++;
825 if (Nost>firstIndex+1)
826 {
827 removedLevels = 0;
828 for (std::size_t i=firstIndex;i<theTable->
829 {
830 G4bool skipLoop = false;
831 G4int shellFlag = (*theTable)[i]->GetS
832 G4double ionEne = (*theTable)[i]->GetI
833 G4double resEne = (*theTable)[i]->GetR
834 G4double resEnePlus1 = (*theTable)[i+1
835 G4double oscStre = (*theTable)[i]->Get
836 G4double oscStrePlus1 = (*theTable)[i+
837 //if (shellFlag < 10 && ionEne>cutEner
838 if (ionEne>cutEnergy) //remove conditi
839 skipLoop = true;
840 if (resEne<1.0*eV || resEnePlus1<1.0*e
841 skipLoop = true;
842 if (resEnePlus1 > Rgroup*resEne)
843 skipLoop = true;
844 if (!skipLoop)
845 {
846 G4double newRes = G4Exp((oscStre*G4Log(r
847 oscStrePlus1*G4Log(resEnePlus1
848 /(oscStre+oscStrePlus1));
849 (*theTable)[i]->SetResonanceEnergy(newRe
850 G4double newIon = (oscStre*ionEne+
851 oscStrePlus1*(*theTable)[i+1]->Ge
852 (oscStre+oscStrePlus1);
853 (*theTable)[i]->SetIonisationEnergy(newI
854 G4double newStre = oscStre+oscStrePlus1;
855 (*theTable)[i]->SetOscillatorStrength(ne
856 G4double newHartree = (oscStre*(*theTabl
857 oscStrePlus1*(*theTable)[i+1]->GetH
858 (oscStre+oscStrePlus1);
859 (*theTable)[i]->SetHartreeFactor(newHart
860 if ((*theTable)[i]->GetParentZ() != (*th
861 (*theTable)[i]->SetParentZ(0.);
862 if (shellFlag < 10 || (*theTable)[i+1]->
863 {
864 G4int newFlag = std::min(shellFlag,(
865 (*theTable)[i]->SetShellFlag(newFlag
866 }
867 else
868 (*theTable)[i]->SetShellFlag(30);
869 //We've lost anyway the track of the ori
870 (*theTable)[i]->SetParentShellID((*theTa
871
872
873 if (i<theTable->size()-2)
874 {
875 for (std::size_t ii=i+1;ii<theTable-
876 (*theTable)[ii] = (*theTable)[ii+1];
877 }
878 //G4cout << theTable->size() << G4endl;
879 theTable->erase(theTable->begin()+theTab
880 removedLevels++;
881 }
882 }
883 }
884 if (removedLevels)
885 {
886 Nost -= removedLevels;
887 loopAgain = true;
888 }
889 if (Rgroup < 1.414213 || Nost > 64)
890 {
891 Rgroup = Rgroup*Rgroup;
892 loopAgain = true;
893 }
894 //Add protection against infinite loops
895 if (nLoops > 100 && !removedLevels)
896 loopAgain = false;
897 }while(loopAgain);
898
899 if (fVerbosityLevel > 1)
900 {
901 G4cout << "Final grouping factor for Ion
902 }
903
904 //Final Electron/Positron model parameters
905 for (std::size_t i=0;i<theTable->size();++i)
906 {
907 //Set cutoff recoil energy for the reson
908 G4double ionEne = (*theTable)[i]->GetIon
909 if (ionEne < 1e-3*eV)
910 {
911 G4double resEne = (*theTable)[i]->GetReson
912 (*theTable)[i]->SetIonisationEnergy(0.*eV)
913 (*theTable)[i]->SetCutoffRecoilResonantEne
914 }
915 else
916 (*theTable)[i]->SetCutoffRecoilResonantEnerg
917 }
918
919 //Last step
920 fOscillatorStoreIonisation->insert(std::make
921
922 /******************************************
923 SAME FOR COMPTON
924 ******************************************/
925 //
926 //Copy helper in the oscillatorTable for Com
927 //
928 //Oscillator table Ionisation for the materi
929 G4PenelopeOscillatorTable* theTableC = new G
930 //order by ionisation energy
931 std::sort(helper->begin(),helper->end());
932 //COPY THE HELPER (vector of object) to theT
933 for (std::size_t i=0;i<helper->size();++i)
934 {
935 //copy content --> one may need it later
936 G4PenelopeOscillator* theOsc = new G4Pen
937 theTableC->push_back(theOsc);
938 }
939 //Oscillators of outer shells with resonance
940 //Rgroup are grouped as a single oscillator
941 Rgroup = 1.5;
942 Nost = theTableC->size();
943
944 firstIndex = (isAConductor) ? 1 : 0; //for c
945 loopAgain = false;
946 removedLevels = 0;
947 do
948 {
949 nLoops++;
950 loopAgain = false;
951 if (Nost>firstIndex+1)
952 {
953 removedLevels = 0;
954 for (std::size_t i=firstIndex;i<theTableC-
955 {
956 G4bool skipLoop = false;
957 G4double ionEne = (*theTableC)[i]->Get
958 G4double ionEnePlus1 = (*theTableC)[i+
959 G4double oscStre = (*theTableC)[i]->Ge
960 G4double oscStrePlus1 = (*theTableC)[i
961 //if (shellFlag < 10 && ionEne>cutEner
962 if (ionEne>cutEnergy)
963 skipLoop = true;
964 if (ionEne<1.0*eV || ionEnePlus1<1.0*e
965 skipLoop = true;
966 if (ionEnePlus1 > Rgroup*ionEne)
967 skipLoop = true;
968
969 if (!skipLoop)
970 {
971 G4double newIon = (oscStre*ionEne+
972 oscStrePlus1*ionEnePlus1)/
973 (oscStre+oscStrePlus1);
974 (*theTableC)[i]->SetIonisationEnergy(new
975 G4double newStre = oscStre+oscStrePlus1;
976 (*theTableC)[i]->SetOscillatorStrength(n
977 G4double newHartree = (oscStre*(*theTabl
978 oscStrePlus1*(*theTableC)[i+1]->Get
979 (oscStre+oscStrePlus1);
980 (*theTableC)[i]->SetHartreeFactor(newHar
981 if ((*theTableC)[i]->GetParentZ() != (*t
982 (*theTableC)[i]->SetParentZ(0.);
983 (*theTableC)[i]->SetShellFlag(30);
984 (*theTableC)[i]->SetParentShellID((*theT
985
986 if (i<theTableC->size()-2)
987 {
988 for (std::size_t ii=i+1;ii<theTableC
989 (*theTableC)[ii] = (*theTableC)[ii+1];
990 }
991 theTableC->erase(theTableC->begin()+theT
992 removedLevels++;
993 }
994 }
995 }
996 if (removedLevels)
997 {
998 Nost -= removedLevels;
999 loopAgain = true;
1000 }
1001 if (Rgroup < 2.0 || Nost > 64)
1002 {
1003 Rgroup = Rgroup*Rgroup;
1004 loopAgain = true;
1005 }
1006 //Add protection against infinite loops
1007 if (nLoops > 100 && !removedLevels)
1008 loopAgain = false;
1009 }while(loopAgain);
1010
1011
1012 if (fVerbosityLevel > 1)
1013 {
1014 G4cout << "Final grouping factor for Co
1015 }
1016
1017 //Last step
1018 fOscillatorStoreCompton->insert(std::make_
1019
1020 //CLEAN UP theHelper and its content
1021 delete helper;
1022 if (fVerbosityLevel > 1)
1023 Dump(material);
1024
1025 return;
1026 }
1027
1028 //....oooOO0OOooo........oooOO0OOooo........o
1029
1030 void G4PenelopeOscillatorManager::ReadElement
1031 {
1032 if (fVerbosityLevel > 0)
1033 {
1034 G4cout << "G4PenelopeOscillatorManager:
1035 G4cout << "Going to read Element Data"
1036 }
1037 const char* path = G4FindDataDir("G4LEDAT
1038 if(!path)
1039 {
1040 G4String excep = "G4PenelopeOscillatorM
1041 G4Exception("G4PenelopeOscillatorManage
1042 "em0006",FatalException,excep);
1043 return;
1044 }
1045 G4String pathString(path);
1046 G4String pathFile = pathString + "/penelope
1047 std::ifstream file(pathFile);
1048
1049 if (!file.is_open())
1050 {
1051 G4String excep = "G4PenelopeOscillatorM
1052 G4Exception("G4PenelopeOscillatorManage
1053 "em0003",FatalException,excep);
1054 }
1055
1056 G4AtomicTransitionManager* theTransitionMan
1057 G4AtomicTransitionManager::Instance();
1058 theTransitionManager->Initialise();
1059
1060 //Read header (22 lines)
1061 G4String theHeader;
1062 for (G4int iline=0;iline<22;iline++)
1063 getline(file,theHeader);
1064 //Done
1065 G4int Z=0;
1066 G4int shellCode = 0;
1067 G4String shellId = "NULL";
1068 G4int occupationNumber = 0;
1069 G4double ionisationEnergy = 0.0*eV;
1070 G4double hartreeProfile = 0.;
1071 G4int shellCounter = 0;
1072 G4int oldZ = -1;
1073 G4int numberOfShells = 0;
1074 //Start reading data
1075 for (G4int i=0;!file.eof();i++)
1076 {
1077 file >> Z >> shellCode >> shellId >> oc
1078 if (Z>0 && i<2000)
1079 {
1080 fElementData[0][i] = Z;
1081 fElementData[1][i] = shellCode;
1082 fElementData[2][i] = occupationNumber;
1083 //reset things
1084 if (Z != oldZ)
1085 {
1086 shellCounter = 0;
1087 oldZ = Z;
1088 numberOfShells = theTransitionManager
1089 }
1090 G4double bindingEnergy = -1*eV;
1091 if (shellCounter<numberOfShells)
1092 {
1093 G4AtomicShell* shell = theTransitionM
1094 bindingEnergy = shell->BindingEnergy(
1095 }
1096 //Valid level found in the G4AtomicTransi
1097 //the ionisation energy found in the Pene
1098 fElementData[3][i] = (bindingEnergy>100*e
1099 fElementData[4][i] = hartreeProfile;
1100 shellCounter++;
1101 }
1102 }
1103 file.close();
1104
1105 if (fVerbosityLevel > 1)
1106 {
1107 G4cout << "G4PenelopeOscillatorManager:
1108 }
1109 fReadElementData = true;
1110 return;
1111 }
1112
1113 //....oooOO0OOooo........oooOO0OOooo........o
1114 G4double G4PenelopeOscillatorManager::GetMean
1115 {
1116 // (1) First time, create fOscillatorStores
1117 CheckForTablesCreated();
1118
1119 // (2) Check if the material has been alrea
1120 if (fExcitationEnergy->count(mat))
1121 return fExcitationEnergy->find(mat)->seco
1122
1123 // (3) If we are here, it means that we hav
1124 BuildOscillatorTable(mat);
1125
1126 // (4) now, the oscillator store should be
1127 if (fExcitationEnergy->count(mat))
1128 return fExcitationEnergy->find(mat)->seco
1129 else
1130 {
1131 G4cout << "G4PenelopeOscillatorManager:
1132 G4cout << "Impossible to retrieve the e
1133 return 0;
1134 }
1135 }
1136
1137 //....oooOO0OOooo........oooOO0OOooo........o
1138 G4double G4PenelopeOscillatorManager::GetPlas
1139 {
1140 // (1) First time, create fOscillatorStores
1141 CheckForTablesCreated();
1142
1143 // (2) Check if the material has been alrea
1144 if (fPlasmaSquared->count(mat))
1145 return fPlasmaSquared->find(mat)->second;
1146
1147 // (3) If we are here, it means that we hav
1148 BuildOscillatorTable(mat);
1149
1150 // (4) now, the oscillator store should be
1151 if (fPlasmaSquared->count(mat))
1152 return fPlasmaSquared->find(mat)->second;
1153 else
1154 {
1155 G4cout << "G4PenelopeOscillatorManager:
1156 G4cout << "Impossible to retrieve the p
1157 return 0;
1158 }
1159 }
1160
1161 //....oooOO0OOooo........oooOO0OOooo........o
1162
1163 G4double G4PenelopeOscillatorManager::GetAtom
1164 {
1165 // (1) First time, create fOscillatorStores
1166 CheckForTablesCreated();
1167
1168 // (2) Check if the material has been alrea
1169 if (fAtomsPerMolecule->count(mat))
1170 return fAtomsPerMolecule->find(mat)->seco
1171
1172 // (3) If we are here, it means that we hav
1173 BuildOscillatorTable(mat);
1174
1175 // (4) now, the oscillator store should be
1176 if (fAtomsPerMolecule->count(mat))
1177 return fAtomsPerMolecule->find(mat)->seco
1178 else
1179 {
1180 G4cout << "G4PenelopeOscillatorManager:
1181 G4cout << "Impossible to retrieve the n
1182 << mat->GetName() << G4endl;
1183 return 0;
1184 }
1185 }
1186
1187 //....oooOO0OOooo........oooOO0OOooo........o
1188
1189 G4double G4PenelopeOscillatorManager::GetNumb
1190 {
1191 // (1) First time, create fOscillatorStores
1192 CheckForTablesCreated();
1193
1194 // (2) Check if the material/Z couple has b
1195 std::pair<const G4Material*,G4int> theKey =
1196 if (fAtomTablePerMolecule->count(theKey))
1197 return fAtomTablePerMolecule->find(theKey
1198
1199 // (3) If we are here, it means that we hav
1200 BuildOscillatorTable(mat);
1201
1202 // (4) now, the oscillator store should be
1203 if (fAtomTablePerMolecule->count(theKey))
1204 return fAtomTablePerMolecule->find(theKey
1205 else
1206 {
1207 G4cout << "G4PenelopeOscillatorManager:
1208 G4cout << "Impossible to retrieve the n
1209 << Z << " in material " << mat->GetNam
1210 return 0;
1211 }
1212 }
1213