119 static const G4double a = 20.0 , b = 230.0 , c = 440.0;
129 G4double p1Z = Z*(d1 + e1*Z + f1*Z*
Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
130 p3Z = Z*(d3 + e3*Z + f3*Z*
Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
133 if (Z < 1.5) { T0 = 40.0*
keV; }
136 xSection = p1Z*
G4Log(1.+2.*X)/X
137 + (p2Z + p3Z*X + p4Z*X*
X)/(1. + a*X + b*X*X + c*X*X*X);
141 if (gammaEnergy < T0) {
144 + (p2Z + p3Z*X + p4Z*X*
X)/(1. + a*X + b*X*X + c*X*X*X);
145 G4double c1 = -T0*(sigma-xSection)/(xSection*dT0);
147 if (Z > 1.5) { c2 = 0.375-0.0556*
G4Log(Z); }
149 xSection *=
G4Exp(-y*(c1+c2*y));
152 if(xSection < 0.0) { xSection = 0.0; }
161 std::vector<G4DynamicParticle*>* fvect,
183 for(i=0; i<nShells; ++i) {
208 for(i=0; i<nShells; ++i) {
if(xprob <=
fProbabilities[i]) {
break; } }
211 lv1.
set(0.0,0.0,energy,energy);
219 eKinEnergy = bindingEnergy*
x;
220 ePotEnergy = bindingEnergy*(1.0 +
x);
226 G4double sintet = sqrt((1 - costet)*(1 + costet));
227 lv2.
set(eTotMomentum*sintet*cos(phi),eTotMomentum*sintet*sin(phi),
232 gamEnergy0 =
lv1.
e();
250 G4double alpha2 = alpha1 + 0.5*(1 - epsilon0sq);
255 if(nloop > nlooplim) {
return; }
260 if ( alpha1 > alpha2*rndm[0] ) {
261 epsilon =
G4Exp(-alpha1*rndm[1]);
265 epsilonsq = epsilon0sq + (1.- epsilon0sq)*rndm[1];
266 epsilon = sqrt(epsilonsq);
269 onecost = (1.-
epsilon)/(epsilon*E0_m);
270 sint2 = onecost*(2.-onecost);
271 greject = 1. - epsilon*sint2/(1.+ epsilonsq);
274 }
while (greject < rndm[2]);
275 gamEnergy1 = epsilon*gamEnergy0;
284 if(sint2 < 0.0) { sint2 = 0.0; }
285 costet = 1. - onecost;
286 sintet = sqrt(sint2);
287 phi = twopi * rndmEngineMod->
flat();
295 lv1.
set(gamEnergy1*v.
x(),gamEnergy1*v.
y(),gamEnergy1*v.
z(),gamEnergy1);
304 }
while ( eKinEnergy < 0.0 );
311 gamEnergy1 =
lv1.
e();
331 fvect->push_back(dp);
332 }
else { eKinEnergy = 0.0; }
345 G4int nbefore = fvect->size();
347 G4int nafter = fvect->size();
349 for (
G4int j=nbefore; j<nafter; ++j) {
350 G4double e = ((*fvect)[j])->GetKineticEnergy();
351 if(esec + e > edep) {
354 ((*fvect)[j])->SetKineticEnergy(e);
367 for (
G4int jj=nafter-1; jj>j; --jj) {
378 if(std::abs(energy - gamEnergy1 - eKinEnergy - esec - edep) >
eV) {
379 G4cout <<
"### G4KleinNishinaModel dE(eV)= "
380 << (energy - gamEnergy1 - eKinEnergy - esec -
edep)/
eV
382 <<
" E(keV)= " << energy/
keV
383 <<
" Ebind(keV)= " << bindingEnergy/
keV
384 <<
" Eg(keV)= " << gamEnergy1/
keV
385 <<
" Ee(keV)= " << eKinEnergy/
keV
386 <<
" Esec(keV)= " << esec/
keV
387 <<
" Edep(keV)= " << edep/
keV
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetAtomicShell(G4int index) const
G4ParticleChangeForGamma * fParticleChange
CLHEP::Hep3Vector G4ThreeVector
std::vector< ExP01TrackerHit * > a
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4VAtomDeexcitation * AtomDeexcitation()
G4VAtomDeexcitation * fAtomDeexcitation
G4int GetNbOfAtomicShells() const
G4double lowestSecondaryEnergy
static constexpr double keV
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual void flatArray(const int size, double *vect)=0
Hep3Vector & rotateUz(const Hep3Vector &)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
std::vector< G4EmElementSelector * > * GetElementSelectors()
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4double G4Log(G4double x)
std::vector< G4double > fProbabilities
G4double LowEnergyLimit() const
virtual ~G4KleinNishinaModel()
void SetDeexcitationFlag(G4bool val)
G4KleinNishinaModel(const G4String &nam="KleinNishina")
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
static const G4double T0[78]
static constexpr double electron_mass_c2
G4ParticleDefinition * theGamma
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
static constexpr double twopi
static const G4int nlooplim
static constexpr double eV
void set(double x, double y, double z, double t)
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) override
static G4Electron * Electron()
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double cut, G4double emax) override
double epsilon(double density, double temperature)
G4ParticleDefinition * theElectron
G4int GetNbOfShellElectrons(G4int index) const
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
Hep3Vector boostVector() const
G4ParticleChangeForGamma * GetParticleChangeForGamma()
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double bindingEnergy(G4int A, G4int Z)
void ProposeTrackStatus(G4TrackStatus status)
HepLorentzVector & boost(double, double, double)
static constexpr double barn