| 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 // 27 Jul 2010 L Pandola First complete i
32 // 18 Jan 2011 L.Pandola Stricter check o
33 // Should never hap
34 // 01 Feb 2011 L Pandola Suppress fake ener
35 // Make sure that flu
36 // above threshold
37 // 25 May 2011 L Pandola Renamed (make v200
38 // 26 Jan 2012 L Pandola Migration of Atomi
39 // 09 Mar 2012 L Pandola Moved the manageme
40 // cross sections to
41 // get normalized she
42 // 07 Oct 2013 L. Pandola Migration to MT
43 // 23 Jun 2015 L. Pandola Keep track of the
44 // atomic de-excitat
45 // 29 Aug 2018 L. Pandola Fix bug causing en
46 //
47
48 #include "G4PenelopeIonisationModel.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4ParticleDefinition.hh"
52 #include "G4MaterialCutsCouple.hh"
53 #include "G4ProductionCutsTable.hh"
54 #include "G4DynamicParticle.hh"
55 #include "G4AtomicTransitionManager.hh"
56 #include "G4AtomicShell.hh"
57 #include "G4Gamma.hh"
58 #include "G4Electron.hh"
59 #include "G4Positron.hh"
60 #include "G4PenelopeOscillatorManager.hh"
61 #include "G4PenelopeOscillator.hh"
62 #include "G4PenelopeCrossSection.hh"
63 #include "G4PhysicsFreeVector.hh"
64 #include "G4PhysicsLogVector.hh"
65 #include "G4LossTableManager.hh"
66 #include "G4PenelopeIonisationXSHandler.hh"
67 #include "G4EmParameters.hh"
68 #include "G4AutoLock.hh"
69
70 //....oooOO0OOooo........oooOO0OOooo........oo
71 namespace { G4Mutex PenelopeIonisationModelMu
72
73 G4PenelopeIonisationModel::G4PenelopeIonisatio
74 const G4String& nam)
75 :G4VEmModel(nam),fParticleChange(nullptr),fP
76 fCrossSectionHandler(nullptr),
77 fAtomDeexcitation(nullptr), fKineticEnergy1
78 fCosThetaPrimary(1.0),fEnergySecondary(0.*e
79 fCosThetaSecondary(0.0),fTargetOscillator(-
80 fIsInitialised(false),fPIXEflag(false),fLoc
81 {
82 fIntrinsicLowEnergyLimit = 100.0*eV;
83 fIntrinsicHighEnergyLimit = 100.0*GeV;
84 // SetLowEnergyLimit(fIntrinsicLowEnergyLim
85 SetHighEnergyLimit(fIntrinsicHighEnergyLimit
86 fNBins = 200;
87
88 if (part)
89 SetParticle(part);
90
91 //
92 fOscManager = G4PenelopeOscillatorManager::G
93 //
94 fVerboseLevel= 0;
95 // Verbosity scale:
96 // 0 = nothing
97 // 1 = warning for energy non-conservation
98 // 2 = details of energy budget
99 // 3 = calculation of cross sections, file o
100 // 4 = entering in methods
101
102 // Atomic deexcitation model activated by de
103 SetDeexcitationFlag(true);
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 G4PenelopeIonisationModel::~G4PenelopeIonisati
109 {
110 if (IsMaster() || fLocalTable)
111 {
112 if (fCrossSectionHandler)
113 delete fCrossSectionHandler;
114 }
115 }
116
117 //....oooOO0OOooo........oooOO0OOooo........oo
118
119 void G4PenelopeIonisationModel::Initialise(con
120 const G4DataVector& theCuts)
121 {
122 if (fVerboseLevel > 3)
123 G4cout << "Calling G4PenelopeIonisationMod
124
125 fAtomDeexcitation = G4LossTableManager::Inst
126 //Issue warning if the AtomicDeexcitation ha
127 if (!fAtomDeexcitation)
128 {
129 G4cout << G4endl;
130 G4cout << "WARNING from G4PenelopeIonisa
131 G4cout << "Atomic de-excitation module i
132 G4cout << "any fluorescence/Auger emissi
133 G4cout << "Please make sure this is inte
134 }
135
136 if (fAtomDeexcitation)
137 fPIXEflag = fAtomDeexcitation->IsPIXEActiv
138
139 //If the PIXE flag is active, the PIXE inter
140 //atomic de-excitation. The model does not n
141 //Issue warnings here
142 if (fPIXEflag && IsMaster() && particle==G4E
143 {
144 G4String theModel = G4EmParameters::Inst
145 G4cout << "=============================
146 G4cout << "The G4PenelopeIonisationModel
147 G4cout << "Atomic de-excitation will be
148 G4cout << "interface by using the shell
149 G4cout << "The built-in model procedure
150 G4cout << "*Please be sure this is inten
151 G4cout << "/process/em/pixe false" << G4
152 G4cout << "=============================
153 }
154
155 SetParticle(particle);
156
157 //Only the master model creates/manages the
158 //pointer to the table, and use it as readon
159 if (IsMaster() && particle == fParticle)
160 {
161 //Set the number of bins for the tables.
162 fNBins = (std::size_t) (20*std::log10(Hi
163 fNBins = std::max(fNBins,(std::size_t)10
164
165 //Clear and re-build the tables
166 if (fCrossSectionHandler)
167 {
168 delete fCrossSectionHandler;
169 fCrossSectionHandler = 0;
170 }
171 fCrossSectionHandler = new G4PenelopeIon
172 fCrossSectionHandler->SetVerboseLevel(fV
173
174 //Build tables for all materials
175 G4ProductionCutsTable* theCoupleTable =
176 G4ProductionCutsTable::GetProductionCutsTabl
177 for (G4int i=0;i<(G4int)theCoupleTable->
178 {
179 const G4Material* theMat =
180 theCoupleTable->GetMaterialCutsCouple(i)
181 //Forces the building of the cross section
182 fCrossSectionHandler->BuildXSTable(theMat,
183 IsMaster());
184 }
185
186 if (fVerboseLevel > 2) {
187 G4cout << "Penelope Ionisation model v2008 i
188 << "Energy range: "
189 << LowEnergyLimit() / keV << " keV -
190 << HighEnergyLimit() / GeV << " GeV.
191 << fNBins << " bins."
192 << G4endl;
193 }
194 }
195
196 if(fIsInitialised)
197 return;
198 fParticleChange = GetParticleChangeForLoss()
199 fIsInitialised = true;
200
201 return;
202 }
203
204 //....oooOO0OOooo........oooOO0OOooo........oo
205
206 void G4PenelopeIonisationModel::InitialiseLoca
207 G4VEmModel *masterModel)
208 {
209 if (fVerboseLevel > 3)
210 G4cout << "Calling G4PenelopeIonisationMo
211 //
212 //Check that particle matches: one might hav
213 //for e+ and e-).
214 //
215 if (part == fParticle)
216 {
217 //Get the const table pointers from the
218 const G4PenelopeIonisationModel* theMode
219 static_cast<G4PenelopeIonisationModel*> (mas
220
221 //Copy pointers to the data tables
222 fCrossSectionHandler = theModel->fCrossS
223
224 //copy data
225 fNBins = theModel->fNBins;
226
227 //Same verbosity for all workers, as the
228 fVerboseLevel = theModel->fVerboseLevel;
229 }
230 return;
231 }
232
233
234 //....oooOO0OOooo........oooOO0OOooo........oo
235
236 G4double G4PenelopeIonisationModel::CrossSecti
237 const G4ParticleDefinition*
238 theParticle,
239 G4double energy,
240 G4double cutEnergy,
241 G4double)
242 {
243 // Penelope model v2008 to calculate the cro
244 // threshold. It makes use of the Generalise
245 // D. Liljequist, J. Phys. D: Appl. Phys. 1
246 //
247 // The total cross section is calculated ana
248 // into account the atomic oscillators comin
249 //
250 // For incident e- the maximum energy allowe
251 // because particles are undistinghishable.
252 //
253 // The contribution is splitted in three par
254 // distant transverse collisions and close c
255 // its own analytical function.
256 // Fermi density correction is calculated an
257 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),
258 //
259 if (fVerboseLevel > 3)
260 G4cout << "Calling CrossSectionPerVolume()
261
262 SetupForMaterial(theParticle, material, ener
263
264 G4double totalCross = 0.0;
265 G4double crossPerMolecule = 0.;
266
267 //Either Initialize() was not called, or we
268 //not invoked
269 if (!fCrossSectionHandler)
270 {
271 //create a **thread-local** version of t
272 //Unit Tests
273 fLocalTable = true;
274 fCrossSectionHandler = new G4PenelopeIon
275 }
276
277 const G4PenelopeCrossSection* theXS =
278 fCrossSectionHandler->GetCrossSectionTable
279 material,
280 cutEnergy);
281 if (!theXS)
282 {
283 //If we are here, it means that Initiali
284 //not filled up. This can happen in a Un
285 if (fVerboseLevel > 0)
286 {
287 //Issue a G4Exception (warning) only in ve
288 G4ExceptionDescription ed;
289 ed << "Unable to retrieve the cross sectio
290 theParticle->GetParticleName() <<
291 " in " << material->GetName() << ", cut
292 ed << "This can happen only in Unit Tests
293 G4Exception("G4PenelopeIonisationModel::Cr
294 "em2038",JustWarning,ed);
295 }
296 //protect file reading via autolock
297 G4AutoLock lock(&PenelopeIonisationModel
298 fCrossSectionHandler->BuildXSTable(mater
299 lock.unlock();
300 //now it should be ok
301 theXS =
302 fCrossSectionHandler->GetCrossSectionTableFo
303 material,
304 cutEnergy);
305 }
306
307 if (theXS)
308 crossPerMolecule = theXS->GetHardCrossSect
309
310 G4double atomDensity = material->GetTotNbOfA
311 G4double atPerMol = fOscManager->GetAtomsPe
312
313 if (fVerboseLevel > 3)
314 G4cout << "Material " << material->GetName
315 "atoms per molecule" << G4endl;
316
317 G4double moleculeDensity = 0.;
318 if (atPerMol)
319 moleculeDensity = atomDensity/atPerMol;
320 G4double crossPerVolume = crossPerMolecule*m
321
322 if (fVerboseLevel > 2)
323 {
324 G4cout << "G4PenelopeIonisationModel " <<
325 G4cout << "Mean free path for delta emissi
326 energy/keV << " keV = " << (1./crossPerV
327 if (theXS)
328 totalCross = (theXS->GetTotalCrossSectio
329 G4cout << "Total free path for ionisation
330 energy/keV << " keV = " << (1./totalCros
331 }
332 return crossPerVolume;
333 }
334
335 //....oooOO0OOooo........oooOO0OOooo........oo
336
337 //This is a dummy method. Never inkoved by the
338 //a warning if one tries to get Cross Sections
339 //G4EmCalculator.
340 G4double G4PenelopeIonisationModel::ComputeCro
341 G4double,
342 G4double,
343 G4double,
344 G4double,
345 G4double)
346 {
347 G4cout << "*** G4PenelopeIonisationModel --
348 G4cout << "Penelope Ionisation model v2008 d
349 G4cout << "so the result is always zero. For
350 G4cout << "GetCrossSectionPerVolume() or Get
351 return 0;
352 }
353
354 //....oooOO0OOooo........oooOO0OOooo........oo
355
356 G4double G4PenelopeIonisationModel::ComputeDED
357 const G4ParticleDefinition* the
358 G4double kineticEnergy,
359 G4double cutEnergy)
360 {
361 // Penelope model v2008 to calculate the sto
362 // below the threshold. It makes use of the
363 // model from
364 // D. Liljequist, J. Phys. D: Appl. Phys. 1
365 //
366 // The stopping power is calculated analytic
367 // section from the GOS models, which includ
368 // distant longitudinal collisions, distant
369 // close collisions. Only the atomic oscilla
370 // (according to the threshold) are consider
371 // Differential cross sections have a differ
372 //
373 // Fermi density correction is calculated an
374 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),
375
376 if (fVerboseLevel > 3)
377 G4cout << "Calling ComputeDEDX() of G4Pene
378
379 //Either Initialize() was not called, or we
380 //not invoked
381 if (!fCrossSectionHandler)
382 {
383 //create a **thread-local** version of t
384 //Unit Tests
385 fLocalTable = true;
386 fCrossSectionHandler = new G4PenelopeIon
387 }
388
389 const G4PenelopeCrossSection* theXS =
390 fCrossSectionHandler->GetCrossSectionTable
391 cutEnergy);
392 if (!theXS)
393 {
394 //If we are here, it means that Initiali
395 //not filled up. This can happen in a Un
396 if (fVerboseLevel > 0)
397 {
398 //Issue a G4Exception (warning) only in ve
399 G4ExceptionDescription ed;
400 ed << "Unable to retrieve the cross sectio
401 theParticle->GetParticleName() <<
402 " in " << material->GetName() << ", cut
403 ed << "This can happen only in Unit Tests
404 G4Exception("G4PenelopeIonisationModel::Co
405 "em2038",JustWarning,ed);
406 }
407 //protect file reading via autolock
408 G4AutoLock lock(&PenelopeIonisationModel
409 fCrossSectionHandler->BuildXSTable(mater
410 lock.unlock();
411 //now it should be ok
412 theXS =
413 fCrossSectionHandler->GetCrossSectionTableFo
414 material,
415 cutEnergy);
416 }
417
418 G4double sPowerPerMolecule = 0.0;
419 if (theXS)
420 sPowerPerMolecule = theXS->GetSoftStopping
421
422 G4double atomDensity = material->GetTotNbOfA
423 G4double atPerMol = fOscManager->GetAtomsPe
424
425 G4double moleculeDensity = 0.;
426 if (atPerMol)
427 moleculeDensity = atomDensity/atPerMol;
428 G4double sPowerPerVolume = sPowerPerMolecule
429
430 if (fVerboseLevel > 2)
431 {
432 G4cout << "G4PenelopeIonisationModel " <
433 G4cout << "Stopping power < " << cutEner
434 kineticEnergy/keV << " keV = " <<
435 sPowerPerVolume/(keV/mm) << " keV/mm" << G4e
436 }
437 return sPowerPerVolume;
438 }
439
440 //....oooOO0OOooo........oooOO0OOooo........oo
441
442 G4double G4PenelopeIonisationModel::MinEnergyC
443 const G4MaterialCutsCouple*)
444 {
445 return fIntrinsicLowEnergyLimit;
446 }
447
448 //....oooOO0OOooo........oooOO0OOooo........oo
449
450 void G4PenelopeIonisationModel::SampleSecondar
451 const G4MaterialCutsCouple* coup
452 const G4DynamicParticle* aDynami
453 G4double cutE, G4double)
454 {
455 // Penelope model v2008 to sample the final
456 // It makes use of the Generalised Oscillato
457 // D. Liljequist, J. Phys. D: Appl. Phys. 1
458 //
459 // The GOS model is used to calculate the in
460 // the atomic oscillators coming in the play
461 // contributions (distant longitudinal colli
462 // close collisions). Then the target shell
463 // sampled. Final state of the delta-ray (en
464 // to the analytical distributions (dSigma/d
465 // channels.
466 // For e-, the maximum energy for the delta-
467 // particles are indistinghusbable), while i
468 // e+.
469 // The efficiency of the random sampling alg
470 // decreases when initial and cutoff energy
471 // and 1 keV threshold, 99% for 10-MeV prima
472 // Differential cross sections have a differ
473 //
474 // WARNING: The model provides an _average_
475 // Anyway, the energy spectrum associated to
476 // atomic shell is approximated as a single
477 // show _unphysical_ narrow peaks at energie
478 // resonance energies. The spurious speaks a
479 // multiple inelastic collisions.
480 //
481 // The model determines also the original sh
482 // in order to produce fluorescence de-excit
483 //
484 // Fermi density correction is calculated an
485 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),
486
487 if (fVerboseLevel > 3)
488 G4cout << "Calling SamplingSecondaries() o
489
490 G4double kineticEnergy0 = aDynamicParticle->
491 const G4ParticleDefinition* theParticle = aD
492
493 if (kineticEnergy0 <= fIntrinsicLowEnergyLim
494 {
495 fParticleChange->SetProposedKineticEnergy(
496 fParticleChange->ProposeLocalEnergyDeposit
497 return ;
498 }
499
500 const G4Material* material = couple->GetMate
501 const G4PenelopeOscillatorTable* theTable =
502
503 G4ParticleMomentum particleDirection0 = aDyn
504
505 //Initialise final-state variables. The prop
506 // SampleFinalStateElectron() and SampleFina
507 fKineticEnergy1=kineticEnergy0;
508 fCosThetaPrimary=1.0;
509 fEnergySecondary=0.0;
510 fCosThetaSecondary=1.0;
511 fTargetOscillator = -1;
512
513 if (theParticle == G4Electron::Electron())
514 SampleFinalStateElectron(material,cutE,kin
515 else if (theParticle == G4Positron::Positron
516 SampleFinalStatePositron(material,cutE,kin
517 else
518 {
519 G4ExceptionDescription ed;
520 ed << "Invalid particle " << theParticle
521 G4Exception("G4PenelopeIonisationModel::
522 "em0001",FatalException,ed);
523
524 }
525 if (fEnergySecondary == 0) return;
526
527 if (fVerboseLevel > 3)
528 {
529 G4cout << "G4PenelopeIonisationModel::Sam
530 theParticle->GetParticleName() << G4endl;
531 G4cout << "Final eKin = " << fKineticEne
532 G4cout << "Final cosTheta = " << fCosThe
533 G4cout << "Delta-ray eKin = " << fEnergy
534 G4cout << "Delta-ray cosTheta = " << fCo
535 G4cout << "Oscillator: " << fTargetOscil
536 }
537
538 //Update the primary particle
539 G4double sint = std::sqrt(1. - fCosThetaPrim
540 G4double phiPrimary = twopi * G4UniformRand
541 G4double dirx = sint * std::cos(phiPrimary);
542 G4double diry = sint * std::sin(phiPrimary);
543 G4double dirz = fCosThetaPrimary;
544
545 G4ThreeVector electronDirection1(dirx,diry,d
546 electronDirection1.rotateUz(particleDirectio
547
548 if (fKineticEnergy1 > 0)
549 {
550 fParticleChange->ProposeMomentumDirectio
551 fParticleChange->SetProposedKineticEnerg
552 }
553 else
554 fParticleChange->SetProposedKineticEnergy(
555
556 //Generate the delta ray
557 G4double ionEnergyInPenelopeDatabase =
558 (*theTable)[fTargetOscillator]->GetIonisat
559
560 //Now, try to handle fluorescence
561 //Notice: merged levels are indicated with Z
562 G4int shFlag = (*theTable)[fTargetOscillator
563 G4int Z = (G4int) (*theTable)[fTargetOscilla
564
565 //initialize here, then check photons create
566 const G4AtomicTransitionManager* transitionM
567 G4double bindingEnergy = 0.*eV;
568
569 const G4AtomicShell* shell = nullptr;
570 //Real level
571 if (Z > 0 && shFlag<30)
572 {
573 shell = transitionManager->Shell(Z,shFla
574 bindingEnergy = shell->BindingEnergy();
575 //shellId = shell->ShellId();
576 }
577
578 //correct the fEnergySecondary to account fo
579 //database of ionisation energies is in gene
580 //from the fluorescence database used in Gea
581 fEnergySecondary += ionEnergyInPenelopeDatab
582
583 G4double localEnergyDeposit = bindingEnergy;
584 //testing purposes only
585 G4double energyInFluorescence = 0;
586 G4double energyInAuger = 0;
587
588 if (fEnergySecondary < 0)
589 {
590 //It means that there was some problem/m
591 //In this case, the available energy is
592 //to the Penelope database, but not acco
593 //Full residual energy is deposited loca
594 localEnergyDeposit += fEnergySecondary;
595 fEnergySecondary = 0.0;
596 }
597
598 //Notice: shell might be nullptr (invalid!)
599 //Disable the built-in de-excitation of the
600 //case, the PIXE interface takes care (stati
601 //de-excitation.
602 //Now, take care of fluorescence, if require
603 if (fAtomDeexcitation && !fPIXEflag && shell
604 {
605 G4int index = couple->GetIndex();
606 if (fAtomDeexcitation->CheckDeexcitation
607 {
608 std::size_t nBefore = fvect->size();
609 fAtomDeexcitation->GenerateParticles(fvect
610 std::size_t nAfter = fvect->size();
611
612 if (nAfter>nBefore) //actual production of
613 {
614 for (std::size_t j=nBefore;j<nAfter;++
615 {
616 G4double itsEnergy = ((*fvect)[j])->GetK
617 if (itsEnergy < localEnergyD
618 {
619 localEnergyDeposit -= itsEnergy;
620 if (((*fvect)[j])->GetParticleDefini
621 energyInFluorescence += itsEnergy;
622 else if (((*fvect)[j])->GetParticleD
623 energyInAuger += itsEnergy;
624 }
625 else //invalid secondary: ta
626 {
627 delete (*fvect)[j];
628 (*fvect)[j] = nullptr;
629 }
630 }
631 }
632 }
633 }
634
635 // Generate the delta ray --> to be done onl
636 if (fEnergySecondary > cutE)
637 {
638 G4DynamicParticle* electron = nullptr;
639 G4double sinThetaE = std::sqrt(1.-fCosTh
640 G4double phiEl = phiPrimary+pi; //pi wit
641 G4double xEl = sinThetaE * std::cos(phiE
642 G4double yEl = sinThetaE * std::sin(phiE
643 G4double zEl = fCosThetaSecondary;
644 G4ThreeVector eDirection(xEl,yEl,zEl); /
645 eDirection.rotateUz(particleDirection0);
646 electron = new G4DynamicParticle (G4Elec
647 eDirection,fEnergySecondary) ;
648 fvect->push_back(electron);
649 }
650 else
651 {
652 localEnergyDeposit += fEnergySecondary;
653 fEnergySecondary = 0;
654 }
655
656 if (localEnergyDeposit < 0) //Should not be:
657 {
658 G4Exception("G4PenelopeIonisationModel::
659 "em2099",JustWarning,"WARNING: Negative
660 localEnergyDeposit=0.;
661 }
662 fParticleChange->ProposeLocalEnergyDeposit(l
663
664 if (fVerboseLevel > 1)
665 {
666 G4cout << "-----------------------------
667 G4cout << "Energy balance from G4Penelop
668 G4cout << "Incoming primary energy: " <<
669 G4cout << "-----------------------------
670 G4cout << "Outgoing primary energy: " <<
671 G4cout << "Delta ray " << fEnergySeconda
672 if (energyInFluorescence)
673 G4cout << "Fluorescence x-rays: " << energyI
674 if (energyInAuger)
675 G4cout << "Auger electrons: " << energyInAug
676 G4cout << "Local energy deposit " << loc
677 G4cout << "Total final state: " << (fEne
678 localEnergyDeposit+energyInAuger)/
679 " keV" << G4endl;
680 G4cout << "-----------------------------
681 }
682
683 if (fVerboseLevel > 0)
684 {
685 G4double energyDiff = std::fabs(fEnergyS
686 localEnergyDeposit+energyInAuger
687 if (energyDiff > 0.05*keV)
688 G4cout << "Warning from G4PenelopeIonisation
689 (fEnergySecondary+energyInFluorescence+fKi
690 " keV (final) vs. " <<
691 kineticEnergy0/keV << " keV (initial)" <<
692 }
693
694 }
695
696 //....oooOO0OOooo........oooOO0OOooo........oo
697 void G4PenelopeIonisationModel::SampleFinalSta
698 G4double cutEnergy,
699 G4double kineticEnergy)
700 {
701 // This method sets the final ionisation par
702 // fKineticEnergy1, fCosThetaPrimary (= upda
703 // fEnergySecondary, fCosThetaSecondary (= i
704 // fTargetOscillator (= ionised oscillator)
705 //
706 // The method implements SUBROUTINE EINa of
707 //
708
709 G4PenelopeOscillatorTable* theTable = fOscMa
710 std::size_t numberOfOscillators = theTable->
711 const G4PenelopeCrossSection* theXS =
712 fCrossSectionHandler->GetCrossSectionTable
713 cutEnergy);
714 G4double delta = fCrossSectionHandler->GetDe
715
716 // Selection of the active oscillator
717 G4double TST = G4UniformRand();
718 fTargetOscillator = G4int(numberOfOscillator
719 G4double XSsum = 0.;
720
721 for (std::size_t i=0;i<numberOfOscillators-1
722 {
723 XSsum += theXS->GetNormalizedShellCrossS
724
725 if (XSsum > TST)
726 {
727 fTargetOscillator = (G4int) i;
728 break;
729 }
730 }
731
732 if (fVerboseLevel > 3)
733 {
734 G4cout << "SampleFinalStateElectron: sam
735 fTargetOscillator << "." << G4endl;
736 G4cout << "Ionisation energy: " <<
737 (*theTable)[fTargetOscillator]->GetIonisatio
738 " eV " << G4endl;
739 G4cout << "Resonance energy: : " <<
740 (*theTable)[fTargetOscillator]->GetResonance
741 << G4endl;
742 }
743 //Constants
744 G4double rb = kineticEnergy + 2.0*electron_m
745 G4double gam = 1.0+kineticEnergy/electron_ma
746 G4double gam2 = gam*gam;
747 G4double beta2 = (gam2-1.0)/gam2;
748 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/g
749
750 //Partial cross section of the active oscill
751 G4double resEne = (*theTable)[fTargetOscilla
752 G4double invResEne = 1.0/resEne;
753 G4double ionEne = (*theTable)[fTargetOscilla
754 G4double cutoffEne = (*theTable)[fTargetOsci
755 G4double XHDL = 0.;
756 G4double XHDT = 0.;
757 G4double QM = 0.;
758 G4double cps = 0.;
759 G4double cp = 0.;
760
761 //Distant excitations
762 if (resEne > cutEnergy && resEne < kineticEn
763 {
764 cps = kineticEnergy*rb;
765 cp = std::sqrt(cps);
766 G4double XHDT0 = std::max(G4Log(gam2)-be
767 if (resEne > 1.0e-6*kineticEnergy)
768 {
769 G4double cpp = std::sqrt((kineticEnergy-re
770 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_
771 }
772 else
773 {
774 QM = resEne*resEne/(beta2*2.0*electron_mas
775 QM *= (1.0-QM*0.5/electron_mass_c2);
776 }
777 if (QM < cutoffEne)
778 {
779 XHDL = G4Log(cutoffEne*(QM+2.0*electron_ma
780 *invResEne;
781 XHDT = XHDT0*invResEne;
782 }
783 else
784 {
785 QM = cutoffEne;
786 XHDL = 0.;
787 XHDT = 0.;
788 }
789 }
790 else
791 {
792 QM = cutoffEne;
793 cps = 0.;
794 cp = 0.;
795 XHDL = 0.;
796 XHDT = 0.;
797 }
798
799 //Close collisions
800 G4double EE = kineticEnergy + ionEne;
801 G4double wmaxc = 0.5*EE;
802 G4double wcl = std::max(cutEnergy,cutoffEne)
803 G4double rcl = wcl/EE;
804 G4double XHC = 0.;
805 if (wcl < wmaxc)
806 {
807 G4double rl1 = 1.0-rcl;
808 G4double rrl1 = 1.0/rl1;
809 XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
810 (1.0-amol)*G4Log(rcl*rrl1))/EE;
811 }
812
813 //Total cross section per molecule for the a
814 G4double XHTOT = XHC + XHDL + XHDT;
815
816 //very small cross section, do nothing
817 if (XHTOT < 1.e-14*barn)
818 {
819 fKineticEnergy1=kineticEnergy;
820 fCosThetaPrimary=1.0;
821 fEnergySecondary=0.0;
822 fCosThetaSecondary=1.0;
823 fTargetOscillator = G4int(numberOfOscill
824 return;
825 }
826
827 //decide which kind of interaction we'll hav
828 TST = XHTOT*G4UniformRand();
829
830 // Hard close collision
831 G4double TS1 = XHC;
832
833 if (TST < TS1)
834 {
835 G4double A = 5.0*amol;
836 G4double ARCL = A*0.5*rcl;
837 G4double rk=0.;
838 G4bool loopAgain = false;
839 do
840 {
841 loopAgain = false;
842 G4double fb = (1.0+ARCL)*G4UniformRand();
843 if (fb < 1)
844 rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));
845 else
846 rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
847 G4double rk2 = rk*rk;
848 G4double rkf = rk/(1.0-rk);
849 G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+r
850 if (G4UniformRand()*(1.0+A*rk2) > phi)
851 loopAgain = true;
852 }while(loopAgain);
853 //energy and scattering angle (primary e
854 G4double deltaE = rk*EE;
855
856 fKineticEnergy1 = kineticEnergy - deltaE
857 fCosThetaPrimary = std::sqrt(fKineticEne
858 //energy and scattering angle of the del
859 fEnergySecondary = deltaE - ionEne; //su
860 fCosThetaSecondary= std::sqrt(deltaE*rb/
861 if (fVerboseLevel > 3)
862 G4cout << "SampleFinalStateElectron: sampled
863 return;
864 }
865
866 //Hard distant longitudinal collisions
867 TS1 += XHDL;
868 G4double deltaE = resEne;
869 fKineticEnergy1 = kineticEnergy - deltaE;
870
871 if (TST < TS1)
872 {
873 G4double QS = QM/(1.0+QM*0.5/electron_ma
874 G4double Q = QS/(std::pow((QS/cutoffEne)
875 - (QS*0.5/electron_mass_c2));
876 G4double QTREV = Q*(Q+2.0*electron_mass_
877 G4double cpps = fKineticEnergy1*(fKineti
878 fCosThetaPrimary = (cpps+cps-QTREV)/(2.0
879 if (fCosThetaPrimary > 1.)
880 fCosThetaPrimary = 1.0;
881 //energy and emission angle of the delta
882 fEnergySecondary = deltaE - ionEne;
883 fCosThetaSecondary = 0.5*(deltaE*(kineti
884 if (fCosThetaSecondary > 1.0)
885 fCosThetaSecondary = 1.0;
886 if (fVerboseLevel > 3)
887 G4cout << "SampleFinalStateElectron: sampled
888 return;
889 }
890
891 //Hard distant transverse collisions
892 fCosThetaPrimary = 1.0;
893 //energy and emission angle of the delta ray
894 fEnergySecondary = deltaE - ionEne;
895 fCosThetaSecondary = 0.5;
896 if (fVerboseLevel > 3)
897 G4cout << "SampleFinalStateElectron: sampl
898
899 return;
900 }
901
902 //....oooOO0OOooo........oooOO0OOooo........oo
903
904 void G4PenelopeIonisationModel::SampleFinalSta
905 G4double cutEnergy,
906 G4double kineticEnergy)
907 {
908 // This method sets the final ionisation par
909 // fKineticEnergy1, fCosThetaPrimary (= upda
910 // fEnergySecondary, fCosThetaSecondary (= i
911 // fTargetOscillator (= ionised oscillator)
912 //
913 // The method implements SUBROUTINE PINa of
914 //
915
916 G4PenelopeOscillatorTable* theTable = fOscMa
917 std::size_t numberOfOscillators = theTable->
918 const G4PenelopeCrossSection* theXS =
919 fCrossSectionHandler->GetCrossSectionTable
920 cutEnergy);
921 G4double delta = fCrossSectionHandler->GetDe
922
923 // Selection of the active oscillator
924 G4double TST = G4UniformRand();
925 fTargetOscillator = G4int(numberOfOscillator
926 G4double XSsum = 0.;
927 for (std::size_t i=0;i<numberOfOscillators-1
928 {
929 XSsum += theXS->GetNormalizedShellCrossS
930 if (XSsum > TST)
931 {
932 fTargetOscillator = (G4int) i;
933 break;
934 }
935 }
936
937 if (fVerboseLevel > 3)
938 {
939 G4cout << "SampleFinalStatePositron: sam
940 fTargetOscillator << "." << G4endl;
941 G4cout << "Ionisation energy: " << (*the
942 << " eV " << G4endl;
943 G4cout << "Resonance energy: : " << (*th
944 << " eV " << G4endl;
945 }
946
947 //Constants
948 G4double rb = kineticEnergy + 2.0*electron_m
949 G4double gam = 1.0+kineticEnergy/electron_ma
950 G4double gam2 = gam*gam;
951 G4double beta2 = (gam2-1.0)/gam2;
952 G4double g12 = (gam+1.0)*(gam+1.0);
953 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/g
954 //Bhabha coefficients
955 G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0
956 G4double bha2 = amol*(3.0+1.0/g12);
957 G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
958 G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12
959
960 //
961 //Partial cross section of the active oscill
962 //
963 G4double resEne = (*theTable)[fTargetOscilla
964 G4double invResEne = 1.0/resEne;
965 G4double ionEne = (*theTable)[fTargetOscilla
966 G4double cutoffEne = (*theTable)[fTargetOsci
967
968 G4double XHDL = 0.;
969 G4double XHDT = 0.;
970 G4double QM = 0.;
971 G4double cps = 0.;
972 G4double cp = 0.;
973
974 //Distant excitations XS (same as for electr
975 if (resEne > cutEnergy && resEne < kineticEn
976 {
977 cps = kineticEnergy*rb;
978 cp = std::sqrt(cps);
979 G4double XHDT0 = std::max(G4Log(gam2)-be
980 if (resEne > 1.0e-6*kineticEnergy)
981 {
982 G4double cpp = std::sqrt((kineticEnergy-re
983 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_
984 }
985 else
986 {
987 QM = resEne*resEne/(beta2*2.0*electron_mas
988 QM *= (1.0-QM*0.5/electron_mass_c2);
989 }
990 if (QM < cutoffEne)
991 {
992 XHDL = G4Log(cutoffEne*(QM+2.0*electron_ma
993 *invResEne;
994 XHDT = XHDT0*invResEne;
995 }
996 else
997 {
998 QM = cutoffEne;
999 XHDL = 0.;
1000 XHDT = 0.;
1001 }
1002 }
1003 else
1004 {
1005 QM = cutoffEne;
1006 cps = 0.;
1007 cp = 0.;
1008 XHDL = 0.;
1009 XHDT = 0.;
1010 }
1011 //Close collisions (Bhabha)
1012 G4double wmaxc = kineticEnergy;
1013 G4double wcl = std::max(cutEnergy,cutoffEne
1014 G4double rcl = wcl/kineticEnergy;
1015 G4double XHC = 0.;
1016 if (wcl < wmaxc)
1017 {
1018 G4double rl1 = 1.0-rcl;
1019 XHC = ((1.0/rcl-1.0)+bha1*G4Log(rcl)+bh
1020 + (bha3/2.0)*(rcl*rcl-1.0)
1021 + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineti
1022 }
1023
1024 //Total cross section per molecule for the
1025 G4double XHTOT = XHC + XHDL + XHDT;
1026
1027 //very small cross section, do nothing
1028 if (XHTOT < 1.e-14*barn)
1029 {
1030 fKineticEnergy1=kineticEnergy;
1031 fCosThetaPrimary=1.0;
1032 fEnergySecondary=0.0;
1033 fCosThetaSecondary=1.0;
1034 fTargetOscillator = G4int(numberOfOscil
1035 return;
1036 }
1037
1038 //decide which kind of interaction we'll ha
1039 TST = XHTOT*G4UniformRand();
1040
1041 // Hard close collision
1042 G4double TS1 = XHC;
1043 if (TST < TS1)
1044 {
1045 G4double rl1 = 1.0-rcl;
1046 G4double rk=0.;
1047 G4bool loopAgain = false;
1048 do
1049 {
1050 loopAgain = false;
1051 rk = rcl/(1.0-G4UniformRand()*rl1);
1052 G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(
1053 if (G4UniformRand() > phi)
1054 loopAgain = true;
1055 }while(loopAgain);
1056 //energy and scattering angle (primary
1057 G4double deltaE = rk*kineticEnergy;
1058 fKineticEnergy1 = kineticEnergy - delta
1059 fCosThetaPrimary = std::sqrt(fKineticEn
1060 //energy and scattering angle of the de
1061 fEnergySecondary = deltaE - ionEne; //s
1062 fCosThetaSecondary= std::sqrt(deltaE*rb
1063 if (fVerboseLevel > 3)
1064 G4cout << "SampleFinalStatePositron: sample
1065 return;
1066 }
1067
1068 //Hard distant longitudinal collisions
1069 TS1 += XHDL;
1070 G4double deltaE = resEne;
1071 fKineticEnergy1 = kineticEnergy - deltaE;
1072 if (TST < TS1)
1073 {
1074 G4double QS = QM/(1.0+QM*0.5/electron_m
1075 G4double Q = QS/(std::pow((QS/cutoffEne
1076 - (QS*0.5/electron_mass_c2));
1077 G4double QTREV = Q*(Q+2.0*electron_mass
1078 G4double cpps = fKineticEnergy1*(fKinet
1079 fCosThetaPrimary = (cpps+cps-QTREV)/(2.
1080 if (fCosThetaPrimary > 1.)
1081 fCosThetaPrimary = 1.0;
1082 //energy and emission angle of the delt
1083 fEnergySecondary = deltaE - ionEne;
1084 fCosThetaSecondary = 0.5*(deltaE*(kinet
1085 if (fCosThetaSecondary > 1.0)
1086 fCosThetaSecondary = 1.0;
1087 if (fVerboseLevel > 3)
1088 G4cout << "SampleFinalStatePositron: sample
1089 return;
1090 }
1091
1092 //Hard distant transverse collisions
1093 fCosThetaPrimary = 1.0;
1094 //energy and emission angle of the delta ra
1095 fEnergySecondary = deltaE - ionEne;
1096 fCosThetaSecondary = 0.5;
1097
1098 if (fVerboseLevel > 3)
1099 G4cout << "SampleFinalStatePositron: samp
1100
1101 return;
1102 }
1103
1104 //....oooOO0OOooo........oooOO0OOooo........o
1105
1106 void G4PenelopeIonisationModel::SetParticle(c
1107 {
1108 if(!fParticle) {
1109 fParticle = p;
1110 }
1111 }
1112