189 if (verboseLevel > 3) {
190 G4cout <<
"G4LivermoreComptonModifiedModel::SampleSecondaries() E(MeV)= "
208 G4double epsilon0Local = 1. / (1. + 2. * e0m);
209 G4double epsilon0Sq = epsilon0Local * epsilon0Local;
210 G4double alpha1 = -std::log(epsilon0Local);
211 G4double alpha2 = 0.5 * (1. - epsilon0Sq);
232 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) *
G4UniformRand();
233 epsilon = std::sqrt(epsilonSq);
236 oneCosT = (1. -
epsilon) / ( epsilon * e0m);
237 sinT2 = oneCosT * (2. - oneCosT);
238 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/
cm);
240 gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
245 G4double sinTheta = std::sqrt (sinT2);
247 G4double dirx = sinTheta * std::cos(phi);
248 G4double diry = sinTheta * std::sin(phi);
257 G4int maxDopplerIterations = 1000;
259 G4double photonEoriginal = epsilon * photonEnergy0;
266 G4double momentum_au_to_nat = 1.992851740*std::pow(10.,-24.);
267 G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.);
292 }
while(Alpha >= (
pi/2.0));
294 ePAU = pSample / std::cos(Alpha);
298 G4double ePSI = ePAU * momentum_au_to_nat;
299 G4double u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c;
303 systemE = eEIncident+photonEnergy0;
307 G4double pDoppler2 = pDoppler * pDoppler;
309 G4double var3 = var2*var2 - pDoppler2;
310 G4double var4 = var2 - pDoppler2 * cosTheta;
311 G4double var = var4*var4 - var3 + pDoppler2 * var3;
315 G4double scale = photonEnergy0 / var3;
317 if (
G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
318 else { photonE = (var4 + varSqrt) * scale; }
324 }
while ( iteration <= maxDopplerIterations &&
325 (photonE < 0. || photonE > eMax ) );
336 if(eKineticEnergy < 0.0) {
337 G4cout <<
"Error, kinetic energy of electron less than zero" <<
G4endl;
344 G4double E_num = photonEnergy0 - photonE*cosTheta;
345 G4double E_dom = sqrt(photonEnergy0*photonEnergy0 + photonE*photonE -2*photonEnergy0*photonE*cosTheta);
347 G4double sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE));
349 eDirX = sinThetaE * std::cos(phi);
350 eDirY = sinThetaE * std::sin(phi);
354 eDirection.rotateUz(photonDirection0);
356 eDirection,eKineticEnergy) ;
357 fvect->push_back(dp);
363 if (iteration >= maxDopplerIterations)
365 photonE = photonEoriginal;
372 photonDirection1.rotateUz(photonDirection0);
377 if (photonEnergy1 > 0.)
381 if (iteration < maxDopplerIterations)
384 eDirection.rotateUz(photonDirection0);
386 eDirection,eKineticEnergy) ;
387 fvect->push_back(dp);
399 if(fAtomDeexcitation && iteration < maxDopplerIterations) {
402 size_t nbefore = fvect->size();
406 size_t nafter = fvect->size();
407 if(nafter > nbefore) {
408 for (
size_t i=nbefore; i<nafter; ++i) {
409 bindingE -= ((*fvect)[i])->GetKineticEnergy();
414 if(bindingE < 0.0) { bindingE = 0.0; }
G4double LowEnergyLimit() const
static constexpr double h_Planck
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
G4double GetKineticEnergy() const
const G4String & GetName() const
G4ParticleChangeForGamma * fParticleChange
G4ParticleDefinition * GetDefinition() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4int SelectRandomShell(G4int Z) const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
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
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static constexpr double c_light
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static constexpr double MeV
static constexpr double pi
void ProposeTrackStatus(G4TrackStatus status)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
static constexpr double fine_structure_const
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)
const G4Material * GetMaterial() const