source: trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08PhotoElectricModel.cc

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

geant4 tag 9.4
File size: 22.6 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4Penelope08PhotoElectricModel.cc,v 1.5 2010/07/28 07:09:16 pandola Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 08 Jan 2010 L Pandola First implementation
34
35#include "G4Penelope08PhotoElectricModel.hh"
36#include "G4ParticleDefinition.hh"
37#include "G4MaterialCutsCouple.hh"
38#include "G4ProductionCutsTable.hh"
39#include "G4DynamicParticle.hh"
40#include "G4PhysicsTable.hh"
41#include "G4PhysicsFreeVector.hh"
42#include "G4ElementTable.hh"
43#include "G4Element.hh"
44#include "G4AtomicTransitionManager.hh"
45#include "G4AtomicShell.hh"
46#include "G4Gamma.hh"
47#include "G4Electron.hh"
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50
51
52G4Penelope08PhotoElectricModel::G4Penelope08PhotoElectricModel(const G4ParticleDefinition*,
53 const G4String& nam)
54 :G4VEmModel(nam),isInitialised(false),logAtomicShellXS(0)
55{
56 fIntrinsicLowEnergyLimit = 100.0*eV;
57 fIntrinsicHighEnergyLimit = 100.0*GeV;
58 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
59 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
60 //
61 verboseLevel= 0;
62 // Verbosity scale:
63 // 0 = nothing
64 // 1 = warning for energy non-conservation
65 // 2 = details of energy budget
66 // 3 = calculation of cross sections, file openings, sampling of atoms
67 // 4 = entering in methods
68
69 //by default the model will inkove the atomic deexcitation
70 SetDeexcitationFlag(true);
71 ActivateAuger(false);
72}
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76G4Penelope08PhotoElectricModel::~G4Penelope08PhotoElectricModel()
77{
78 std::map <const G4int,G4PhysicsTable*>::iterator i;
79 if (logAtomicShellXS)
80 {
81 for (i=logAtomicShellXS->begin();i != logAtomicShellXS->end();i++)
82 {
83 G4PhysicsTable* tab = i->second;
84 tab->clearAndDestroy();
85 delete tab;
86 }
87 }
88 delete logAtomicShellXS;
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
93void G4Penelope08PhotoElectricModel::Initialise(const G4ParticleDefinition* particle,
94 const G4DataVector& cuts)
95{
96 if (verboseLevel > 3)
97 G4cout << "Calling G4Penelope08PhotoElectricModel::Initialise()" << G4endl;
98
99 // logAtomicShellXS is created only once, since it is never cleared
100 if (!logAtomicShellXS)
101 logAtomicShellXS = new std::map<const G4int,G4PhysicsTable*>;
102
103 InitialiseElementSelectors(particle,cuts);
104
105 if (verboseLevel > 0) {
106 G4cout << "Penelope Photo-Electric model is initialized " << G4endl
107 << "Energy range: "
108 << LowEnergyLimit() / MeV << " MeV - "
109 << HighEnergyLimit() / GeV << " GeV"
110 << G4endl;
111 }
112
113 if(isInitialised) return;
114 fParticleChange = GetParticleChangeForGamma();
115 isInitialised = true;
116}
117
118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119
120G4double G4Penelope08PhotoElectricModel::ComputeCrossSectionPerAtom(
121 const G4ParticleDefinition*,
122 G4double energy,
123 G4double Z, G4double,
124 G4double, G4double)
125{
126 //
127 // Penelope model.
128 //
129
130 if (verboseLevel > 3)
131 G4cout << "Calling ComputeCrossSectionPerAtom() of G4Penelope08PhotoElectricModel" << G4endl;
132
133 G4int iZ = (G4int) Z;
134
135 //read data files
136 if (!logAtomicShellXS->count(iZ))
137 ReadDataFile(iZ);
138 //now it should be ok
139 if (!logAtomicShellXS->count(iZ))
140 {
141 G4cout << "Problem in G4Penelope08PhotoElectricModel::ComputeCrossSectionPerAtom"
142 << G4endl;
143 G4Exception();
144 }
145
146 G4double cross = 0;
147
148 G4PhysicsTable* theTable = logAtomicShellXS->find(iZ)->second;
149 G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0];
150
151 if (!totalXSLog)
152 {
153 G4cout << "Problem in G4Penelope08PhotoElectricModel::ComputeCrossSectionPerAtom"
154 << G4endl;
155 G4Exception();
156 }
157 G4double logene = std::log(energy);
158 G4double logXS = totalXSLog->Value(logene);
159 cross = std::exp(logXS);
160
161 if (verboseLevel > 2)
162 G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z <<
163 " = " << cross/barn << " barn" << G4endl;
164 return cross;
165}
166
167//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168
169void G4Penelope08PhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
170 const G4MaterialCutsCouple* couple,
171 const G4DynamicParticle* aDynamicGamma,
172 G4double,
173 G4double)
174{
175 //
176 // Photoelectric effect, Penelope model
177 //
178 // The target atom and the target shell are sampled according to the Livermore
179 // database
180 // D.E. Cullen et al., Report UCRL-50400 (1989)
181 // The angular distribution of the electron in the final state is sampled
182 // according to the Sauter distribution from
183 // F. Sauter, Ann. Phys. 11 (1931) 454
184 // The energy of the final electron is given by the initial photon energy minus
185 // the binding energy. Fluorescence de-excitation is subsequently produced
186 // (to fill the vacancy) according to the general Geant4 G4DeexcitationManager:
187 // J. Stepanek, Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
188
189 if (verboseLevel > 3)
190 G4cout << "Calling SamplingSecondaries() of G4Penelope08PhotoElectricModel" << G4endl;
191
192 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
193
194 // always kill primary
195 fParticleChange->ProposeTrackStatus(fStopAndKill);
196 fParticleChange->SetProposedKineticEnergy(0.);
197
198 if (photonEnergy <= fIntrinsicLowEnergyLimit)
199 {
200 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
201 return ;
202 }
203
204 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
205
206 // Select randomly one element in the current material
207 if (verboseLevel > 2)
208 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
209
210 // atom can be selected efficiently if element selectors are initialised
211 const G4Element* anElement =
212 SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy);
213 G4int Z = (G4int) anElement->GetZ();
214 if (verboseLevel > 2)
215 G4cout << "Selected " << anElement->GetName() << G4endl;
216
217 // Select the ionised shell in the current atom according to shell cross sections
218 //shellIndex = 0 --> K shell
219 // 1-3 --> L shells
220 // 4-8 --> M shells
221 // 9 --> outer shells cumulatively
222 //
223 size_t shellIndex = SelectRandomShell(Z,photonEnergy);
224
225 if (verboseLevel > 2)
226 G4cout << "Selected shell " << shellIndex << " of element " << anElement->GetName() << G4endl;
227
228 // Retrieve the corresponding identifier and binding energy of the selected shell
229 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
230
231 //The number of shell cross section possibly reported in the Penelope database
232 //might be different from the number of shells in the G4AtomicTransitionManager
233 //(namely, Penelope may contain more shell, especially for very light elements).
234 //In order to avoid a warning message from the G4AtomicTransitionManager, I
235 //add this protection. Results are anyway changed, because when G4AtomicTransitionManager
236 //has a shellID>maxID, it sets the shellID to the last valid shell.
237 size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
238 if (shellIndex >= numberOfShells)
239 shellIndex = numberOfShells-1;
240
241 const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
242 G4double bindingEnergy = shell->BindingEnergy();
243 G4int shellId = shell->ShellId();
244
245 //Penelope considers only K, L and M shells. Cross sections of outer shells are
246 //not included in the Penelope database. If SelectRandomShell() returns
247 //shellIndex = 9, it means that an outer shell was ionized. In this case the
248 //Penelope recipe is to set bindingEnergy = 0 (the energy is entirely assigned
249 //to the electron) and to disregard fluorescence.
250 if (shellIndex == 9)
251 bindingEnergy = 0.*eV;
252
253
254 G4double localEnergyDeposit = 0.0;
255 G4double cosTheta = 1.0;
256
257 // Primary outcoming electron
258 G4double eKineticEnergy = photonEnergy - bindingEnergy;
259
260 // There may be cases where the binding energy of the selected shell is > photon energy
261 // In such cases do not generate secondaries
262 if (eKineticEnergy > 0.)
263 {
264 // The electron is created
265 // Direction sampled from the Sauter distribution
266 cosTheta = SampleElectronDirection(eKineticEnergy);
267 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
268 G4double phi = twopi * G4UniformRand() ;
269 G4double dirx = sinTheta * std::cos(phi);
270 G4double diry = sinTheta * std::sin(phi);
271 G4double dirz = cosTheta ;
272 G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
273 electronDirection.rotateUz(photonDirection);
274 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
275 electronDirection,
276 eKineticEnergy);
277 fvect->push_back(electron);
278 }
279 else
280 {
281 bindingEnergy = photonEnergy;
282 }
283
284 G4double energyInFluorescence = 0; //testing purposes
285
286 //Now, take care of fluorescence, if required
287 if(DeexcitationFlag() && Z > 5)
288 {
289 const G4ProductionCutsTable* theCoupleTable=
290 G4ProductionCutsTable::GetProductionCutsTable();
291 size_t indx = couple->GetIndex();
292 G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
293 G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
294
295 // Protection to avoid generating photons in the unphysical case of
296 // shell binding energy > photon energy
297 if (bindingEnergy > cutG || bindingEnergy > cutE)
298 {
299 deexcitationManager.SetCutForSecondaryPhotons(cutG);
300 deexcitationManager.SetCutForAugerElectrons(cutE);
301 std::vector<G4DynamicParticle*>* photonVector =
302 deexcitationManager.GenerateParticles(Z,shellId);
303 //Check for secondaries
304 if(photonVector)
305 {
306 for (size_t k=0; k< photonVector->size(); k++)
307 {
308 G4DynamicParticle* aPhoton = (*photonVector)[k];
309 if (aPhoton)
310 {
311 G4double itsEnergy = aPhoton->GetKineticEnergy();
312 if (itsEnergy <= bindingEnergy)
313 {
314 if(aPhoton->GetDefinition() == G4Gamma::Gamma())
315 energyInFluorescence += itsEnergy;
316 bindingEnergy -= itsEnergy;
317 fvect->push_back(aPhoton);
318 }
319 else
320 {
321 delete aPhoton;
322 (*photonVector)[k] = 0;
323 }
324 }
325 }
326 delete photonVector;
327 }
328 }
329 }
330 //Residual energy is deposited locally
331 localEnergyDeposit += bindingEnergy;
332
333 if (localEnergyDeposit < 0)
334 {
335 G4cout << "WARNING - "
336 << "G4Penelope08PhotoElectric::PostStepDoIt - Negative energy deposit"
337 << G4endl;
338 localEnergyDeposit = 0;
339 }
340
341 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
342
343 if (verboseLevel > 1)
344 {
345 G4cout << "-----------------------------------------------------------" << G4endl;
346 G4cout << "Energy balance from G4Penelope08PhotoElectric" << G4endl;
347 G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " <<
348 anElement->GetName() << G4endl;
349 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
350 G4cout << "-----------------------------------------------------------" << G4endl;
351 if (eKineticEnergy)
352 G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
353 G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
354 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
355 G4cout << "Total final state: " << (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV <<
356 " keV" << G4endl;
357 G4cout << "-----------------------------------------------------------" << G4endl;
358 }
359 if (verboseLevel > 0)
360 {
361 G4double energyDiff =
362 std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit-photonEnergy);
363 if (energyDiff > 0.05*keV)
364 G4cout << "Warning from G4Penelope08PhotoElectric: problem with energy conservation: " <<
365 (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV
366 << " keV (final) vs. " <<
367 photonEnergy/keV << " keV (initial)" << G4endl;
368 }
369}
370
371//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
372
373void G4Penelope08PhotoElectricModel::ActivateAuger(G4bool augerbool)
374{
375 if (!DeexcitationFlag() && augerbool)
376 {
377 G4cout << "WARNING - G4Penelope08PhotoElectricModel" << G4endl;
378 G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
379 G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
380 }
381 deexcitationManager.ActivateAugerElectronProduction(augerbool);
382 if (verboseLevel > 1)
383 G4cout << "Auger production set to " << augerbool << G4endl;
384}
385
386//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
387
388G4double G4Penelope08PhotoElectricModel::SampleElectronDirection(G4double energy)
389{
390 G4double costheta = 1.0;
391 if (energy>1*GeV) return costheta;
392
393 //1) initialize energy-dependent variables
394 // Variable naming according to Eq. (2.24) of Penelope Manual
395 // (pag. 44)
396 G4double gamma = 1.0 + energy/electron_mass_c2;
397 G4double gamma2 = gamma*gamma;
398 G4double beta = std::sqrt((gamma2-1.0)/gamma2);
399
400 // ac corresponds to "A" of Eq. (2.31)
401 //
402 G4double ac = (1.0/beta) - 1.0;
403 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
404 G4double a2 = ac + 2.0;
405 G4double gtmax = 2.0*(a1 + 1.0/ac);
406
407 G4double tsam = 0;
408 G4double gtr = 0;
409
410 //2) sampling. Eq. (2.31) of Penelope Manual
411 // tsam = 1-std::cos(theta)
412 // gtr = rejection function according to Eq. (2.28)
413 do{
414 G4double rand = G4UniformRand();
415 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
416 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
417 }while(G4UniformRand()*gtmax > gtr);
418 costheta = 1.0-tsam;
419
420
421 return costheta;
422}
423
424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
425
426void G4Penelope08PhotoElectricModel::ReadDataFile(G4int Z)
427{
428 if (verboseLevel > 2)
429 {
430 G4cout << "G4Penelope08PhotoElectricModel::ReadDataFile()" << G4endl;
431 G4cout << "Going to read PhotoElectric data files for Z=" << Z << G4endl;
432 }
433
434 char* path = getenv("G4LEDATA");
435 if (!path)
436 {
437 G4String excep = "G4Penelope08PhotoElectricModel - G4LEDATA environment variable not set!";
438 G4Exception(excep);
439 }
440
441 /*
442 Read the cross section file
443 */
444 std::ostringstream ost;
445 if (Z>9)
446 ost << path << "/penelope/photoelectric/pdgph" << Z << ".p08";
447 else
448 ost << path << "/penelope/photoelectric/pdgph0" << Z << ".p08";
449 std::ifstream file(ost.str().c_str());
450 if (!file.is_open())
451 {
452 G4String excep = "G4Penelope08PhotoElectricModel - data file " + G4String(ost.str()) + " not found!";
453 G4Exception(excep);
454 }
455 //I have to know in advance how many points are in the data list
456 //to initialize the G4PhysicsFreeVector()
457 size_t ndata=0;
458 G4String line;
459 while( getline(file, line) )
460 ndata++;
461 ndata -= 1;
462 //G4cout << "Found: " << ndata << " lines" << G4endl;
463
464 file.clear();
465 file.close();
466 file.open(ost.str().c_str());
467
468 G4int readZ =0;
469 size_t nShells= 0;
470 file >> readZ >> nShells;
471
472 if (verboseLevel > 3)
473 G4cout << "Element Z=" << Z << " , nShells = " << nShells << G4endl;
474
475 //check the right file is opened.
476 if (readZ != Z || nShells <= 0)
477 {
478 G4cout << "G4Penelope08PhotoElectricModel::ReadDataFile()" << G4endl;
479 G4cout << "Corrupted data file for Z=" << Z << G4endl;
480 G4Exception();
481 }
482 G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
483
484 //the table has to contain nShell+1 G4PhysicsFreeVectors,
485 //(theTable)[0] --> total cross section
486 //(theTable)[ishell] --> cross section for shell (ishell-1)
487
488 //reserve space for the vectors
489 //everything is log-log
490 for (size_t i=0;i<nShells+1;i++)
491 thePhysicsTable->push_back(new G4PhysicsFreeVector(ndata));
492
493 size_t k =0;
494 for (k=0;k<ndata && !file.eof();k++)
495 {
496 G4double energy = 0;
497 G4double aValue = 0;
498 file >> energy ;
499 energy *= eV;
500 G4double logene = std::log(energy);
501 //loop on the columns
502 for (size_t i=0;i<nShells+1;i++)
503 {
504 file >> aValue;
505 aValue *= barn;
506 G4PhysicsFreeVector* theVec = (G4PhysicsFreeVector*) ((*thePhysicsTable)[i]);
507 if (aValue < 1e-40*cm2) //protection against log(0)
508 aValue = 1e-40*cm2;
509 theVec->PutValue(k,logene,std::log(aValue));
510 }
511 }
512
513 if (verboseLevel > 2)
514 {
515 G4cout << "G4Penelope08PhotoElectricModel: read " << k << " points for element Z = "
516 << Z << G4endl;
517 }
518
519 logAtomicShellXS->insert(std::make_pair(Z,thePhysicsTable));
520
521 file.close();
522 return;
523}
524
525//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
526
527size_t G4Penelope08PhotoElectricModel::SelectRandomShell(G4int Z,G4double energy)
528{
529 G4double logEnergy = std::log(energy);
530
531 //Check if data have been read (it should be!)
532 if (!logAtomicShellXS->count(Z))
533 {
534 G4cout << "Problem in G4Penelope08PhotoElectricModel::SelectRandomShell" << G4endl;
535 G4cout << "Cannot find data for Z=" << Z << G4endl;
536 G4Exception();
537 }
538
539 size_t shellIndex = 0;
540
541 G4PhysicsTable* theTable = logAtomicShellXS->find(Z)->second;
542
543 G4DataVector* tempVector = new G4DataVector();
544
545 G4double sum = 0;
546 //loop on shell partial XS, retrieve the value for the given energy and store on
547 //a temporary vector
548 tempVector->push_back(sum); //first element is zero
549
550 G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0];
551 G4double logXS = totalXSLog->Value(logEnergy);
552 G4double totalXS = std::exp(logXS);
553
554 //Notice: totalXS is the total cross section and it does *not* correspond to
555 //the sum of partialXS's, since these include only K, L and M shells.
556 //
557 // Therefore, here one have to consider the possibility of ionisation of
558 // an outer shell. Conventionally, it is indicated with id=10 in Penelope
559 //
560
561 for (size_t k=1;k<theTable->entries();k++)
562 {
563 G4PhysicsFreeVector* partialXSLog = (G4PhysicsFreeVector*) (*theTable)[k];
564 G4double logXS = partialXSLog->Value(logEnergy);
565 G4double partialXS = std::exp(logXS);
566 sum += partialXS;
567 tempVector->push_back(sum);
568 }
569
570 tempVector->push_back(totalXS); //last element
571
572 G4double random = G4UniformRand()*totalXS;
573
574 /*
575 for (size_t i=0;i<tempVector->size(); i++)
576 G4cout << i << " " << (*tempVector)[i]/totalXS << G4endl;
577 */
578
579 //locate bin of tempVector
580 //Now one has to sample according to the elements in tempVector
581 //This gives the left edge of the interval...
582 size_t lowerBound = 0;
583 size_t upperBound = tempVector->size()-1;
584 while (lowerBound <= upperBound)
585 {
586 size_t midBin = (lowerBound + upperBound)/2;
587 if( random < (*tempVector)[midBin])
588 upperBound = midBin-1;
589 else
590 lowerBound = midBin+1;
591 }
592
593 shellIndex = upperBound;
594
595 delete tempVector;
596 return shellIndex;
597}
598
599//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
600
601size_t G4Penelope08PhotoElectricModel::GetNumberOfShellXS(G4int Z)
602{
603 //read data files
604 if (!logAtomicShellXS->count(Z))
605 ReadDataFile(Z);
606 //now it should be ok
607 if (!logAtomicShellXS->count(Z))
608 {
609 G4cout << "Problem in G4Penelope08PhotoElectricModel::GetNumberOfShellXS()"
610 << G4endl;
611 G4Exception();
612 }
613 //one vector is allocated for the _total_ cross section
614 size_t nEntries = logAtomicShellXS->find(Z)->second->entries();
615 return (nEntries-1);
616}
617
618//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
619
620G4double G4Penelope08PhotoElectricModel::GetShellCrossSection(G4int Z,size_t shellID,G4double energy)
621{
622 //this forces also the loading of the data
623 size_t entries = GetNumberOfShellXS(Z);
624
625 if (shellID >= entries)
626 {
627 G4cout << "Element Z=" << Z << " has data for " << entries << " shells only" << G4endl;
628 G4cout << "so shellID should be from 0 to " << entries-1 << G4endl;
629 return 0;
630 }
631
632 G4PhysicsTable* theTable = logAtomicShellXS->find(Z)->second;
633 //[0] is the total XS, shellID is in the element [shellID+1]
634 G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[shellID+1];
635
636 if (!totalXSLog)
637 {
638 G4cout << "Problem in G4Penelope08PhotoElectricModel::GetShellCrossSection()"
639 << G4endl;
640 G4Exception();
641 }
642 G4double logene = std::log(energy);
643 G4double logXS = totalXSLog->Value(logene);
644 G4double cross = std::exp(logXS);
645 if (cross < 2e-40*cm2) cross = 0;
646 return cross;
647}
648
649//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
650
651G4String G4Penelope08PhotoElectricModel::WriteTargetShell(size_t shellID)
652{
653 G4String theShell = "outer shell";
654 if (shellID == 0)
655 theShell = "K";
656 else if (shellID == 1)
657 theShell = "L1";
658 else if (shellID == 2)
659 theShell = "L2";
660 else if (shellID == 3)
661 theShell = "L3";
662 else if (shellID == 4)
663 theShell = "M1";
664 else if (shellID == 5)
665 theShell = "M2";
666 else if (shellID == 6)
667 theShell = "M3";
668 else if (shellID == 7)
669 theShell = "M4";
670 else if (shellID == 8)
671 theShell = "M5";
672
673 return theShell;
674}
Note: See TracBrowser for help on using the repository browser.