Geant4 9.6.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// $Id$
27//
28//
29// Author: Sebastien Incerti
30// 30 October 2008
31// on base of G4LowEnergyPhotoElectric developed by A.Forti and M.G.Pia
32//
33// 22 Oct 2012 A & V Ivanchenko Migration data structure to G4PhysicsVector
34//
35
37#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
50using namespace std;
51
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
55 const G4String& nam)
56 : G4VEmModel(nam),fParticleChange(0),maxZ(99),
57 nShellLimit(100),fDeexcitationActive(false),isInitialised(false),
58 fAtomDeexcitation(0)
59{
60 verboseLevel= 0;
61 // Verbosity scale:
62 // 0 = nothing
63 // 1 = warning for energy non-conservation
64 // 2 = details of energy budget
65 // 3 = calculation of cross sections, file openings, sampling of atoms
66 // 4 = entering in methods
67
68 theGamma = G4Gamma::Gamma();
69 theElectron = G4Electron::Electron();
70
71 // default generator
73
74 for(G4int i=0; i<maxZ; ++i) {
75 fCrossSection[i] = 0;
76 fCrossSectionLE[i] = 0;
77 fNShells[i] = 0;
78 fNShellsUsed[i] = 0;
79 }
80
81 if(verboseLevel>0) {
82 G4cout << "Livermore PhotoElectric is constructed "
83 << " nShellLimit= " << nShellLimit << G4endl;
84 }
85
86 //Mark this model as "applicable" for atomic deexcitation
88}
89
90//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91
93{
94 for(G4int i=0; i<maxZ; ++i) {
95 delete fCrossSection[i];
96 delete fCrossSectionLE[i];
97 }
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
102void
104 const G4DataVector&)
105{
106 if (verboseLevel > 2) {
107 G4cout << "Calling G4LivermorePhotoElectricModel::Initialise()" << G4endl;
108 }
109
110 char* path = getenv("G4LEDATA");
111
112 G4ProductionCutsTable* theCoupleTable =
114 G4int numOfCouples = theCoupleTable->GetTableSize();
115
116 for(G4int i=0; i<numOfCouples; ++i)
117 {
118 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
119 const G4Material* material = couple->GetMaterial();
120 const G4ElementVector* theElementVector = material->GetElementVector();
121 G4int nelm = material->GetNumberOfElements();
122
123 for (G4int j=0; j<nelm; ++j)
124 {
125 G4int Z = (G4int)(*theElementVector)[j]->GetZ();
126 if(Z < 1) { Z = 1; }
127 else if(Z > maxZ) { Z = maxZ; }
128 if(!fCrossSection[Z]) { ReadData(Z, path); }
129 }
130 }
131 //
132 if (verboseLevel > 2) {
133 G4cout << "Loaded cross section files for LivermorePhotoElectric model"
134 << G4endl;
135 }
136 if(!isInitialised) {
137 isInitialised = true;
139
140 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
141 }
142 fDeexcitationActive = false;
143 if(fAtomDeexcitation) {
144 fDeexcitationActive = fAtomDeexcitation->IsFluoActive();
145 }
146
147 if (verboseLevel > 0) {
148 G4cout << "LivermorePhotoElectric model is initialized " << G4endl
149 << G4endl;
150 }
151}
152
153//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154
157 G4double energy,
158 G4double ZZ, G4double,
160{
161 if (verboseLevel > 3) {
162 G4cout << "G4LivermorePhotoElectricModel::Calling ComputeCrossSectionPerAtom()"
163 << " Z= " << ZZ << " R(keV)= " << energy/keV << G4endl;
164 }
165 G4double cs = 0.0;
166 G4double gammaEnergy = energy;
167
168 G4int Z = G4lrint(ZZ);
169 if(Z < 1 || Z >= maxZ) { return cs; }
170
171 // element was not initialised
172 if(!fCrossSection[Z]) {
173 char* path = getenv("G4LEDATA");
174 ReadData(Z, path);
175 if(!fCrossSection[Z]) { return cs; }
176 }
177
178 G4int idx = fNShells[Z]*6 - 4;
179 if (gammaEnergy <= (fParam[Z])[idx-1]) { return cs; }
180
181 G4double x1 = 1.0/gammaEnergy;
182 G4double x2 = x1*x1;
183 G4double x3 = x2*x1;
184
185 // parameterisation
186 if(gammaEnergy >= (fParam[Z])[0]) {
187 G4double x4 = x2*x2;
188 cs = x1*((fParam[Z])[idx] + x1*(fParam[Z])[idx+1]
189 + x2*(fParam[Z])[idx+2] + x3*(fParam[Z])[idx+3]
190 + x4*(fParam[Z])[idx+4]);
191 // high energy part
192 } else if(gammaEnergy >= (fParam[Z])[1]) {
193 cs = x3*(fCrossSection[Z])->Value(gammaEnergy);
194
195 // low energy part
196 } else {
197 cs = x3*(fCrossSectionLE[Z])->Value(gammaEnergy);
198 }
199 if (verboseLevel > 1) {
200 G4cout << "LivermorePhotoElectricModel: E(keV)= " << gammaEnergy/keV
201 << " Z= " << Z << " cross(barn)= " << cs/barn << G4endl;
202 }
203 return cs;
204}
205
206//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
207
208void
210 std::vector<G4DynamicParticle*>* fvect,
211 const G4MaterialCutsCouple* couple,
212 const G4DynamicParticle* aDynamicGamma,
213 G4double,
214 G4double)
215{
216 G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
217 if (verboseLevel > 3) {
218 G4cout << "G4LivermorePhotoElectricModel::SampleSecondaries() Egamma(keV)= "
219 << gammaEnergy/keV << G4endl;
220 }
221
222 // kill incident photon
225
226 // Returns the normalized direction of the momentum
227 G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
228
229 // Select randomly one element in the current material
230 //G4cout << "Select random atom Egamma(keV)= " << gammaEnergy/keV << G4endl;
231 const G4Element* elm = SelectRandomAtom(couple->GetMaterial(),theGamma,
232 gammaEnergy);
233 G4int Z = G4lrint(elm->GetZ());
234
235 // Select the ionised shell in the current atom according to shell
236 // cross sections
237 // G4cout << "Select random shell Z= " << Z << G4endl;
238
239 if(Z >= maxZ) { Z = maxZ-1; }
240
241 // element was not initialised
242 if(!fCrossSection[Z]) {
243 char* path = getenv("G4LEDATA");
244 ReadData(Z, path);
245 if(!fCrossSection[Z]) {
247 return;
248 }
249 }
250
251 // shell index
252 size_t shellIdx = 0;
253 size_t nn = fNShellsUsed[Z];
254
255 if(nn > 1) {
256 if(gammaEnergy >= (fParam[Z])[0]) {
257 G4double x1 = 1.0/gammaEnergy;
258 G4double x2 = x1*x1;
259 G4double x3 = x2*x1;
260 G4double x4 = x3*x1;
261 G4int idx = nn*6 - 4;
262 // when do sampling common factors are not taken into account
263 // so cross section is not real
264 G4double cs0 = G4UniformRand()*((fParam[Z])[idx] + x1*(fParam[Z])[idx+1]
265 + x2*(fParam[Z])[idx+2]
266 + x3*(fParam[Z])[idx+3]
267 + x4*(fParam[Z])[idx+4]);
268 for(shellIdx=0; shellIdx<nn; ++shellIdx) {
269 idx = shellIdx*6 + 2;
270 if(gammaEnergy > (fParam[Z])[idx-1]) {
271 G4double cs = (fParam[Z])[idx] + x1*(fParam[Z])[idx+1]
272 + x2*(fParam[Z])[idx+2] + x3*(fParam[Z])[idx+3]
273 + x4*(fParam[Z])[idx+4];
274 if(cs >= cs0) { break; }
275 }
276 }
277 if(shellIdx >= nn) { shellIdx = nn-1; }
278
279 } else {
280
281 // when do sampling common factors are not taken into account
282 // so cross section is not real
283 G4double cs = G4UniformRand();
284
285 if(gammaEnergy >= (fParam[Z])[1]) {
286 cs *= (fCrossSection[Z])->Value(gammaEnergy);
287 } else {
288 cs *= (fCrossSectionLE[Z])->Value(gammaEnergy);
289 }
290
291 for(size_t j=0; j<nn; ++j) {
292 shellIdx = (size_t)fShellCrossSection.GetComponentID(Z, j);
293 if(gammaEnergy > (fParam[Z])[6*shellIdx+1]) {
294 cs -= fShellCrossSection.GetValueForComponent(Z, j, gammaEnergy);
295 }
296 if(cs <= 0.0 || j+1 == nn) { break; }
297 }
298 }
299 }
300
301 G4double bindingEnergy = (fParam[Z])[shellIdx*6 + 1];
302 //G4cout << "Z= " << Z << " shellIdx= " << shellIdx
303 // << " nShells= " << fNShells[Z]
304 // << " Ebind(keV)= " << bindingEnergy/keV
305 // << " Egamma(keV)= " << gammaEnergy/keV << G4endl;
306
307 const G4AtomicShell* shell = 0;
308
309 // no de-excitation from the last shell
310 if(fDeexcitationActive && shellIdx + 1 == nn) {
312 shell = fAtomDeexcitation->GetAtomicShell(Z, as);
313 }
314
315 // If binding energy of the selected shell is larger than photon energy
316 // do not generate secondaries
317 if(gammaEnergy < bindingEnergy) {
319 return;
320 }
321
322 // Primary outcoming electron
323 G4double eKineticEnergy = gammaEnergy - bindingEnergy;
324 G4double edep = bindingEnergy;
325
326 // Calculate direction of the photoelectron
327 G4ThreeVector electronDirection =
328 GetAngularDistribution()->SampleDirection(aDynamicGamma,
329 eKineticEnergy,
330 shellIdx,
331 couple->GetMaterial());
332
333 // The electron is created
334 G4DynamicParticle* electron = new G4DynamicParticle (theElectron,
335 electronDirection,
336 eKineticEnergy);
337 fvect->push_back(electron);
338
339 // Sample deexcitation
340 if(shell) {
341 G4int index = couple->GetIndex();
342 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
343 size_t nbefore = fvect->size();
344
345 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
346 size_t nafter = fvect->size();
347 if(nafter > nbefore) {
348 for (size_t j=nbefore; j<nafter; ++j) {
349 edep -= ((*fvect)[j])->GetKineticEnergy();
350 }
351 }
352 }
353 }
354 // energy balance - excitation energy left
355 if(edep > 0.0) {
357 }
358}
359
360//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
361
362void
363G4LivermorePhotoElectricModel::ReadData(G4int Z, const char* path)
364{
365 if (verboseLevel > 1)
366 {
367 G4cout << "Calling ReadData() of G4LivermoreGammaConversionModel"
368 << G4endl;
369 }
370
371 if(fCrossSection[Z]) { return; }
372
373 const char* datadir = path;
374
375 if(!datadir)
376 {
377 datadir = getenv("G4LEDATA");
378 if(!datadir)
379 {
380 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
381 "em0006",FatalException,
382 "Environment variable G4LEDATA not defined");
383 return;
384 }
385 }
386
387 // spline for photoeffect total x-section above K-shell
388 fCrossSection[Z] = new G4LPhysicsFreeVector();
389 fCrossSection[Z]->SetSpline(true);
390
391 std::ostringstream ost;
392 ost << datadir << "/livermore/phot/pe-cs-" << Z <<".dat";
393 std::ifstream fin(ost.str().c_str());
394 if( !fin.is_open()) {
396 ed << "G4LivermorePhotoElectricModel data file <" << ost.str().c_str()
397 << "> is not opened!" << G4endl;
398 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
399 "em0003",FatalException,
400 ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
401 return;
402 } else {
403 if(verboseLevel > 3) { G4cout << "File " << ost.str().c_str()
404 << " is opened by G4LivermorePhotoElectricModel" << G4endl;}
405 fCrossSection[Z]->Retrieve(fin, true);
406 fCrossSection[Z]->ScaleVector(MeV, barn);
407 fin.close();
408 }
409
410 // read fit parameters
411 G4int n1 = 0;
412 G4int n2 = 0;
413 G4double x;
414 std::ostringstream ost1;
415 ost1 << datadir << "/livermore/phot/pe-" << Z <<".dat";
416 std::ifstream fin1(ost1.str().c_str());
417 if( !fin1.is_open()) {
419 ed << "G4LivermorePhotoElectricModel data file <" << ost1.str().c_str()
420 << "> is not opened!" << G4endl;
421 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
422 "em0003",FatalException,
423 ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
424 return;
425 } else {
426 if(verboseLevel > 3) {
427 G4cout << "File " << ost1.str().c_str()
428 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
429 }
430 fin1 >> n1 >> n2 >> x;
431 fNShells[Z] = n1;
432 (fParam[Z]).reserve(6*n1+1);
433 (fParam[Z]).push_back(x*MeV);
434 for(G4int i=0; i<n1; ++i) {
435 for(G4int j=0; j<6; ++j) {
436 fin1 >> x;
437 if(0 == j) { x *= MeV; }
438 else { x *= barn; }
439 (fParam[Z]).push_back(x);
440 }
441 }
442 fin1.close();
443 }
444 // there is a possibility to used only main shells
445 if(nShellLimit < n2) { n2 = nShellLimit; }
446 fShellCrossSection.InitialiseForComponent(Z, n2);
447 fNShellsUsed[Z] = n2;
448
449 if(1 < n2) {
450 std::ostringstream ost2;
451 ost2 << datadir << "/livermore/phot/pe-ss-cs-" << Z <<".dat";
452 std::ifstream fin2(ost2.str().c_str());
453 if( !fin2.is_open()) {
455 ed << "G4LivermorePhotoElectricModel data file <" << ost2.str().c_str()
456 << "> is not opened!" << G4endl;
457 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
458 "em0003",FatalException,
459 ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
460 return;
461 } else {
462 if(verboseLevel > 3) {
463 G4cout << "File " << ost2.str().c_str()
464 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
465 }
466
467 G4int n3, n4;
468 G4double y;
469 for(G4int i=0; i<n2; ++i) {
470 fin2 >> x >> y >> n3 >> n4;
472 for(G4int j=0; j<n3; ++j) {
473 fin2 >> x >> y;
474 v->PutValues(j, x*MeV, y*barn);
475 }
476 fShellCrossSection.AddComponent(Z, n4, v);
477 }
478 fin2.close();
479 }
480 }
481
482 // no spline for photoeffect total x-section below K-shell
483 if(1 < fNShells[Z]) {
484 fCrossSectionLE[Z] = new G4LPhysicsFreeVector();
485
486 std::ostringstream ost3;
487 ost3 << datadir << "/livermore/phot/pe-le-cs-" << Z <<".dat";
488 std::ifstream fin3(ost3.str().c_str());
489 if( !fin3.is_open()) {
491 ed << "G4LivermorePhotoElectricModel data file <" << ost3.str().c_str()
492 << "> is not opened!" << G4endl;
493 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
494 "em0003",FatalException,
495 ed,"G4LEDATA version should be G4EMLOW6.32 or later.");
496 return;
497 } else {
498 if(verboseLevel > 3) {
499 G4cout << "File " << ost3.str().c_str()
500 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
501 }
502 fCrossSectionLE[Z]->Retrieve(fin3, true);
503 fCrossSectionLE[Z]->ScaleVector(MeV, barn);
504 fin3.close();
505 }
506 }
507}
508
509//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4AtomicShellEnumerator
std::vector< G4Element * > G4ElementVector
@ FatalException
@ fStopAndKill
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
void InitialiseForComponent(G4int Z, G4int nComponents=0)
G4int GetComponentID(G4int Z, size_t idx)
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
G4double GetValueForComponent(G4int Z, size_t idx, G4double kinEnergy)
G4double GetZ() const
Definition: G4Element.hh:131
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
void PutValues(size_t binNumber, G4double binValue, G4double dataValue)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4ParticleChangeForGamma * fParticleChange
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4LivermorePhotoElectricModel(const G4String &nam="LivermorePhElectric")
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:189
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
void SetProposedKineticEnergy(G4double proposedKinEnergy)
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
static G4ProductionCutsTable * GetProductionCutsTable()
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:508
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:109
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:286
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:459
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:641
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:515
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
int G4lrint(double ad)
Definition: templates.hh:163