Loading...
Searching...
No Matches
G4PenelopeComptonModel Class Reference
#include <G4PenelopeComptonModel.hh>
Inheritance diagram for G4PenelopeComptonModel:
Protected Attributes | |
| G4ParticleChangeForGamma * | fParticleChange |
| const G4ParticleDefinition * | fParticle |
| Protected Attributes inherited from G4VEmModel | |
| G4ElementData * | fElementData = nullptr |
| G4VParticleChange * | pParticleChange = nullptr |
| G4PhysicsTable * | xSectionTable = nullptr |
| const G4Material * | pBaseMaterial = nullptr |
| const std::vector< G4double > * | theDensityFactor = nullptr |
| const std::vector< G4int > * | theDensityIdx = nullptr |
| G4double | inveplus |
| G4double | pFactor = 1.0 |
| std::size_t | currentCoupleIndex = 0 |
| std::size_t | basedCoupleIndex = 0 |
| G4bool | lossFlucFlag = true |
Additional Inherited Members | |
| Protected Member Functions inherited from G4VEmModel | |
| G4ParticleChangeForLoss * | GetParticleChangeForLoss () |
| G4ParticleChangeForGamma * | GetParticleChangeForGamma () |
| virtual G4double | MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy) |
| const G4MaterialCutsCouple * | CurrentCouple () const |
| void | SetCurrentElement (const G4Element *) |
Detailed Description
Definition at line 62 of file G4PenelopeComptonModel.hh.
Constructor & Destructor Documentation
◆ G4PenelopeComptonModel() [1/2]
|
explicit |
Definition at line 63 of file G4PenelopeComptonModel.cc.
66 fAtomDeexcitation(nullptr),
67 fOscManager(nullptr),fIsInitialised(false)
68{
69 fIntrinsicLowEnergyLimit = 100.0*eV;
70 fIntrinsicHighEnergyLimit = 100.0*GeV;
71 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
72 //
73 fOscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
74
75 if (part)
76 SetParticle(part);
77
78 fVerboseLevel= 0;
79 // Verbosity scale:
80 // 0 = nothing
81 // 1 = warning for energy non-conservation
82 // 2 = details of energy budget
83 // 3 = calculation of cross sections, file openings, sampling of atoms
84 // 4 = entering in methods
85
86 //Mark this model as "applicable" for atomic deexcitation
88
89 fTransitionManager = G4AtomicTransitionManager::Instance();
90}
static G4AtomicTransitionManager * Instance()
Definition G4AtomicTransitionManager.cc:49
G4ParticleChangeForGamma * fParticleChange
Definition G4PenelopeComptonModel.hh:103
const G4ParticleDefinition * fParticle
Definition G4PenelopeComptonModel.hh:104
static G4PenelopeOscillatorManager * GetOscillatorManager()
Definition G4PenelopeOscillatorManager.cc:91
Referenced by G4PenelopeComptonModel(), InitialiseLocal(), and operator=().
◆ ~G4PenelopeComptonModel()
|
virtual |
Definition at line 94 of file G4PenelopeComptonModel.cc.
95{;}
◆ G4PenelopeComptonModel() [2/2]
|
delete |
Member Function Documentation
◆ ComputeCrossSectionPerAtom()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 249 of file G4PenelopeComptonModel.cc.
◆ CrossSectionPerVolume()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 176 of file G4PenelopeComptonModel.cc.
181{
182 // Penelope model v2008 to calculate the Compton scattering cross section:
183 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
184 //
185 // The cross section for Compton scattering is calculated according to the Klein-Nishina
186 // formula for energy > 5 MeV.
187 // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega,
188 // which is integrated numerically in cos(theta), G4PenelopeComptonModel::DifferentialCrossSection().
189 // The parametrization includes the J(p)
190 // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations
191 // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
192 //
193 if (fVerboseLevel > 3)
195 SetupForMaterial(p, material, energy);
196
197 G4double cs = 0;
198 //Force null cross-section if below the low-energy edge of the table
200 return cs;
201
202 //Retrieve the oscillator table for this material
203 G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableCompton(material);
204
205 if (energy < 5*MeV) //explicit calculation for E < 5 MeV
206 {
207 size_t numberOfOscillators = theTable->size();
208 for (size_t i=0;i<numberOfOscillators;i++)
209 {
210 G4PenelopeOscillator* theOsc = (*theTable)[i];
211 //sum contributions coming from each oscillator
212 cs += OscillatorTotalCrossSection(energy,theOsc);
213 }
214 }
215 else //use Klein-Nishina for E>5 MeV
216 cs = KleinNishinaCrossSection(energy,material);
217
218 //cross sections are in units of pi*classic_electr_radius^2
219 cs *= pi*classic_electr_radius*classic_electr_radius;
220
221 //Now, cs is the cross section *per molecule*, let's calculate the
222 //cross section per volume
224 G4double atPerMol = fOscManager->GetAtomsPerMolecule(material);
225
226 if (fVerboseLevel > 3)
229
230 G4double moleculeDensity = 0.;
231
232 if (atPerMol)
233 moleculeDensity = atomDensity/atPerMol;
234
235 G4double csvolume = cs*moleculeDensity;
236
237 if (fVerboseLevel > 2)
241 return csvolume;
242}
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
Definition G4PenelopeOscillatorManager.hh:65
G4double GetTotNbOfAtomsPerVolume() const
Definition G4Material.hh:197
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition G4VEmModel.cc:380
G4double energy(const ThreeVector &p, const G4double m)
Definition G4INCLKinematicsUtils.cc:185
◆ GetVerbosityLevel()
|
inline |
Definition at line 97 of file G4PenelopeComptonModel.hh.
97{return fVerboseLevel;};
◆ Initialise()
|
overridevirtual |
Implements G4VEmModel.
Definition at line 99 of file G4PenelopeComptonModel.cc.
101{
102 if (fVerboseLevel > 3)
104
106 //Issue warning if the AtomicDeexcitation has not been declared
107 if (!fAtomDeexcitation)
108 {
114 }
115
116 SetParticle(part);
117
119 {
120
121 if (fVerboseLevel > 0)
122 {
124 << "Energy range: "
127 }
128 //Issue a warning, if the model is going to be used down to a
129 //energy which is outside the validity of the model itself
131 {
132 G4ExceptionDescription ed;
133 ed << "Using the Penelope Compton model outside its intrinsic validity range. "
134 << G4endl;
136 ed << "-> Instrinsic low-energy limit = " << fIntrinsicLowEnergyLimit/keV << "keV "
137 << G4endl;
141 }
142 }
143
144 if(fIsInitialised) return;
146 fIsInitialised = true;
147
148}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
static G4LossTableManager * Instance()
Definition G4LossTableManager.cc:90
G4VAtomDeexcitation * AtomDeexcitation()
Definition G4LossTableManager.hh:396
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition G4VEmModel.cc:121
◆ InitialiseLocal()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 152 of file G4PenelopeComptonModel.cc.
154{
155 if (fVerboseLevel > 3)
157 //
158 //Check that particle matches: one might have multiple master models (e.g.
159 //for e+ and e-).
160 //
162 {
163 //Get the const table pointers from the master to the workers
166
167 //Same verbosity for all workers, as the master
168 fVerboseLevel = theModel->fVerboseLevel;
169 }
170 return;
171}
G4PenelopeComptonModel(const G4ParticleDefinition *p=nullptr, const G4String &processName="PenCompton")
Definition G4PenelopeComptonModel.cc:63
◆ operator=()
|
delete |
◆ SampleSecondaries()
|
overridevirtual |
Implements G4VEmModel.
Definition at line 265 of file G4PenelopeComptonModel.cc.
270{
271 // Penelope model v2008 to sample the Compton scattering final state.
272 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
273 // The model determines also the original shell from which the electron is expelled,
274 // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
275 //
276 // The final state for Compton scattering is calculated according to the Klein-Nishina
277 // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and
278 // one can assume that the target electron is at rest.
279 // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega,
280 // to sample the scattering angle and the energy of the emerging electron, which is
281 // G4PenelopeComptonModel::DifferentialCrossSection(). The rejection method is
282 // used to sample cos(theta). The efficiency increases monotonically with photon energy and is
283 // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV,
284 // respectively.
285 // The parametrization includes the J(p) distribution profiles for the atomic shells, that are
286 // tabulated
287 // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201.
288 // Doppler broadening is included.
289 //
290
291 if (fVerboseLevel > 3)
293
295
296 // do nothing below the threshold
297 // should never get here because the XS is zero below the limit
299 return;
300
303
304 G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableCompton(material);
305
307 G4double rn[nmax]={0.0};
308 G4double pac[nmax]={0.0};
309
312 G4double cosTheta = 1.0;
313 G4double hartreeFunc = 0.0;
314 G4double oscStren = 0.0;
315 size_t numberOfOscillators = theTable->size();
316 size_t targetOscillator = 0;
317 G4double ionEnergy = 0.0*eV;
318
319 G4double ek = photonEnergy0/electron_mass_c2;
320 G4double ek2 = 2.*ek+1.0;
321 G4double eks = ek*ek;
322 G4double ek1 = eks-ek2-1.0;
323
324 G4double taumin = 1.0/ek2;
325 //This is meant to fix a possible (rare) floating point exception in the sampling of tau below,
326 //causing an infinite loop. The maximum of tau is not 1., but the closest double which can
327 //be represented (i.e. ~ 1. - 1e-16). Fix by Domenico Iuso
329 if (fVerboseLevel > 3)
331 //To here.
333 G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2);
334
335 G4double TST = 0;
336 G4double tau = 0.;
337
338 //If the incoming photon is above 5 MeV, the quicker approach based on the
339 //pure Klein-Nishina formula is used
340 if (photonEnergy0 > 5*MeV)
341 {
342 do{
343 do{
345 tau = std::pow(taumin,G4UniformRand());
346 else
347 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
348 //rejection function
349 TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
351 if (tau > taumax) tau = taumax; //prevent FP exception causing infinite loop
352 epsilon=tau;
353 cosTheta = 1.0 - (1.0-tau)/(ek*tau);
354
355 //Target shell electrons
356 TST = fOscManager->GetTotalZ(material)*G4UniformRand();
357 targetOscillator = numberOfOscillators-1; //last level
358 S=0.0;
360 for (size_t j=0;j<numberOfOscillators && !levelFound; j++)
361 {
362 S += (*theTable)[j]->GetOscillatorStrength();
364 {
365 targetOscillator = j;
366 levelFound = true;
367 }
368 }
369 //check whether the level is valid
370 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
372 }
373 else //photonEnergy0 < 5 MeV
374 {
375 //Incoherent scattering function for theta=PI
377 G4double pzomc=0.0;
378 G4double rni=0.0;
379 G4double aux=0.0;
380 for (size_t i=0;i<numberOfOscillators;i++)
381 {
382 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
383 if (photonEnergy0 > ionEnergy)
384 {
385 G4double aux2 = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
386 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
387 oscStren = (*theTable)[i]->GetOscillatorStrength();
388 pzomc = hartreeFunc*(aux2-electron_mass_c2*ionEnergy)/
389 (electron_mass_c2*std::sqrt(2.0*aux2+ionEnergy*ionEnergy));
390 if (pzomc > 0)
391 rni = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
392 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
393 else
394 rni = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
395 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
396 s0 += oscStren*rni;
397 }
398 }
399 //Sampling tau
400 G4double cdt1 = 0.;
401 do
402 {
404 tau = std::pow(taumin,G4UniformRand());
405 else
406 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
407 if (tau > taumax) tau = taumax; //prevent FP exception causing infinite loop
408 cdt1 = (1.0-tau)/(ek*tau);
409 //Incoherent scattering function
410 S = 0.;
411 for (size_t i=0;i<numberOfOscillators;i++)
412 {
413 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
414 if (photonEnergy0 > ionEnergy) //sum only on excitable levels
415 {
416 aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
417 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
418 oscStren = (*theTable)[i]->GetOscillatorStrength();
419 pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/
420 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
421 if (pzomc > 0)
422 rn[i] = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
423 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
424 else
425 rn[i] = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
426 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
427 S += oscStren*rn[i];
428 pac[i] = S;
429 }
430 else
431 pac[i] = S-1e-6;
432 }
433 //Rejection function
434 TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
436
437 cosTheta = 1.0 - cdt1;
438 G4double fpzmax=0.0,fpz=0.0;
440 //Target electron shell
441 do
442 {
443 do
444 {
446 targetOscillator = numberOfOscillators-1; //last level
448 for (size_t i=0;i<numberOfOscillators && !levelFound;i++)
449 {
450 if (pac[i]>TST)
451 {
452 targetOscillator = i;
453 levelFound = true;
454 }
455 }
457 hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor();
458 oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength();
461 (std::sqrt(2.0)*hartreeFunc);
462 else
464 (std::sqrt(2.0)*hartreeFunc);
465 } while (pzomc < -1);
466
467 // F(EP) rejection
468 G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
469 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
470 if (AF > 0)
471 fpzmax = 1.0+AF*0.2;
472 else
473 fpzmax = 1.0-AF*0.2;
474 fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
476
477 //Energy of the scattered photon
478 G4double T = pzomc*pzomc;
479 G4double b1 = 1.0-T*tau*tau;
480 G4double b2 = 1.0-T*tau*cosTheta;
481 if (pzomc > 0.0)
482 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
483 else
484 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
485 } //energy < 5 MeV
486
487 //Ok, the kinematics has been calculated.
488 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
490 G4double dirx = sinTheta * std::cos(phi);
491 G4double diry = sinTheta * std::sin(phi);
492 G4double dirz = cosTheta ;
493
494 // Update G4VParticleChange for the scattered photon
495 G4ThreeVector photonDirection1(dirx,diry,dirz);
496 photonDirection1.rotateUz(photonDirection0);
497 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
498
500
501 if (photonEnergy1 > 0.)
502 fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
503 else
504 {
505 fParticleChange->SetProposedKineticEnergy(0.) ;
507 }
508
509 //Create scattered electron
511 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
512
513 G4double Q2 =
514 photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
516
517 if (Q2 > 1.0e-12)
518 cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
519 else
520 cosThetaE = 1.0;
521 G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
522
523 //Now, try to handle fluorescence
524 //Notice: merged levels are indicated with Z=0 and flag=30
525 G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
527
528 //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
531
532 //Real level
533 if (Z > 0 && shFlag<30)
534 {
535 shell = fTransitionManager->Shell(Z,shFlag-1);
537 }
538
539 G4double ionEnergyInPenelopeDatabase = ionEnergy;
540 //protection against energy non-conservation
541 ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase);
542
543 //subtract the excitation energy. If not emitted by fluorescence
544 //the ionization energy is deposited as local energy deposition
545 G4double eKineticEnergy = diffEnergy - ionEnergy;
546 G4double localEnergyDeposit = ionEnergy;
549
550 if (eKineticEnergy < 0)
551 {
552 //It means that there was some problem/mismatch between the two databases.
553 //Try to make it work
554 //In this case available Energy (diffEnergy) < ionEnergy
555 //Full residual energy is deposited locally
556 localEnergyDeposit = diffEnergy;
557 eKineticEnergy = 0.0;
558 }
559
560 //the local energy deposit is what remains: part of this may be spent for fluorescence.
561 //Notice: shell might be nullptr (invalid!) if shFlag=30. Must be protected
562 //Now, take care of fluorescence, if required
563 if (fAtomDeexcitation && shell)
564 {
566 if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index))
567 {
568 size_t nBefore = fvect->size();
569 fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
570 size_t nAfter = fvect->size();
571
572 if (nAfter > nBefore) //actual production of fluorescence
573 {
574 for (size_t j=nBefore;j<nAfter;j++) //loop on products
575 {
576 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
577 if (itsEnergy < localEnergyDeposit) // valid secondary, generate it
578 {
579 localEnergyDeposit -= itsEnergy;
581 energyInFluorescence += itsEnergy;
582 else if (((*fvect)[j])->GetParticleDefinition() ==
583 G4Electron::Definition())
584 energyInAuger += itsEnergy;
585 }
586 else //invalid secondary: takes more than the available energy: delete it
587 {
588 delete (*fvect)[j];
589 (*fvect)[j] = nullptr;
590 }
591 }
592 }
593
594 }
595 }
596
597 //Always produce explicitly the electron
598 G4DynamicParticle* electron = 0;
599
600 G4double xEl = sinThetaE * std::cos(phi+pi);
601 G4double yEl = sinThetaE * std::sin(phi+pi);
602 G4double zEl = cosThetaE;
604 eDirection.rotateUz(photonDirection0);
606 eDirection,eKineticEnergy) ;
607 fvect->push_back(electron);
608
609 if (localEnergyDeposit < 0) //Should not be: issue a G4Exception (warning)
610 {
613 localEnergyDeposit=0.;
614 }
615 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
616
617 G4double electronEnergy = 0.;
618 if (electron)
619 electronEnergy = eKineticEnergy;
620 if (fVerboseLevel > 1)
621 {
628 if (energyInFluorescence)
630 if (energyInAuger)
634 localEnergyDeposit+energyInAuger)/keV <<
637 }
638 if (fVerboseLevel > 0)
639 {
640 G4double energyDiff = std::fabs(photonEnergy1+
641 electronEnergy+energyInFluorescence+
642 localEnergyDeposit+energyInAuger-photonEnergy0);
643 if (energyDiff > 0.05*keV)
645 (photonEnergy1+electronEnergy+energyInFluorescence+energyInAuger+localEnergyDeposit)/keV <<
646 " keV (final) vs. " <<
648 }
649}
G4double epsilon(G4double density, G4double temperature)
Definition G4DNAElectronHoleRecombination.cc:98
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
const G4Material * GetMaterial() const
Definition G4MaterialCutsCouple.hh:166
G4double bindingEnergy(G4int A, G4int Z)
Definition bindingEnergy.cc:36
◆ SetVerbosityLevel()
|
inline |
Definition at line 96 of file G4PenelopeComptonModel.hh.
96{fVerboseLevel = lev;};
Member Data Documentation
◆ fParticle
|
protected |
Definition at line 104 of file G4PenelopeComptonModel.hh.
Referenced by G4PenelopeComptonModel(), Initialise(), and InitialiseLocal().
◆ fParticleChange
|
protected |
Definition at line 103 of file G4PenelopeComptonModel.hh.
Referenced by G4PenelopeComptonModel(), Initialise(), and SampleSecondaries().
The documentation for this class was generated from the following files:
