| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Author: Luciano Pandola
28 //
29 // History:
30 // --------
31 // 08 Jan 2010 L Pandola First implementati
32 // 01 Feb 2011 L Pandola Suppress fake ener
33 // is active.
34 // Make sure that flu
35 // only if above thre
36 // 25 May 2011 L Pandola Renamed (make v200
37 // 10 Jun 2011 L Pandola Migrate atomic dee
38 // 07 Oct 2011 L Pandola Bug fix (potential
39 // 27 Sep 2013 L Pandola Migrate to MT para
40 // tables.
41 // 02 Oct 2013 L Pandola Rewrite sampling a
42 // to improve CPU per
43 //
44
45 #include "G4PenelopePhotoElectricModel.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4ParticleDefinition.hh"
49 #include "G4MaterialCutsCouple.hh"
50 #include "G4DynamicParticle.hh"
51 #include "G4PhysicsTable.hh"
52 #include "G4PhysicsFreeVector.hh"
53 #include "G4ElementTable.hh"
54 #include "G4Element.hh"
55 #include "G4AtomicTransitionManager.hh"
56 #include "G4AtomicShell.hh"
57 #include "G4Gamma.hh"
58 #include "G4Electron.hh"
59 #include "G4AutoLock.hh"
60 #include "G4LossTableManager.hh"
61 #include "G4Exp.hh"
62
63 //....oooOO0OOooo........oooOO0OOooo........oo
64
65 const G4int G4PenelopePhotoElectricModel::fMax
66 G4PhysicsTable* G4PenelopePhotoElectricModel::
67
68 //....oooOO0OOooo........oooOO0OOooo........oo
69
70 G4PenelopePhotoElectricModel::G4PenelopePhotoE
71 const G4String& nam)
72 :G4VEmModel(nam),fParticleChange(nullptr),fP
73 fAtomDeexcitation(nullptr),fIsInitialised(f
74 {
75 fIntrinsicLowEnergyLimit = 100.0*eV;
76 fIntrinsicHighEnergyLimit = 100.0*GeV;
77 // SetLowEnergyLimit(fIntrinsicLowEnergyLim
78 SetHighEnergyLimit(fIntrinsicHighEnergyLimit
79 //
80
81 if (part)
82 SetParticle(part);
83
84 fVerboseLevel= 0;
85 // Verbosity scale:
86 // 0 = nothing
87 // 1 = warning for energy non-conservation
88 // 2 = details of energy budget
89 // 3 = calculation of cross sections, file o
90 // 4 = entering in methods
91
92 //Mark this model as "applicable" for atomic
93 SetDeexcitationFlag(true);
94
95 fTransitionManager = G4AtomicTransitionManag
96 }
97
98 //....oooOO0OOooo........oooOO0OOooo........oo
99
100 G4PenelopePhotoElectricModel::~G4PenelopePhoto
101 {
102 if (IsMaster() || fLocalTable)
103 {
104 for(G4int i=0; i<=fMaxZ; ++i)
105 {
106 if(fLogAtomicShellXS[i]) {
107 fLogAtomicShellXS[i]->clearAndDestroy();
108 delete fLogAtomicShellXS[i];
109 fLogAtomicShellXS[i] = nullptr;
110 }
111 }
112 }
113 }
114
115 //....oooOO0OOooo........oooOO0OOooo........oo
116
117 void G4PenelopePhotoElectricModel::Initialise(
118 const G4DataVector& cuts)
119 {
120 if (fVerboseLevel > 3)
121 G4cout << "Calling G4PenelopePhotoElectri
122
123 fAtomDeexcitation = G4LossTableManager::Inst
124 //Issue warning if the AtomicDeexcitation ha
125 if (!fAtomDeexcitation)
126 {
127 G4cout << G4endl;
128 G4cout << "WARNING from G4PenelopePhotoE
129 G4cout << "Atomic de-excitation module i
130 G4cout << "any fluorescence/Auger emissi
131 G4cout << "Please make sure this is inte
132 }
133
134 SetParticle(particle);
135
136 //Only the master model creates/fills/destro
137 if (IsMaster() && particle == fParticle)
138 {
139 G4ProductionCutsTable* theCoupleTable =
140 G4ProductionCutsTable::GetProductionCutsTabl
141
142 for (G4int i=0;i<(G4int)theCoupleTable->
143 {
144 const G4Material* material =
145 theCoupleTable->GetMaterialCutsCouple(i)
146 const G4ElementVector* theElementVector =
147
148 for (std::size_t j=0;j<material->GetNumber
149 {
150 G4int iZ = theElementVector->at(j)->Ge
151 //read data files only in the master
152 if (!fLogAtomicShellXS[iZ])
153 ReadDataFile(iZ);
154 }
155 }
156
157 InitialiseElementSelectors(particle,cuts
158
159 if (fVerboseLevel > 0) {
160 G4cout << "Penelope Photo-Electric model v20
161 << "Energy range: "
162 << LowEnergyLimit() / MeV << " MeV -
163 << HighEnergyLimit() / GeV << " GeV";
164 }
165 }
166
167 if(fIsInitialised) return;
168 fParticleChange = GetParticleChangeForGamma(
169 fIsInitialised = true;
170
171 }
172
173 void G4PenelopePhotoElectricModel::InitialiseL
174 G4VEmModel *masterModel)
175 {
176 if (fVerboseLevel > 3)
177 G4cout << "Calling G4PenelopePhotoElectri
178 //
179 //Check that particle matches: one might hav
180 //for e+ and e-).
181 //
182 if (part == fParticle)
183 {
184 SetElementSelectors(masterModel->GetElem
185
186 //Get the const table pointers from the
187 const G4PenelopePhotoElectricModel* theM
188 static_cast<G4PenelopePhotoElectricModel*> (
189 for(G4int i=0; i<=fMaxZ; ++i)
190 fLogAtomicShellXS[i] = theModel->fLogAtomicS
191 //Same verbosity for all workers, as the
192 fVerboseLevel = theModel->fVerboseLevel;
193 }
194
195 return;
196 }
197
198 //....oooOO0OOooo........oooOO0OOooo........oo
199 namespace { G4Mutex PenelopePhotoElectricMode
200 G4double G4PenelopePhotoElectricModel::Compute
201 const G4ParticleDefinition*,
202 G4double energy,
203 G4double Z, G4double,
204 G4double, G4double)
205 {
206 //
207 // Penelope model v2008
208 //
209 if (fVerboseLevel > 3)
210 G4cout << "Calling ComputeCrossSectionPerA
211
212 G4int iZ = G4int(Z);
213
214 if (!fLogAtomicShellXS[iZ])
215 {
216 //If we are here, it means that Initiali
217 //not filled up. This can happen in a Un
218 if (fVerboseLevel > 0)
219 {
220 //Issue a G4Exception (warning) only in ve
221 G4ExceptionDescription ed;
222 ed << "Unable to retrieve the shell cross
223 ed << "This can happen only in Unit Tests
224 G4Exception("G4PenelopePhotoElectricModel:
225 "em2038",JustWarning,ed);
226 }
227 //protect file reading via autolock
228 G4AutoLock lock(&PenelopePhotoElectricMo
229 ReadDataFile(iZ);
230 lock.unlock();
231 }
232
233 G4double cross = 0;
234 G4PhysicsTable* theTable = fLogAtomicShellX
235 G4PhysicsFreeVector* totalXSLog = (G4Physics
236
237 if (!totalXSLog)
238 {
239 G4Exception("G4PenelopePhotoElectricMod
240 "em2039",FatalException,
241 "Unable to retrieve the total cross sec
242 return 0;
243 }
244 G4double logene = G4Log(energy);
245 G4double logXS = totalXSLog->Value(logene);
246 cross = G4Exp(logXS);
247
248 if (fVerboseLevel > 2)
249 G4cout << "Photoelectric cross section at
250 " = " << cross/barn << " barn" << G4endl
251 return cross;
252 }
253
254 //....oooOO0OOooo........oooOO0OOooo........oo
255
256 void G4PenelopePhotoElectricModel::SampleSecon
257 const G4MaterialCutsCouple* c
258 const G4DynamicParticle* aDyn
259 G4double,
260 G4double)
261 {
262 //
263 // Photoelectric effect, Penelope model v200
264 //
265 // The target atom and the target shell are
266 // database
267 // D.E. Cullen et al., Report UCRL-50400 (1
268 // The angular distribution of the electron
269 // according to the Sauter distribution from
270 // F. Sauter, Ann. Phys. 11 (1931) 454
271 // The energy of the final electron is given
272 // the binding energy. Fluorescence de-excit
273 // (to fill the vacancy) according to the ge
274 // J. Stepanek, Comp. Phys. Comm. 1206 pp 1
275
276 if (fVerboseLevel > 3)
277 G4cout << "Calling SamplingSecondaries() o
278
279 G4double photonEnergy = aDynamicGamma->GetKi
280
281 // always kill primary
282 fParticleChange->ProposeTrackStatus(fStopAnd
283 fParticleChange->SetProposedKineticEnergy(0.
284
285 if (photonEnergy <= fIntrinsicLowEnergyLimit
286 {
287 fParticleChange->ProposeLocalEnergyDepos
288 return ;
289 }
290
291 G4ParticleMomentum photonDirection = aDynami
292
293 // Select randomly one element in the curren
294 if (fVerboseLevel > 2)
295 G4cout << "Going to select element in " <<
296
297 // atom can be selected efficiently if eleme
298 const G4Element* anElement =
299 SelectRandomAtom(couple,G4Gamma::GammaDefi
300 G4int Z = anElement->GetZasInt();
301 if (fVerboseLevel > 2)
302 G4cout << "Selected " << anElement->GetNam
303
304 // Select the ionised shell in the current a
305 //shellIndex = 0 --> K shell
306 // 1-3 --> L shells
307 // 4-8 --> M shells
308 // 9 --> outer shells cumulative
309 //
310 std::size_t shellIndex = SelectRandomShell(Z
311
312 if (fVerboseLevel > 2)
313 G4cout << "Selected shell " << shellIndex
314
315 // Retrieve the corresponding identifier and
316 const G4AtomicTransitionManager* transitionM
317
318 //The number of shell cross section possibly
319 //might be different from the number of shel
320 //(namely, Penelope may contain more shell,
321 //In order to avoid a warning message from t
322 //add this protection. Results are anyway ch
323 //has a shellID>maxID, it sets the shellID t
324 std::size_t numberOfShells = (std::size_t) t
325 if (shellIndex >= numberOfShells)
326 shellIndex = numberOfShells-1;
327
328 const G4AtomicShell* shell = fTransitionMana
329 G4double bindingEnergy = shell->BindingEnerg
330
331 //Penelope considers only K, L and M shells.
332 //not included in the Penelope database. If
333 //shellIndex = 9, it means that an outer she
334 //Penelope recipe is to set bindingEnergy =
335 //to the electron) and to disregard fluoresc
336 if (shellIndex == 9)
337 bindingEnergy = 0.*eV;
338
339 G4double localEnergyDeposit = 0.0;
340 G4double cosTheta = 1.0;
341
342 // Primary outcoming electron
343 G4double eKineticEnergy = photonEnergy - bin
344
345 // There may be cases where the binding ener
346 // In such cases do not generate secondaries
347 if (eKineticEnergy > 0.)
348 {
349 // The electron is created
350 // Direction sampled from the Sauter dis
351 cosTheta = SampleElectronDirection(eKine
352 G4double sinTheta = std::sqrt(1-cosTheta
353 G4double phi = twopi * G4UniformRand() ;
354 G4double dirx = sinTheta * std::cos(phi)
355 G4double diry = sinTheta * std::sin(phi)
356 G4double dirz = cosTheta ;
357 G4ThreeVector electronDirection(dirx,dir
358 electronDirection.rotateUz(photonDirecti
359 G4DynamicParticle* electron = new G4Dyna
360 electronDirection,
361 eKineticEnergy);
362 fvect->push_back(electron);
363 }
364 else
365 bindingEnergy = photonEnergy;
366
367 G4double energyInFluorescence = 0; //testing
368 G4double energyInAuger = 0; //testing purpos
369
370 //Now, take care of fluorescence, if require
371 //recipe, I have to skip fluoresence complet
372 //(= sampling of a shell outer than K,L,M)
373 if (fAtomDeexcitation && shellIndex<9)
374 {
375 G4int index = couple->GetIndex();
376 if (fAtomDeexcitation->CheckDeexcitation
377 {
378 std::size_t nBefore = fvect->size();
379 fAtomDeexcitation->GenerateParticles(fvect
380 std::size_t nAfter = fvect->size();
381
382 if (nAfter > nBefore) //actual production
383 {
384 for (std::size_t j=nBefore;j<nAfter;++
385 {
386 G4double itsEnergy = ((*fvect)[j])->GetK
387 if (itsEnergy < bindingEnergy) // valid
388 {
389 bindingEnergy -= itsEnergy;
390 if (((*fvect)[j])->GetParticleDefini
391 energyInFluorescence += itsEnergy;
392 else if (((*fvect)[j])->GetParticleD
393 energyInAuger += itsEnergy;
394 }
395 else //invalid secondary: takes more tha
396 {
397 delete (*fvect)[j];
398 (*fvect)[j] = nullptr;
399 }
400 }
401 }
402 }
403 }
404
405 //Residual energy is deposited locally
406 localEnergyDeposit += bindingEnergy;
407
408 if (localEnergyDeposit < 0) //Should not be:
409 {
410 G4Exception("G4PenelopePhotoElectricMode
411 "em2099",JustWarning,"WARNING: Negative
412 localEnergyDeposit = 0;
413 }
414
415 fParticleChange->ProposeLocalEnergyDeposit(l
416
417 if (fVerboseLevel > 1)
418 {
419 G4cout << "-----------------------------
420 G4cout << "Energy balance from G4Penelop
421 G4cout << "Selected shell: " << WriteTar
422 anElement->GetName() << G4endl;
423 G4cout << "Incoming photon energy: " <<
424 G4cout << "-----------------------------
425 if (eKineticEnergy)
426 G4cout << "Outgoing electron " << eKineticEn
427 if (energyInFluorescence)
428 G4cout << "Fluorescence x-rays: " << energyI
429 if (energyInAuger)
430 G4cout << "Auger electrons: " << energyInAug
431 G4cout << "Local energy deposit " << loc
432 G4cout << "Total final state: " <<
433 (eKineticEnergy+energyInFluorescence+localEn
434 " keV" << G4endl;
435 G4cout << "-----------------------------
436 }
437 if (fVerboseLevel > 0)
438 {
439 G4double energyDiff =
440 std::fabs(eKineticEnergy+energyInFluorescenc
441 if (energyDiff > 0.05*keV)
442 {
443 G4cout << "Warning from G4PenelopePhotoEle
444 (eKineticEnergy+energyInFluorescence+loc
445 << " keV (final) vs. " <<
446 photonEnergy/keV << " keV (initial)" <<
447 G4cout << "-------------------------------
448 G4cout << "Energy balance from G4PenelopeP
449 G4cout << "Selected shell: " << WriteTarge
450 anElement->GetName() << G4endl;
451 G4cout << "Incoming photon energy: " << ph
452 G4cout << "-------------------------------
453 if (eKineticEnergy)
454 G4cout << "Outgoing electron " << eKinet
455 if (energyInFluorescence)
456 G4cout << "Fluorescence x-rays: " << ene
457 if (energyInAuger)
458 G4cout << "Auger electrons: " << energyI
459 G4cout << "Local energy deposit " << local
460 G4cout << "Total final state: " <<
461 (eKineticEnergy+energyInFluorescence+loc
462 " keV" << G4endl;
463 G4cout << "-------------------------------
464 }
465 }
466 }
467
468 //....oooOO0OOooo........oooOO0OOooo........oo
469
470 G4double G4PenelopePhotoElectricModel::SampleE
471 {
472 G4double costheta = 1.0;
473 if (energy>1*GeV) return costheta;
474
475 //1) initialize energy-dependent variables
476 // Variable naming according to Eq. (2.24) o
477 // (pag. 44)
478 G4double gamma = 1.0 + energy/electron_mass_
479 G4double gamma2 = gamma*gamma;
480 G4double beta = std::sqrt((gamma2-1.0)/gamma
481
482 // ac corresponds to "A" of Eq. (2.31)
483 //
484 G4double ac = (1.0/beta) - 1.0;
485 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(ga
486 G4double a2 = ac + 2.0;
487 G4double gtmax = 2.0*(a1 + 1.0/ac);
488
489 G4double tsam = 0;
490 G4double gtr = 0;
491
492 //2) sampling. Eq. (2.31) of Penelope Manual
493 // tsam = 1-std::cos(theta)
494 // gtr = rejection function according to Eq.
495 do{
496 G4double rand = G4UniformRand();
497 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(r
498 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
499 }while(G4UniformRand()*gtmax > gtr);
500 costheta = 1.0-tsam;
501
502 return costheta;
503 }
504
505 //....oooOO0OOooo........oooOO0OOooo........oo
506
507 void G4PenelopePhotoElectricModel::ReadDataFil
508 {
509 if (!IsMaster())
510 //Should not be here!
511 G4Exception("G4PenelopePhotoElectricModel:
512 "em0100",FatalException,"Worker thread in
513
514 if (fVerboseLevel > 2)
515 {
516 G4cout << "G4PenelopePhotoElectricModel:
517 G4cout << "Going to read PhotoElectric d
518 }
519
520 const char* path = G4FindDataDir("G4LEDATA
521 if(!path)
522 {
523 G4String excep = "G4PenelopePhotoElectri
524 G4Exception("G4PenelopePhotoElectricMode
525 "em0006",FatalException,excep);
526 return;
527 }
528
529 /*
530 Read the cross section file
531 */
532 std::ostringstream ost;
533 if (Z>9)
534 ost << path << "/penelope/photoelectric/pd
535 else
536 ost << path << "/penelope/photoelectric/pd
537 std::ifstream file(ost.str().c_str());
538 if (!file.is_open())
539 {
540 G4String excep = "G4PenelopePhotoElectri
541 G4Exception("G4PenelopePhotoElectricMode
542 "em0003",FatalException,excep);
543 }
544 //I have to know in advance how many points
545 //to initialize the G4PhysicsFreeVector()
546 std::size_t ndata=0;
547 G4String line;
548 while( getline(file, line) )
549 ndata++;
550 ndata -= 1;
551 //G4cout << "Found: " << ndata << " lines" <
552
553 file.clear();
554 file.close();
555 file.open(ost.str().c_str());
556
557 G4int readZ =0;
558 std::size_t nShells= 0;
559 file >> readZ >> nShells;
560
561 if (fVerboseLevel > 3)
562 G4cout << "Element Z=" << Z << " , nShells
563
564 //check the right file is opened.
565 if (readZ != Z || nShells <= 0 || nShells >
566 {
567 G4ExceptionDescription ed;
568 ed << "Corrupted data file for Z=" << Z
569 G4Exception("G4PenelopePhotoElectricMode
570 "em0005",FatalException,ed);
571 return;
572 }
573 G4PhysicsTable* thePhysicsTable = new G4Phys
574
575 //the table has to contain nShell+1 G4Physic
576 //(theTable)[0] --> total cross section
577 //(theTable)[ishell] --> cross section for s
578
579 //reserve space for the vectors
580 //everything is log-log
581 for (std::size_t i=0;i<nShells+1;++i)
582 thePhysicsTable->push_back(new G4PhysicsFr
583
584 std::size_t k =0;
585 for (k=0;k<ndata && !file.eof();++k)
586 {
587 G4double energy = 0;
588 G4double aValue = 0;
589 file >> energy ;
590 energy *= eV;
591 G4double logene = G4Log(energy);
592 //loop on the columns
593 for (std::size_t i=0;i<nShells+1;++i)
594 {
595 file >> aValue;
596 aValue *= barn;
597 G4PhysicsFreeVector* theVec = (G4PhysicsFr
598 if (aValue < 1e-40*cm2) //protection again
599 aValue = 1e-40*cm2;
600 theVec->PutValue(k,logene,G4Log(aValue));
601 }
602 }
603
604 if (fVerboseLevel > 2)
605 {
606 G4cout << "G4PenelopePhotoElectricModel:
607 << Z << G4endl;
608 }
609
610 fLogAtomicShellXS[Z] = thePhysicsTable;
611
612 file.close();
613 return;
614 }
615
616 //....oooOO0OOooo........oooOO0OOooo........oo
617
618 std::size_t G4PenelopePhotoElectricModel::GetN
619 {
620 if (!IsMaster())
621 //Should not be here!
622 G4Exception("G4PenelopePhotoElectricModel:
623 "em0100",FatalException,"Worker thread in
624
625 //read data files
626 if (!fLogAtomicShellXS[Z])
627 ReadDataFile(Z);
628 //now it should be ok
629 if (!fLogAtomicShellXS[Z])
630 {
631 G4ExceptionDescription ed;
632 ed << "Cannot find shell cross section
633 G4Exception("G4PenelopePhotoElectricMod
634 "em2038",FatalException,ed);
635 }
636 //one vector is allocated for the _total_ cr
637 std::size_t nEntries = fLogAtomicShellXS[Z]-
638 return (nEntries-1);
639 }
640
641 //....oooOO0OOooo........oooOO0OOooo........oo
642
643 G4double G4PenelopePhotoElectricModel::GetShel
644 {
645 //this forces also the loading of the data
646 std::size_t entries = GetNumberOfShellXS(Z);
647
648 if (shellID >= entries)
649 {
650 G4cout << "Element Z=" << Z << " has dat
651 G4cout << "so shellID should be from 0 t
652 return 0;
653 }
654
655 G4PhysicsTable* theTable = fLogAtomicShellX
656 //[0] is the total XS, shellID is in the ele
657 G4PhysicsFreeVector* totalXSLog = (G4Physics
658
659 if (!totalXSLog)
660 {
661 G4Exception("G4PenelopePhotoElectricMod
662 "em2039",FatalException,
663 "Unable to retrieve the total cross sec
664 return 0;
665 }
666 G4double logene = G4Log(energy);
667 G4double logXS = totalXSLog->Value(logene);
668 G4double cross = G4Exp(logXS);
669 if (cross < 2e-40*cm2) cross = 0;
670 return cross;
671 }
672
673 //....oooOO0OOooo........oooOO0OOooo........oo
674
675 G4String G4PenelopePhotoElectricModel::WriteTa
676 {
677 G4String theShell = "outer shell";
678 if (shellID == 0)
679 theShell = "K";
680 else if (shellID == 1)
681 theShell = "L1";
682 else if (shellID == 2)
683 theShell = "L2";
684 else if (shellID == 3)
685 theShell = "L3";
686 else if (shellID == 4)
687 theShell = "M1";
688 else if (shellID == 5)
689 theShell = "M2";
690 else if (shellID == 6)
691 theShell = "M3";
692 else if (shellID == 7)
693 theShell = "M4";
694 else if (shellID == 8)
695 theShell = "M5";
696
697 return theShell;
698 }
699
700 //....oooOO0OOooo........oooOO0OOooo........oo
701
702 void G4PenelopePhotoElectricModel::SetParticle
703 {
704 if(!fParticle) {
705 fParticle = p;
706 }
707 }
708
709 //....oooOO0OOooo........oooOO0OOooo........oo
710
711 std::size_t G4PenelopePhotoElectricModel::Sele
712 {
713 G4double logEnergy = G4Log(energy);
714
715 //Check if data have been read (it should be
716 if (!fLogAtomicShellXS[Z])
717 {
718 G4ExceptionDescription ed;
719 ed << "Cannot find shell cross section
720 G4Exception("G4PenelopePhotoElectricMod
721 "em2038",FatalException,ed);
722 }
723
724 G4PhysicsTable* theTable = fLogAtomicShellX
725
726 G4double sum = 0;
727 G4PhysicsFreeVector* totalXSLog = (G4Physics
728 G4double logXS = totalXSLog->Value(logEnergy
729 G4double totalXS = G4Exp(logXS);
730
731 //Notice: totalXS is the total cross section
732 //the sum of partialXS's, since these includ
733 //
734 // Therefore, here one have to consider the
735 // an outer shell. Conventionally, it is ind
736 //
737 G4double random = G4UniformRand()*totalXS;
738
739 for (std::size_t k=1;k<theTable->entries();+
740 {
741 //Add one shell
742 G4PhysicsFreeVector* partialXSLog = (G4P
743 G4double logXSLocal = partialXSLog->Valu
744 G4double partialXS = G4Exp(logXSLocal);
745 sum += partialXS;
746 if (random <= sum)
747 return k-1;
748 }
749 //none of the shells K, L, M: return outer s
750 return 9;
751 }
752