Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4LivermorePhotoElectricModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// Author: Sebastien Incerti
28// 30 October 2008
29// on base of G4LowEnergyPhotoElectric developed by A.Forti and M.G.Pia
30//
31// 22 Oct 2012 A & V Ivanchenko Migration data structure to G4PhysicsVector
32// 1 June 2017 M Bandieramonte: New model based on livermore/epics2014
33// evaluated data - parameterization fits in two ranges
34
36
37#include "G4AtomicShell.hh"
38#include "G4AutoLock.hh"
40#include "G4Electron.hh"
41#include "G4Gamma.hh"
42#include "G4LossTableManager.hh"
46#include "G4SystemOfUnits.hh"
48#include "G4EmParameters.hh"
49#include <thread>
50
51//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52
53G4ElementData* G4LivermorePhotoElectricModel::fCrossSection = nullptr;
54G4ElementData* G4LivermorePhotoElectricModel::fCrossSectionLE = nullptr;
55std::vector<G4double>* G4LivermorePhotoElectricModel::fParamHigh[] = {nullptr};
56std::vector<G4double>* G4LivermorePhotoElectricModel::fParamLow[] = {nullptr};
57G4int G4LivermorePhotoElectricModel::fNShells[] = {0};
58G4int G4LivermorePhotoElectricModel::fNShellsUsed[] = {0};
59G4Material* G4LivermorePhotoElectricModel::fWater = nullptr;
60G4double G4LivermorePhotoElectricModel::fWaterEnergyLimit = 0.0;
61G4String G4LivermorePhotoElectricModel::fDataDirectory = "";
62
63static std::once_flag applyOnce;
64
65namespace
66{
67 G4Mutex livPhotoeffMutex = G4MUTEX_INITIALIZER;
68}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
73{
74 verboseLevel = 0;
75 // Verbosity scale:
76 // 0 = nothing
77 // 1 = warning for energy non-conservation
78 // 2 = details of energy budget
79 // 3 = calculation of cross sections, file openings, sampling of atoms
80 // 4 = entering in methods
81
82 theGamma = G4Gamma::Gamma();
83 theElectron = G4Electron::Electron();
84
85 // default generator
87
88 if (verboseLevel > 0) {
89 G4cout << "Livermore PhotoElectric is constructed "
90 << " nShellLimit= " << nShellLimit << G4endl;
91 }
92
93 // Mark this model as "applicable" for atomic deexcitation
95
96 // For water
97 fSandiaCof.resize(4, 0.0);
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
103{
104 if (isInitializer) {
105 for (G4int i = 0; i < ZMAXPE; ++i) {
106 if (fParamHigh[i]) {
107 delete fParamHigh[i];
108 fParamHigh[i] = nullptr;
109 }
110 if (fParamLow[i]) {
111 delete fParamLow[i];
112 fParamLow[i] = nullptr;
113 }
114 }
115 }
116}
117
118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119
121 const G4DataVector&)
122{
123 if (verboseLevel > 1) {
124 G4cout << "Calling G4LivermorePhotoElectricModel::Initialise() " << G4endl;
125 }
126
127 // initialise static tables only once
128 std::call_once(applyOnce, [this]() { isInitializer = true; });
129
130 if (isInitializer) {
131 G4AutoLock l(&livPhotoeffMutex);
132 FindDirectoryPath();
133 if (fWater == nullptr) {
134 fWater = G4Material::GetMaterial("G4_WATER", false);
135 if (fWater == nullptr) {
136 fWater = G4Material::GetMaterial("Water", false);
137 }
138 if (fWater != nullptr) {
139 fWaterEnergyLimit = 13.6 * CLHEP::eV;
140 }
141 }
142
143 if (fCrossSection == nullptr) {
144 fCrossSection = new G4ElementData(ZMAXPE);
145 fCrossSection->SetName("PhotoEffXS");
146 fCrossSectionLE = new G4ElementData(ZMAXPE);
147 fCrossSectionLE->SetName("PhotoEffLowXS");
148 }
149
150 const G4ElementTable* elemTable = G4Element::GetElementTable();
151 std::size_t numElems = (*elemTable).size();
152 for (std::size_t ie = 0; ie < numElems; ++ie) {
153 const G4Element* elem = (*elemTable)[ie];
154 G4int Z = elem->GetZasInt();
155 if (Z < ZMAXPE) {
156 if (fCrossSection->GetElementData(Z) == nullptr) {
157 ReadData(Z);
158 }
159 }
160 }
161 l.unlock();
162 }
163
164 if (verboseLevel > 1) {
165 G4cout << "Loaded cross section files for new LivermorePhotoElectric model" << G4endl;
166 }
167 if (nullptr == fParticleChange) {
169 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
170 }
171
172 fDeexcitationActive = false;
173 if (nullptr != fAtomDeexcitation) {
174 fDeexcitationActive = fAtomDeexcitation->IsFluoActive();
175 }
176
177 if (verboseLevel > 1) {
178 G4cout << "LivermorePhotoElectric model is initialized " << G4endl;
179 }
180}
181
182//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183
186 const G4ParticleDefinition* p,
187 G4double energy,
189{
190 fCurrSection = 0.0;
191 if (fWater && (material == fWater || material->GetBaseMaterial() == fWater)) {
192 if (energy <= fWaterEnergyLimit) {
193 fWater->GetSandiaTable()->GetSandiaCofWater(energy, fSandiaCof);
194
195 G4double energy2 = energy * energy;
196 G4double energy3 = energy * energy2;
197 G4double energy4 = energy2 * energy2;
198
199 fCurrSection = material->GetDensity()
200 * (fSandiaCof[0] / energy + fSandiaCof[1] / energy2 +
201 fSandiaCof[2] / energy3 + fSandiaCof[3] / energy4);
202 }
203 }
204 if (0.0 == fCurrSection) {
205 fCurrSection = G4VEmModel::CrossSectionPerVolume(material, p, energy);
206 }
207 return fCurrSection;
208}
209
210//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211
214 G4double energy,
215 G4double ZZ, G4double,
217{
218 if (verboseLevel > 3) {
219 G4cout << "\n G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom():"
220 << " Z= " << ZZ << " R(keV)= " << energy / keV << G4endl;
221 }
222 G4double cs = 0.0;
223 G4int Z = G4lrint(ZZ);
224 if (Z >= ZMAXPE || Z <= 0) {
225 return cs;
226 }
227 // if element was not initialised
228
229 // do initialisation safely for MT mode
230 if (fCrossSection->GetElementData(Z) == nullptr) {
231 InitialiseOnFly(Z);
232 if (fCrossSection->GetElementData(Z) == nullptr) { return cs; }
233 }
234
235 // 7: rows in the parameterization file; 5: number of parameters
236 G4int idx = fNShells[Z] * 7 - 5;
237
238 energy = std::max(energy, (*(fParamHigh[Z]))[idx - 1]);
239
240 G4double x1 = 1.0 / energy;
241 G4double x2 = x1 * x1;
242 G4double x3 = x2 * x1;
243
244 // high energy parameterisation
245 if (energy >= (*(fParamHigh[Z]))[0]) {
246 G4double x4 = x2 * x2;
247 G4double x5 = x4 * x1;
248
249 cs = x1 *
250 ((*(fParamHigh[Z]))[idx] + x1 * (*(fParamHigh[Z]))[idx + 1]
251 + x2 * (*(fParamHigh[Z]))[idx + 2] + x3 * (*(fParamHigh[Z]))[idx + 3]
252 + x4 * (*(fParamHigh[Z]))[idx + 4] + x5 * (*(fParamHigh[Z]))[idx + 5]);
253 }
254 // low energy parameterisation
255 else if (energy >= (*(fParamLow[Z]))[0]) {
256 G4double x4 = x2 * x2;
257 G4double x5 = x4 * x1; // this variable usage can probably be optimized
258 cs = x1 *
259 ((*(fParamLow[Z]))[idx] + x1 * (*(fParamLow[Z]))[idx + 1]
260 + x2 * (*(fParamLow[Z]))[idx + 2] + x3 * (*(fParamLow[Z]))[idx + 3]
261 + x4 * (*(fParamLow[Z]))[idx + 4] + x5 * (*(fParamLow[Z]))[idx + 5]);
262 }
263 // Tabulated values above k-shell ionization energy
264 else if (energy >= (*(fParamHigh[Z]))[1]) {
265 cs = x3 * (fCrossSection->GetElementData(Z))->Value(energy);
266 }
267 // Tabulated values below k-shell ionization energy
268 else {
269 cs = x3 * (fCrossSectionLE->GetElementData(Z))->Value(energy);
270 }
271 if (verboseLevel > 1) {
272 G4cout << "G4LivermorePhotoElectricModel: E(keV)= " << energy / keV << " Z= " << Z
273 << " cross(barn)= " << cs / barn << G4endl;
274 }
275 return cs;
276}
277
278//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279
280void G4LivermorePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
281 const G4MaterialCutsCouple* couple,
282 const G4DynamicParticle* aDynamicGamma,
284{
285 G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
286 if (verboseLevel > 3) {
287 G4cout << "G4LivermorePhotoElectricModel::SampleSecondaries() Egamma(keV)= "
288 << gammaEnergy / keV << G4endl;
289 }
290
291 // kill incident photon
294
295 // low-energy photo-effect in water - full absorption
296 const G4Material* material = couple->GetMaterial();
297 if (fWater && (material == fWater || material->GetBaseMaterial() == fWater)) {
298 if (gammaEnergy <= fWaterEnergyLimit) {
300 return;
301 }
302 }
303
304 // Returns the normalized direction of the momentum
305 G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
306
307 // Select randomly one element in the current material
308 const G4Element* elm = SelectRandomAtom(material, theGamma, gammaEnergy);
309 G4int Z = elm->GetZasInt();
310
311 // Select the ionised shell in the current atom according to shell
312 // cross sections
313
314 // If element was not initialised gamma should be absorbed
315 if (Z >= ZMAXPE || fCrossSection->GetElementData(Z) == nullptr) {
317 return;
318 }
319
320 // SAMPLING OF THE SHELL INDEX
321 std::size_t shellIdx = 0;
322 std::size_t nn = fNShellsUsed[Z];
323 if (nn > 1) {
324 if (gammaEnergy >= (*(fParamHigh[Z]))[0]) {
325 G4double x1 = 1.0 / gammaEnergy;
326 G4double x2 = x1 * x1;
327 G4double x3 = x2 * x1;
328 G4double x4 = x3 * x1;
329 G4double x5 = x4 * x1;
330 std::size_t idx = nn * 7 - 5;
331 // when do sampling common factors are not taken into account
332 // so cross section is not real
333
334 G4double rand = G4UniformRand();
335 G4double cs0 = rand *
336 ((*(fParamHigh[Z]))[idx] + x1 * (*(fParamHigh[Z]))[idx + 1]
337 + x2 * (*(fParamHigh[Z]))[idx + 2] + x3 * (*(fParamHigh[Z]))[idx + 3]
338 + x4 * (*(fParamHigh[Z]))[idx + 4] + x5 * (*(fParamHigh[Z]))[idx + 5]);
339
340 for (shellIdx = 0; shellIdx < nn; ++shellIdx) {
341 idx = shellIdx * 7 + 2;
342 if (gammaEnergy > (*(fParamHigh[Z]))[idx - 1]) {
343 G4double cs = (*(fParamHigh[Z]))[idx] + x1 * (*(fParamHigh[Z]))[idx + 1]
344 + x2 * (*(fParamHigh[Z]))[idx + 2] + x3 * (*(fParamHigh[Z]))[idx + 3]
345 + x4 * (*(fParamHigh[Z]))[idx + 4] + x5 * (*(fParamHigh[Z]))[idx + 5];
346
347 if (cs >= cs0) {
348 break;
349 }
350 }
351 }
352 if (shellIdx >= nn) {
353 shellIdx = nn - 1;
354 }
355 }
356 else if (gammaEnergy >= (*(fParamLow[Z]))[0]) {
357 G4double x1 = 1.0 / gammaEnergy;
358 G4double x2 = x1 * x1;
359 G4double x3 = x2 * x1;
360 G4double x4 = x3 * x1;
361 G4double x5 = x4 * x1;
362 std::size_t idx = nn * 7 - 5;
363 // when do sampling common factors are not taken into account
364 // so cross section is not real
365 G4double cs0 = G4UniformRand() *
366 ((*(fParamLow[Z]))[idx] + x1 * (*(fParamLow[Z]))[idx + 1]
367 + x2 * (*(fParamLow[Z]))[idx + 2] + x3 * (*(fParamLow[Z]))[idx + 3]
368 + x4 * (*(fParamLow[Z]))[idx + 4] + x5 * (*(fParamLow[Z]))[idx + 5]);
369 for (shellIdx = 0; shellIdx < nn; ++shellIdx) {
370 idx = shellIdx * 7 + 2;
371 if (gammaEnergy > (*(fParamLow[Z]))[idx - 1]) {
372 G4double cs = (*(fParamLow[Z]))[idx] + x1 * (*(fParamLow[Z]))[idx + 1]
373 + x2 * (*(fParamLow[Z]))[idx + 2] + x3 * (*(fParamLow[Z]))[idx + 3]
374 + x4 * (*(fParamLow[Z]))[idx + 4] + x5 * (*(fParamLow[Z]))[idx + 5];
375 if (cs >= cs0) {
376 break;
377 }
378 }
379 }
380 if (shellIdx >= nn) {
381 shellIdx = nn - 1;
382 }
383 }
384 else {
385 // when do sampling common factors are not taken into account
386 // so cross section is not real
387 G4double cs = G4UniformRand();
388
389 if (gammaEnergy >= (*(fParamHigh[Z]))[1]) {
390 // above K-shell binding energy
391 cs *= fCrossSection->GetElementData(Z)->Value(gammaEnergy);
392 }
393 else {
394 // below K-shell binding energy
395 cs *= fCrossSectionLE->GetElementData(Z)->Value(gammaEnergy);
396 }
397
398 for (G4int j = 0; j < (G4int)nn; ++j) {
399 shellIdx = (std::size_t)fCrossSection->GetComponentID(Z, j);
400 if (gammaEnergy > (*(fParamLow[Z]))[7 * shellIdx + 1]) {
401 cs -= fCrossSection->GetValueForComponent(Z, j, gammaEnergy);
402 }
403 if (cs <= 0.0 || j + 1 == (G4int)nn) {
404 break;
405 }
406 }
407 }
408 }
409 // END: SAMPLING OF THE SHELL
410
411 G4double bindingEnergy = (*(fParamHigh[Z]))[shellIdx * 7 + 1];
412 const G4AtomicShell* shell = nullptr;
413
414 // no de-excitation from the last shell
415 if (fDeexcitationActive && shellIdx + 1 < nn) {
416 auto as = G4AtomicShellEnumerator(shellIdx);
417 shell = fAtomDeexcitation->GetAtomicShell(Z, as);
418 }
419
420 // If binding energy of the selected shell is larger than photon energy
421 // do not generate secondaries
422 if (gammaEnergy < bindingEnergy) {
424 return;
425 }
426
427 // Primary outcoming electron
428 G4double eKineticEnergy = gammaEnergy - bindingEnergy;
429 G4double edep = bindingEnergy;
430
431 // Calculate direction of the photoelectron
433 aDynamicGamma, eKineticEnergy, (G4int)shellIdx, couple->GetMaterial());
434
435 // The electron is created
436 auto electron =
437 new G4DynamicParticle(theElectron, electronDirection, eKineticEnergy);
438 fvect->push_back(electron);
439
440 // Sample deexcitation
441 if (nullptr != shell) {
442 G4int index = couple->GetIndex();
443 if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
444 std::size_t nbefore = fvect->size();
445
446 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
447 std::size_t nafter = fvect->size();
448 if (nafter > nbefore) {
449 G4double esec = 0.0;
450 for (std::size_t j = nbefore; j < nafter; ++j) {
451 G4double e = ((*fvect)[j])->GetKineticEnergy();
452 if (esec + e > edep) {
453 // correct energy in order to have energy balance
454 e = edep - esec;
455 ((*fvect)[j])->SetKineticEnergy(e);
456 esec += e;
457 // delete the rest of secondaries (should not happens)
458 for (std::size_t jj = nafter - 1; jj > j; --jj) {
459 delete (*fvect)[jj];
460 fvect->pop_back();
461 }
462 break;
463 }
464 esec += e;
465 }
466 edep -= esec;
467 }
468 }
469 }
470 // energy balance - excitation energy left
471 if (edep > 0.0) {
473 }
474}
475
476//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
477
478const G4String& G4LivermorePhotoElectricModel::FindDirectoryPath()
479{
480 // no check in this method - environment variable is check by utility
481 if (fDataDirectory.empty()) {
482 auto param = G4EmParameters::Instance();
483 std::ostringstream ost;
484 if (param->LivermoreDataDir() == "livermore") {
485 ost << param->GetDirLEDATA() << "/livermore/phot_epics2014/";
486 }
487 else {
488 ost << param->GetDirLEDATA() << "/epics2017/phot/";
489 }
490 fDataDirectory = ost.str();
491 }
492 return fDataDirectory;
493}
494
495//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
496
497void G4LivermorePhotoElectricModel::ReadData(G4int Z)
498{
499 if (Z <= 0 || Z >= ZMAXPE) {
500 G4cout << "G4LivermorePhotoElectricModel::ReadData() Warning: Z="
501 << Z << " is out of range - request ignored" << G4endl;
502 return;
503 }
504 if (verboseLevel > 1) {
505 G4cout << "G4LivermorePhotoElectricModel::ReadData() for Z=" << Z << G4endl;
506 }
507
508 if (fCrossSection->GetElementData(Z) != nullptr) {
509 return;
510 }
511
512 // spline for photoeffect total x-section above K-shell when using EPDL97
513 // but below the parameterized ones
514
515 G4bool spline = (G4EmParameters::Instance()->LivermoreDataDir() == "livermore");
517 auto pv = new G4PhysicsFreeVector(spline);
518
519 // fDataDirectory will be defined after these lines
520 std::ostringstream ost;
521 ost << FindDirectoryPath() << "pe-cs-" << Z << ".dat";
522 std::ifstream fin(ost.str().c_str());
523 if (!fin.is_open()) {
525 ed << "G4LivermorePhotoElectricModel data file <"
526 << ost.str().c_str() << "> is not opened!" << G4endl;
527 G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
528 "G4LEDATA version should be G4EMLOW8.0 or later.");
529 return;
530 }
531 if (verboseLevel > 3) {
532 G4cout << "File " << ost.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
533 << G4endl;
534 }
535 pv->Retrieve(fin, true);
536 pv->ScaleVector(MeV, barn);
537 pv->FillSecondDerivatives();
538 pv->EnableLogBinSearch(number);
539 fCrossSection->InitialiseForElement(Z, pv);
540 fin.close();
541
542 // read high-energy fit parameters
543 fParamHigh[Z] = new std::vector<G4double>;
544 G4int n1 = 0;
545 G4int n2 = 0;
546 G4double x;
547 std::ostringstream ost1;
548 ost1 << fDataDirectory << "pe-high-" << Z << ".dat";
549 std::ifstream fin1(ost1.str().c_str());
550 if (!fin1.is_open()) {
552 ed << "G4LivermorePhotoElectricModel data file <"
553 << ost1.str().c_str() << "> is not opened!" << G4endl;
554 G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
555 "G4LEDATA version should be G4EMLOW7.2 or later.");
556 return;
557 }
558 if (verboseLevel > 3) {
559 G4cout << "File " << ost1.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
560 << G4endl;
561 }
562 fin1 >> n1;
563 if (fin1.fail()) {
564 return;
565 }
566 if (0 > n1 || n1 >= INT_MAX) {
567 n1 = 0;
568 }
569
570 fin1 >> n2;
571 if (fin1.fail()) {
572 return;
573 }
574 if (0 > n2 || n2 >= INT_MAX) {
575 n2 = 0;
576 }
577
578 fin1 >> x;
579 if (fin1.fail()) {
580 return;
581 }
582
583 fNShells[Z] = n1;
584 fParamHigh[Z]->reserve(7 * n1 + 1);
585 fParamHigh[Z]->push_back(x * MeV);
586 for (G4int i = 0; i < n1; ++i) {
587 for (G4int j = 0; j < 7; ++j) {
588 fin1 >> x;
589 if (0 == j) {
590 x *= MeV;
591 }
592 else {
593 x *= barn;
594 }
595 fParamHigh[Z]->push_back(x);
596 }
597 }
598 fin1.close();
599
600 // read low-energy fit parameters
601 fParamLow[Z] = new std::vector<G4double>;
602 G4int n1_low = 0;
603 G4int n2_low = 0;
604 G4double x_low;
605 std::ostringstream ost1_low;
606 ost1_low << fDataDirectory << "pe-low-" << Z << ".dat";
607 std::ifstream fin1_low(ost1_low.str().c_str());
608 if (!fin1_low.is_open()) {
610 ed << "G4LivermorePhotoElectricModel data file <" << ost1_low.str().c_str()
611 << "> is not opened!" << G4endl;
612 G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
613 "G4LEDATA version should be G4EMLOW8.0 or later.");
614 return;
615 }
616 if (verboseLevel > 3) {
617 G4cout << "File " << ost1_low.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
618 << G4endl;
619 }
620 fin1_low >> n1_low;
621 if (fin1_low.fail()) {
622 return;
623 }
624 if (0 > n1_low || n1_low >= INT_MAX) {
625 n1_low = 0;
626 }
627
628 fin1_low >> n2_low;
629 if (fin1_low.fail()) {
630 return;
631 }
632 if (0 > n2_low || n2_low >= INT_MAX) {
633 n2_low = 0;
634 }
635
636 fin1_low >> x_low;
637 if (fin1_low.fail()) {
638 return;
639 }
640
641 fNShells[Z] = n1_low;
642 fParamLow[Z]->reserve(7 * n1_low + 1);
643 fParamLow[Z]->push_back(x_low * MeV);
644 for (G4int i = 0; i < n1_low; ++i) {
645 for (G4int j = 0; j < 7; ++j) {
646 fin1_low >> x_low;
647 if (0 == j) {
648 x_low *= MeV;
649 }
650 else {
651 x_low *= barn;
652 }
653 fParamLow[Z]->push_back(x_low);
654 }
655 }
656 fin1_low.close();
657
658 // there is a possibility to use only main shells
659 if (nShellLimit < n2) {
660 n2 = nShellLimit;
661 }
662 fCrossSection->InitialiseForComponent(Z, n2); // number of shells
663 fNShellsUsed[Z] = n2;
664
665 if (1 < n2) {
666 std::ostringstream ost2;
667 ost2 << fDataDirectory << "pe-ss-cs-" << Z << ".dat";
668 std::ifstream fin2(ost2.str().c_str());
669 if (!fin2.is_open()) {
671 ed << "G4LivermorePhotoElectricModel data file <" << ost2.str().c_str() << "> is not opened!"
672 << G4endl;
673 G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
674 "G4LEDATA version should be G4EMLOW7.2 or later.");
675 return;
676 }
677 if (verboseLevel > 3) {
678 G4cout << "File " << ost2.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
679 << G4endl;
680 }
681
682 G4int n3, n4;
683 G4double y;
684 for (G4int i = 0; i < n2; ++i) {
685 fin2 >> x >> y >> n3 >> n4;
686 auto v = new G4PhysicsFreeVector(n3, x, y);
687 for (G4int j = 0; j < n3; ++j) {
688 fin2 >> x >> y;
689 v->PutValues(j, x * CLHEP::MeV, y * CLHEP::barn);
690 }
691 v->EnableLogBinSearch(number);
692 fCrossSection->AddComponent(Z, n4, v);
693 }
694 fin2.close();
695 }
696
697 // no spline for photoeffect total x-section below K-shell
698 if (1 < fNShells[Z]) {
699 auto pv1 = new G4PhysicsFreeVector(false);
700 std::ostringstream ost3;
701 ost3 << fDataDirectory << "pe-le-cs-" << Z << ".dat";
702 std::ifstream fin3(ost3.str().c_str());
703 if (!fin3.is_open()) {
705 ed << "G4LivermorePhotoElectricModel data file <" << ost3.str().c_str() << "> is not opened!"
706 << G4endl;
707 G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
708 "G4LEDATA version should be G4EMLOW8.0 or later.");
709 return;
710 }
711 if (verboseLevel > 3) {
712 G4cout << "File " << ost3.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
713 << G4endl;
714 }
715 pv1->Retrieve(fin3, true);
716 pv1->ScaleVector(CLHEP::MeV, CLHEP::barn);
717 pv1->EnableLogBinSearch(number);
718 fCrossSectionLE->InitialiseForElement(Z, pv1);
719 fin3.close();
720 }
721}
722
723//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
724
726{
727 if (Z < 1 || Z >= ZMAXPE) {
728 return -1;
729 } // If Z is out of the supported return 0
730
731 InitialiseOnFly(Z);
732 if (fCrossSection->GetElementData(Z) == nullptr ||
733 shell < 0 || shell >= fNShellsUsed[Z]) {
734 return -1;
735 }
736
737 if (Z > 2) {
738 return fCrossSection->GetComponentDataByIndex(Z, shell)->Energy(0);
739 }
740 else {
741 return fCrossSection->GetElementData(Z)->Energy(0);
742 }
743}
744
745//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
746
747void
749{
750 if (fCrossSection == nullptr) {
751 fCrossSection = new G4ElementData(ZMAXPE);
752 fCrossSection->SetName("PhotoEffXS");
753 fCrossSectionLE = new G4ElementData(ZMAXPE);
754 fCrossSectionLE->SetName("PhotoEffLowXS");
755 }
756 ReadData(Z);
757}
758
759//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
760
761void G4LivermorePhotoElectricModel::InitialiseOnFly(G4int Z)
762{
763 if (fCrossSection->GetElementData(Z) == nullptr && Z > 0 && Z < ZMAXPE) {
764 G4AutoLock l(&livPhotoeffMutex);
765 if (fCrossSection->GetElementData(Z) == nullptr) {
766 ReadData(Z);
767 }
768 l.unlock();
769 }
770}
771
772//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
std::vector< G4Element * > G4ElementTable
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
#define elem(i, j)
#define G4MUTEX_INITIALIZER
std::mutex G4Mutex
@ fStopAndKill
double G4double
Definition G4Types.hh:83
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition Randomize.hh:52
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition G4Electron.cc:91
G4PhysicsVector * GetElementData(G4int Z) const
void InitialiseForComponent(G4int Z, G4int nComponents=0)
G4PhysicsVector * GetComponentDataByIndex(G4int Z, std::size_t idx) const
void InitialiseForElement(G4int Z, G4PhysicsVector *v)
G4int GetComponentID(G4int Z, std::size_t idx) const
G4double GetValueForComponent(G4int Z, std::size_t idx, G4double kinEnergy) const
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
void SetName(const G4String &nam)
static G4ElementTable * GetElementTable()
Definition G4Element.cc:389
G4int GetZasInt() const
Definition G4Element.hh:120
static G4EmParameters * Instance()
const G4String & LivermoreDataDir()
G4int NumberForFreeVector() const
static G4Gamma * Gamma()
Definition G4Gamma.cc:81
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double energy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) override
G4double GetBindingEnergy(G4int Z, G4int shell)
void InitialiseForElement(const G4ParticleDefinition *, G4int Z) override
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4LivermorePhotoElectricModel(const G4String &nam="LivermorePhElectric")
G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double energy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX) override
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
G4double GetDensity() const
const G4Material * GetBaseMaterial() const
G4SandiaTable * GetSandiaTable() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double Energy(const std::size_t index) const
G4double Value(const G4double energy, std::size_t &lastidx) const
void GetSandiaCofWater(G4double energy, std::vector< G4double > &coeff) const
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
G4VEmAngularDistribution * GetAngularDistribution()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void SetDeexcitationFlag(G4bool val)
void SetAngularDistribution(G4VEmAngularDistribution *)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
int G4lrint(double ad)
Definition templates.hh:134
#define INT_MAX
Definition templates.hh:90