Geant4 10.7.0
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#include "G4SystemOfUnits.hh"
38#include "G4LossTableManager.hh"
39#include "G4Electron.hh"
40#include "G4Gamma.hh"
46#include "G4AtomicShell.hh"
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50G4LPhysicsFreeVector* G4LivermorePhotoElectricModel::fCrossSection[] = {nullptr};
51G4LPhysicsFreeVector* G4LivermorePhotoElectricModel::fCrossSectionLE[] = {nullptr};
52std::vector<G4double>* G4LivermorePhotoElectricModel::fParamHigh[] = {nullptr};
53std::vector<G4double>* G4LivermorePhotoElectricModel::fParamLow[] = {nullptr};
54G4int G4LivermorePhotoElectricModel::fNShells[] = {0};
55G4int G4LivermorePhotoElectricModel::fNShellsUsed[] = {0};
56G4ElementData* G4LivermorePhotoElectricModel::fShellCrossSection = nullptr;
57G4Material* G4LivermorePhotoElectricModel::fWater = nullptr;
58G4double G4LivermorePhotoElectricModel::fWaterEnergyLimit = 0.0;
59G4String G4LivermorePhotoElectricModel::fDataDirectory = "";
60
61#ifdef G4MULTITHREADED
62 G4Mutex G4LivermorePhotoElectricModel::livPhotoeffMutex = G4MUTEX_INITIALIZER;
63#endif
64
65using namespace std;
66
67//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68
70 : G4VEmModel(nam),fParticleChange(nullptr),maxZ(99),
71 nShellLimit(100),fDeexcitationActive(false),isInitialised(false),
72 fAtomDeexcitation(nullptr)
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 fSandiaCof.resize(4,0.0);
96 fCurrSection = 0.0;
97}
98
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100
102{
103 if(IsMaster()) {
104 delete fShellCrossSection;
105 fShellCrossSection = nullptr;
106 for(G4int i=0; i<maxZ; ++i) {
107 delete fParamHigh[i];
108 fParamHigh[i] = nullptr;
109 delete fParamLow[i];
110 fParamLow[i] = nullptr;
111 delete fCrossSection[i];
112 fCrossSection[i] = nullptr;
113 delete fCrossSectionLE[i];
114 fCrossSectionLE[i] = nullptr;
115 }
116 }
117}
118
119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120
121void
123 const G4DataVector&)
124{
125 if (verboseLevel > 2) {
126 G4cout << "Calling G4LivermorePhotoElectricModel::Initialise() " << G4endl;
127 }
128
129 if(IsMaster()) {
130
131 if(!fWater) {
132 fWater = G4Material::GetMaterial("G4_WATER", false);
133 if(!fWater) { fWater = G4Material::GetMaterial("Water", false); }
134 if(fWater) { fWaterEnergyLimit = 13.6*eV; }
135 }
136
137 if(!fShellCrossSection) { fShellCrossSection = new G4ElementData(); }
138
139 G4ProductionCutsTable* theCoupleTable =
141 G4int numOfCouples = theCoupleTable->GetTableSize();
142
143 for(G4int i=0; i<numOfCouples; ++i) {
144 const G4MaterialCutsCouple* couple =
145 theCoupleTable->GetMaterialCutsCouple(i);
146 const G4Material* material = couple->GetMaterial();
147 const G4ElementVector* theElementVector = material->GetElementVector();
148 G4int nelm = material->GetNumberOfElements();
149
150 for (G4int j=0; j<nelm; ++j) {
151 G4int Z = std::min(maxZ, (*theElementVector)[j]->GetZasInt());
152 if(!fCrossSection[Z]) { ReadData(Z); }
153 }
154 }
155 }
156
157 if (verboseLevel > 2) {
158 G4cout << "Loaded cross section files for new LivermorePhotoElectric model"
159 << G4endl;
160 }
161 if(!isInitialised) {
162 isInitialised = true;
164
165 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
166 }
167 fDeexcitationActive = false;
168 if(fAtomDeexcitation) {
169 fDeexcitationActive = fAtomDeexcitation->IsFluoActive();
170 }
171
172 if (verboseLevel > 0) {
173 G4cout << "LivermorePhotoElectric model is initialized " << G4endl
174 << G4endl;
175 }
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
181 const G4Material* material,
182 const G4ParticleDefinition* p,
183 G4double energy,
185{
186 fCurrSection = 0.0;
187 if(fWater && (material == fWater ||
188 material->GetBaseMaterial() == fWater)) {
189 if(energy <= fWaterEnergyLimit) {
190 fWater->GetSandiaTable()->GetSandiaCofWater(energy, fSandiaCof);
191
192 G4double energy2 = energy*energy;
193 G4double energy3 = energy*energy2;
194 G4double energy4 = energy2*energy2;
195
196 fCurrSection = material->GetDensity()*
197 (fSandiaCof[0]/energy + fSandiaCof[1]/energy2 +
198 fSandiaCof[2]/energy3 + fSandiaCof[3]/energy4);
199 }
200 }
201 if(0.0 == fCurrSection) {
202 fCurrSection = G4VEmModel::CrossSectionPerVolume(material, p, energy);
203 }
204 return fCurrSection;
205}
206
207//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
208
211 G4double energy,
212 G4double ZZ, G4double,
214{
215 if (verboseLevel > 3) {
216 G4cout << "\n G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom():"
217 << " Z= " << ZZ << " R(keV)= " << energy/keV << G4endl;
218 }
219 G4double cs = 0.0;
220 G4int Z = G4lrint(ZZ);
221 if(Z >= maxZ) { return cs; }
222 // if element was not initialised
223
224 // do initialisation safely for MT mode
225 if(!fCrossSection[Z]) { InitialiseForElement(theGamma, Z); }
226
227 //7: rows in the parameterization file; 5: number of parameters
228 G4int idx = fNShells[Z]*7 - 5;
229
230 energy = std::max(energy, (*(fParamHigh[Z]))[idx-1]);
231
232 G4double x1 = 1.0/energy;
233 G4double x2 = x1*x1;
234 G4double x3 = x2*x1;
235
236 // high energy parameterisation
237 if(energy >= (*(fParamHigh[Z]))[0]) {
238
239 G4double x4 = x2*x2;
240 G4double x5 = x4*x1;
241
242 cs = x1*((*(fParamHigh[Z]))[idx] + x1*(*(fParamHigh[Z]))[idx+1]
243 + x2*(*(fParamHigh[Z]))[idx+2] + x3*(*(fParamHigh[Z]))[idx+3]
244 + x4*(*(fParamHigh[Z]))[idx+4]+ x5*(*(fParamHigh[Z]))[idx+5]);
245
246 }
247 // low energy parameterisation
248 else if(energy >= (*(fParamLow[Z]))[0]) {
249
250 G4double x4 = x2*x2;
251 G4double x5 = x4*x1; //this variable usage can probably be optimized
252 cs = x1*((*(fParamLow[Z]))[idx] + x1*(*(fParamLow[Z]))[idx+1]
253 + x2*(*(fParamLow[Z]))[idx+2] + x3*(*(fParamLow[Z]))[idx+3]
254 + x4*(*(fParamLow[Z]))[idx+4]+ x5*(*(fParamLow[Z]))[idx+5]);
255
256 }
257 // Tabulated values above k-shell ionization energy
258 else if(energy >= (*(fParamHigh[Z]))[1]) {
259 cs = x3*(fCrossSection[Z])->Value(energy);
260 }
261 // Tabulated values below k-shell ionization energy
262 else
263 {
264 cs = x3*(fCrossSectionLE[Z])->Value(energy);
265 }
266 if (verboseLevel > 1) {
267 G4cout << "G4LivermorePhotoElectricModel: E(keV)= " << energy/keV
268 << " Z= " << Z << " cross(barn)= " << cs/barn << G4endl;
269 }
270 return cs;
271}
272
273//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
274
275void
277 std::vector<G4DynamicParticle*>* fvect,
278 const G4MaterialCutsCouple* couple,
279 const G4DynamicParticle* aDynamicGamma,
281{
282 G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
283 if (verboseLevel > 3) {
284 G4cout << "G4LivermorePhotoElectricModel::SampleSecondaries() Egamma(keV)= "
285 << gammaEnergy/keV << G4endl;
286 }
287
288 // kill incident photon
291
292 // low-energy photo-effect in water - full absorption
293 const G4Material* material = couple->GetMaterial();
294 if(fWater && (material == fWater ||
295 material->GetBaseMaterial() == fWater)) {
296 if(gammaEnergy <= fWaterEnergyLimit) {
298 return;
299 }
300 }
301
302 // Returns the normalized direction of the momentum
303 G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
304
305 // Select randomly one element in the current material
306 //G4cout << "Select random atom Egamma(keV)= " << gammaEnergy/keV << G4endl;
307 const G4Element* elm = SelectRandomAtom(material, theGamma, gammaEnergy);
308 G4int Z = elm->GetZasInt();
309
310 // Select the ionised shell in the current atom according to shell
311 // cross sections
312 // G4cout << "Select random shell Z= " << Z << G4endl;
313
314 if(Z >= maxZ) { Z = maxZ-1; }
315
316 // element was not initialised gamma should be absorbed
317 if(!fCrossSection[Z]) {
319 return;
320 }
321
322 // SAMPLING OF THE SHELL INDEX
323 size_t shellIdx = 0;
324 size_t nn = fNShellsUsed[Z];
325 if(nn > 1)
326 {
327 if(gammaEnergy >= (*(fParamHigh[Z]))[0])
328 {
329 G4double x1 = 1.0/gammaEnergy;
330 G4double x2 = x1*x1;
331 G4double x3 = x2*x1;
332 G4double x4 = x3*x1;
333 G4double x5 = x4*x1;
334 G4int idx = nn*7 - 5;
335 // when do sampling common factors are not taken into account
336 // so cross section is not real
337
338 G4double rand=G4UniformRand();
339 G4double cs0 = rand*( (*(fParamHigh[Z]))[idx]
340 + x1*(*(fParamHigh[Z]))[idx+1]
341 + x2*(*(fParamHigh[Z]))[idx+2]
342 + x3*(*(fParamHigh[Z]))[idx+3]
343 + x4*(*(fParamHigh[Z]))[idx+4]
344 + x5*(*(fParamHigh[Z]))[idx+5]);
345
346 for(shellIdx=0; shellIdx<nn; ++shellIdx)
347 {
348 idx = shellIdx*7 + 2;
349 if(gammaEnergy > (*(fParamHigh[Z]))[idx-1])
350 {
351 G4double cs =
352 (*(fParamHigh[Z]))[idx]
353 + x1*(*(fParamHigh[Z]))[idx+1]
354 + x2*(*(fParamHigh[Z]))[idx+2]
355 + x3*(*(fParamHigh[Z]))[idx+3]
356 + x4*(*(fParamHigh[Z]))[idx+4]
357 + x5*(*(fParamHigh[Z]))[idx+5];
358
359 if(cs >= cs0) { break; }
360 }
361 }
362 if(shellIdx >= nn) { shellIdx = nn-1; }
363 }
364 else if(gammaEnergy >= (*(fParamLow[Z]))[0])
365 {
366 G4double x1 = 1.0/gammaEnergy;
367 G4double x2 = x1*x1;
368 G4double x3 = x2*x1;
369 G4double x4 = x3*x1;
370 G4double x5 = x4*x1;
371 G4int idx = nn*7 - 5;
372 // when do sampling common factors are not taken into account
373 // so cross section is not real
374 G4double cs0 = G4UniformRand()*((*(fParamLow[Z]))[idx]
375 + x1*(*(fParamLow[Z]))[idx+1]
376 + x2*(*(fParamLow[Z]))[idx+2]
377 + x3*(*(fParamLow[Z]))[idx+3]
378 + x4*(*(fParamLow[Z]))[idx+4]
379 + x5*(*(fParamLow[Z]))[idx+5]);
380 for(shellIdx=0; shellIdx<nn; ++shellIdx)
381 {
382 idx = shellIdx*7 + 2;
383 if(gammaEnergy > (*(fParamLow[Z]))[idx-1])
384 {
385 G4double cs = (*(fParamLow[Z]))[idx] + x1*(*(fParamLow[Z]))[idx+1]
386 + x2*(*(fParamLow[Z]))[idx+2] + x3*(*(fParamLow[Z]))[idx+3]
387 + x4*(*(fParamLow[Z]))[idx+4]+ x5*(*(fParamLow[Z]))[idx+5];
388 if(cs >= cs0) { break; }
389 }
390 }
391 if(shellIdx >= nn) {shellIdx = nn-1;}
392 }
393 else
394 {
395 // when do sampling common factors are not taken into account
396 // so cross section is not real
397 G4double cs = G4UniformRand();
398
399 if(gammaEnergy >= (*(fParamHigh[Z]))[1]) {
400 //above K-shell binding energy
401 cs*= (fCrossSection[Z])->Value(gammaEnergy);
402 }
403 else
404 {
405 //below K-shell binding energy
406 cs *= (fCrossSectionLE[Z])->Value(gammaEnergy);
407 }
408
409 for(size_t j=0; j<nn; ++j) {
410
411 shellIdx = (size_t)fShellCrossSection->GetComponentID(Z, j);
412 if(gammaEnergy > (*(fParamLow[Z]))[7*shellIdx+1]) {
413 cs -= fShellCrossSection->GetValueForComponent(Z, j, gammaEnergy);
414 }
415 if(cs <= 0.0 || j+1 == nn) {break;}
416 }
417 }
418 }
419 // END: SAMPLING OF THE SHELL
420
421 G4double bindingEnergy = (*(fParamHigh[Z]))[shellIdx*7 + 1];
422 //G4cout << "Z= " << Z << " shellIdx= " << shellIdx
423 // << " nShells= " << fNShells[Z]
424 // << " Ebind(keV)= " << bindingEnergy/keV
425 // << " Egamma(keV)= " << gammaEnergy/keV << G4endl;
426
427 const G4AtomicShell* shell = 0;
428
429 // no de-excitation from the last shell
430 if(fDeexcitationActive && shellIdx + 1 < nn) {
432 shell = fAtomDeexcitation->GetAtomicShell(Z, as);
433 }
434
435 // If binding energy of the selected shell is larger than photon energy
436 // do not generate secondaries
437 if(gammaEnergy < bindingEnergy) {
439 return;
440 }
441
442 // Primary outcoming electron
443 G4double eKineticEnergy = gammaEnergy - bindingEnergy;
444 G4double edep = bindingEnergy;
445
446 // Calculate direction of the photoelectron
447 G4ThreeVector electronDirection =
449 eKineticEnergy,
450 shellIdx,
451 couple->GetMaterial());
452
453 // The electron is created
454 G4DynamicParticle* electron = new G4DynamicParticle (theElectron,
455 electronDirection,
456 eKineticEnergy);
457 fvect->push_back(electron);
458
459 // Sample deexcitation
460 if(shell) {
461 G4int index = couple->GetIndex();
462 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
463 G4int nbefore = fvect->size();
464
465 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
466 G4int nafter = fvect->size();
467 if(nafter > nbefore) {
468 G4double esec = 0.0;
469 for (G4int j=nbefore; j<nafter; ++j) {
470
471 G4double e = ((*fvect)[j])->GetKineticEnergy();
472 if(esec + e > edep) {
473 // correct energy in order to have energy balance
474 e = edep - esec;
475 ((*fvect)[j])->SetKineticEnergy(e);
476 esec += e;
477 // delete the rest of secondaries (should not happens)
478 for (G4int jj=nafter-1; jj>j; --jj) {
479 delete (*fvect)[jj];
480 fvect->pop_back();
481 }
482 break;
483 }
484 esec += e;
485 }
486 edep -= esec;
487 }
488 }
489 }
490 // energy balance - excitation energy left
491 if(edep > 0.0) {
493 }
494}
495
496//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
497
498const G4String& G4LivermorePhotoElectricModel::FindDirectoryPath()
499{
500 // check environment variable
501 // build the complete string identifying the file with the data set
502 if(fDataDirectory.empty()) {
503 const char* path = std::getenv("G4LEDATA");
504 if (path) {
505 std::ostringstream ost;
506 ost << path << "/livermore/phot_epics2014/";
507 fDataDirectory = ost.str();
508 } else {
509 G4Exception("G4SeltzerBergerModel::FindDirectoryPath()","em0006",
511 "Environment variable G4LEDATA not defined");
512 }
513 }
514 return fDataDirectory;
515}
516
517//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
518
519void G4LivermorePhotoElectricModel::ReadData(G4int Z)
520{
521 if (verboseLevel > 1)
522 {
523 G4cout << "Calling ReadData() of G4LivermorePhotoElectricModel"
524 << G4endl;
525 }
526
527 if(fCrossSection[Z]) { return; }
528
529 // spline for photoeffect total x-section above K-shell
530 // but below the parameterized ones
531 fCrossSection[Z] = new G4LPhysicsFreeVector();
532 fCrossSection[Z]->SetSpline(true);
533 std::ostringstream ost;
534 ost << FindDirectoryPath() << "pe-cs-" << Z <<".dat";
535 std::ifstream fin(ost.str().c_str());
536 if( !fin.is_open()) {
538 ed << "G4LivermorePhotoElectricModel data file <" << ost.str().c_str()
539 << "> is not opened!" << G4endl;
540 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
541 "em0003",FatalException,
542 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
543 return;
544 }
545 if(verboseLevel > 3) {
546 G4cout << "File " << ost.str().c_str()
547 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
548 }
549 fCrossSection[Z]->Retrieve(fin, true);
550 fCrossSection[Z]->ScaleVector(MeV, barn);
551 fin.close();
552
553 // read high-energy fit parameters
554 fParamHigh[Z] = new std::vector<G4double>;
555 G4int n1 = 0;
556 G4int n2 = 0;
557 G4double x;
558 std::ostringstream ost1;
559 ost1 << fDataDirectory << "pe-high-" << Z <<".dat";
560 std::ifstream fin1(ost1.str().c_str());
561 if( !fin1.is_open()) {
563 ed << "G4LivermorePhotoElectricModel data file <" << ost1.str().c_str()
564 << "> is not opened!" << G4endl;
565 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
566 "em0003",FatalException,
567 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
568 return;
569 }
570 if(verboseLevel > 3) {
571 G4cout << "File " << ost1.str().c_str()
572 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
573 }
574 fin1 >> n1;
575 if(fin1.fail()) { return; }
576 if(0 > n1 || n1 >= INT_MAX) { n1 = 0; }
577
578 fin1 >> n2;
579 if(fin1.fail()) { return; }
580 if(0 > n2 || n2 >= INT_MAX) { n2 = 0; }
581
582 fin1 >> x;
583 if(fin1.fail()) { return; }
584
585 fNShells[Z] = n1;
586 fParamHigh[Z]->reserve(7*n1+1);
587 fParamHigh[Z]->push_back(x*MeV);
588 for(G4int i=0; i<n1; ++i) {
589 for(G4int j=0; j<7; ++j) {
590 fin1 >> x;
591 if(0 == j) { x *= MeV; }
592 else { x *= barn; }
593 fParamHigh[Z]->push_back(x);
594 }
595 }
596 fin1.close();
597
598 // read low-energy fit parameters
599 fParamLow[Z] = new std::vector<G4double>;
600 G4int n1_low = 0;
601 G4int n2_low = 0;
602 G4double x_low;
603 std::ostringstream ost1_low;
604 ost1_low << fDataDirectory << "pe-low-" << Z <<".dat";
605 std::ifstream fin1_low(ost1_low.str().c_str());
606 if( !fin1_low.is_open()) {
608 ed << "G4LivermorePhotoElectricModel data file <" << ost1_low.str().c_str()
609 << "> is not opened!" << G4endl;
610 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
611 "em0003",FatalException,
612 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
613 return;
614 }
615 if(verboseLevel > 3) {
616 G4cout << "File " << ost1_low.str().c_str()
617 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
618 }
619 fin1_low >> n1_low;
620 if(fin1_low.fail()) { return; }
621 if(0 > n1_low || n1_low >= INT_MAX) { n1_low = 0; }
622
623 fin1_low >> n2_low;
624 if(fin1_low.fail()) { return; }
625 if(0 > n2_low || n2_low >= INT_MAX) { n2_low = 0; }
626
627 fin1_low >> x_low;
628 if(fin1_low.fail()) { return; }
629
630 fNShells[Z] = n1_low;
631 fParamLow[Z]->reserve(7*n1_low+1);
632 fParamLow[Z]->push_back(x_low*MeV);
633 for(G4int i=0; i<n1_low; ++i) {
634 for(G4int j=0; j<7; ++j) {
635 fin1_low >> x_low;
636 if(0 == j) { x_low *= MeV; }
637 else { x_low *= barn; }
638 fParamLow[Z]->push_back(x_low);
639 }
640 }
641 fin1_low.close();
642
643 // there is a possibility to use only main shells
644 if(nShellLimit < n2) { n2 = nShellLimit; }
645 fShellCrossSection->InitialiseForComponent(Z, n2);//number of shells
646 fNShellsUsed[Z] = n2;
647
648 if(1 < n2) {
649 std::ostringstream ost2;
650 ost2 << fDataDirectory << "pe-ss-cs-" << Z <<".dat";
651 std::ifstream fin2(ost2.str().c_str());
652 if( !fin2.is_open()) {
654 ed << "G4LivermorePhotoElectricModel data file <" << ost2.str().c_str()
655 << "> is not opened!" << G4endl;
656 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
657 "em0003",FatalException,
658 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
659 return;
660 }
661 if(verboseLevel > 3) {
662 G4cout << "File " << ost2.str().c_str()
663 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
664 }
665
666 G4int n3, n4;
667 G4double y;
668 for(G4int i=0; i<n2; ++i) {
669 fin2 >> x >> y >> n3 >> n4;
671 for(G4int j=0; j<n3; ++j) {
672 fin2 >> x >> y;
673 v->PutValues(j, x*MeV, y*barn);
674 }
675 fShellCrossSection->AddComponent(Z, n4, v);
676 }
677 fin2.close();
678 }
679
680 // no spline for photoeffect total x-section below K-shell
681 if(1 < fNShells[Z]) {
682 fCrossSectionLE[Z] = new G4LPhysicsFreeVector();
683 std::ostringstream ost3;
684 ost3 << fDataDirectory << "pe-le-cs-" << Z <<".dat";
685 std::ifstream fin3(ost3.str().c_str());
686 if( !fin3.is_open()) {
688 ed << "G4LivermorePhotoElectricModel data file <" << ost3.str().c_str()
689 << "> is not opened!" << G4endl;
690 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
691 "em0003",FatalException,
692 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
693 return;
694 }
695 if(verboseLevel > 3) {
696 G4cout << "File " << ost3.str().c_str()
697 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
698 }
699 fCrossSectionLE[Z]->Retrieve(fin3, true);
700 fCrossSectionLE[Z]->ScaleVector(MeV, barn);
701 fin3.close();
702 }
703}
704
705//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
706
708{
709 if(Z < 1 || Z >= maxZ) { return -1;} //If Z is out of the supported return 0
710 //If necessary load data for Z
711 InitialiseForElement(theGamma, Z);
712 if(!fCrossSection[Z] || shell < 0 || shell >= fNShellsUsed[Z]) { return -1; }
713
714 if(Z>2)
715 return fShellCrossSection->GetComponentDataByIndex(Z, shell)->Energy(0);
716 else
717 return fCrossSection[Z]->Energy(0);
718}
719
720//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
721
723 const G4ParticleDefinition*, G4int Z)
724{
725 if (!fCrossSection[Z]) {
726#ifdef G4MULTITHREADED
727 G4MUTEXLOCK(&livPhotoeffMutex);
728 if (!fCrossSection[Z]) {
729#endif
730 ReadData(Z);
731#ifdef G4MULTITHREADED
732 }
733 G4MUTEXUNLOCK(&livPhotoeffMutex);
734#endif
735 }
736}
737
738//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4AtomicShellEnumerator
std::vector< G4Element * > G4ElementVector
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:85
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:251
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:254
std::mutex G4Mutex
Definition: G4Threading.hh:81
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:93
void InitialiseForComponent(G4int Z, G4int nComponents=0)
G4int GetComponentID(G4int Z, size_t idx)
G4PhysicsVector * GetComponentDataByIndex(G4int Z, size_t idx)
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
G4double GetValueForComponent(G4int Z, size_t idx, G4double kinEnergy)
G4int GetZasInt() const
Definition: G4Element.hh:131
static G4Gamma * Gamma()
Definition: G4Gamma.cc:85
void PutValues(std::size_t index, G4double e, G4double dataValue)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double energy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4ParticleChangeForGamma * fParticleChange
G4double GetBindingEnergy(G4int Z, G4int shell)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4LivermorePhotoElectricModel(const G4String &nam="LivermorePhElectric")
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double energy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
G4double GetDensity() const
Definition: G4Material.hh:178
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:231
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
G4SandiaTable * GetSandiaTable() const
Definition: G4Material.hh:227
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:651
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double Energy(std::size_t index) const
virtual void ScaleVector(G4double factorE, G4double factorV)
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
void SetSpline(G4bool)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
static G4ProductionCutsTable * GetProductionCutsTable()
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()
Definition: G4VEmModel.hh:611
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:133
G4bool IsMaster() const
Definition: G4VEmModel.hh:736
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:254
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:415
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:570
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:813
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:618
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
int G4lrint(double ad)
Definition: templates.hh:134
#define INT_MAX
Definition: templates.hh:90