307 if (verboseLevel > 3)
308 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" <<
G4endl;
329 if(!(gammaPolarization0.
isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.
mag()==0))
331 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
337 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
355 G4double epsilon0Local = 1./(1. + 2*E0_m);
356 G4double epsilon0Sq = epsilon0Local*epsilon0Local;
357 G4double alpha1 = - std::log(epsilon0Local);
358 G4double alpha2 = 0.5*(1.- epsilon0Sq);
373 epsilon = std::sqrt(epsilonSq);
376 onecost = (1.-
epsilon)/(epsilon*E0_m);
377 sinThetaSqr = onecost*(2.-onecost);
380 if (sinThetaSqr > 1.)
383 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
384 <<
"sin(theta)**2 = "
390 if (sinThetaSqr < 0.)
393 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
394 <<
"sin(theta)**2 = "
402 G4double x = std::sqrt(onecost/2.) / (wlGamma/
cm);;
404 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
413 G4double phi = SetPhi(epsilon,sinThetaSqr);
426 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
436 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
446 G4double sinTheta = std::sqrt (sinThetaSqr);
452 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
462 <<
" -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
472 G4double dirx = sinTheta*std::cos(phi);
473 G4double diry = sinTheta*std::sin(phi);
486 static G4int maxDopplerIterations = 1000;
488 G4double photonEoriginal = epsilon * gammaEnergy0;
495 if (verboseLevel > 3) {
506 if (verboseLevel > 3) {
507 G4cout <<
"Shell ID= " << shellIdx
508 <<
" Ebind(keV)= " << bindingE/
keV <<
G4endl;
510 eMax = gammaEnergy0 - bindingE;
515 if (verboseLevel > 3) {
520 G4double pDoppler2 = pDoppler * pDoppler;
521 G4double var2 = 1. + onecost * E0_m;
522 G4double var3 = var2*var2 - pDoppler2;
523 G4double var4 = var2 - pDoppler2 * cosTheta;
524 G4double var = var4*var4 - var3 + pDoppler2 * var3;
528 G4double scale = gammaEnergy0 / var3;
530 if (
G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
531 else photonE = (var4 + varSqrt) * scale;
537 }
while ( iteration <= maxDopplerIterations &&
538 (photonE < 0. || photonE > eMax || photonE < eMax*
G4UniformRand()) );
544 if (iteration >= maxDopplerIterations)
546 photonE = photonEoriginal;
550 gammaEnergy1 = photonE;
563 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
570 gammaDirection1 = tmpDirection1;
574 SystemOfRefChange(gammaDirection0,gammaDirection1,
575 gammaPolarization0,gammaPolarization1);
577 if (gammaEnergy1 > 0.)
594 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
598 if(ElecKineEnergy < 0.0) {
607 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
608 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
611 fvect->push_back(dp);
616 if (verboseLevel > 3) {
617 G4cout <<
"Started atomic de-excitation " << fAtomDeexcitation <<
G4endl;
620 if(fAtomDeexcitation && iteration < maxDopplerIterations) {
623 size_t nbefore = fvect->size();
627 size_t nafter = fvect->size();
628 if(nafter > nbefore) {
629 for (
size_t i=nbefore; i<nafter; ++i) {
631 if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
634 bindingE -= ((*fvect)[i])->GetKineticEnergy();
649 G4Exception(
"G4LivermoreComptonModel::SampleSecondaries()",
G4double LowEnergyLimit() const
static constexpr double h_Planck
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
G4double GetKineticEnergy() const
G4ParticleDefinition * GetDefinition() const
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4int SelectRandomShell(G4int Z) const
double howOrthogonal(const Hep3Vector &v) const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double electron_mass_c2
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
G4double BindingEnergy(G4int Z, G4int shellIndex) const
static constexpr double cm
void ProposePolarization(const G4ThreeVector &dir)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
virtual G4double FindValue(G4double x, G4int componentId=0) const
static constexpr double c_light
const G4ThreeVector & GetPolarization() const
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeTrackStatus(G4TrackStatus status)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
static constexpr double fine_structure_const
static constexpr double keV
double epsilon(double density, double temperature)
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)