288{
289
290
291
292
293
294
295
296
297
298
299
300
301
302
304
305 if (verboseLevel > 3) {
306 G4cout <<
"G4LivermoreComptonModel::SampleSecondaries() E(MeV)= " << photonEnergy0 / MeV
308 }
309
310
311
313
314 G4double e0m = photonEnergy0 / electron_mass_c2;
316
317
320
322
323 G4double epsilon0Local = 1. / (1. + 2. * e0m);
324 G4double epsilon0Sq = epsilon0Local * epsilon0Local;
327
328 G4double wlPhoton = h_Planck * c_light / photonEnergy0;
329
330
336
337 if (verboseLevel > 3) {
338 G4cout <<
"Started loop to sample gamma energy" <<
G4endl;
339 }
340
341 do {
345 }
346 else {
347 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) *
G4UniformRand();
348 epsilon = std::sqrt(epsilonSq);
349 }
350
352 sinT2 = oneCosT * (2. - oneCosT);
353 G4double x = std::sqrt(oneCosT / 2.) * cm / wlPhoton;
354 G4double scatteringFunction = ComputeScatteringFunction(x, Z);
355 gReject = (1. -
epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
356
358
360 G4double sinTheta = std::sqrt(sinT2);
362 G4double dirx = sinTheta * std::cos(phi);
363 G4double diry = sinTheta * std::sin(phi);
365
366
367
368
369
370
371
372 static G4int maxDopplerIterations = 1000;
378
380
381 if (verboseLevel > 3) {
383 }
384
385 do {
386 ++iteration;
387
390
391 if (verboseLevel > 3) {
392 G4cout <<
"Shell ID= " << shellIdx <<
" Ebind(keV)= " << bindingE / keV <<
G4endl;
393 }
394
395 eMax = photonEnergy0 - bindingE;
396
397
398
400 if (verboseLevel > 3) {
402 }
403
404 G4double pDoppler = pSample * fine_structure_const;
405 G4double pDoppler2 = pDoppler * pDoppler;
407 G4double var3 = var2 * var2 - pDoppler2;
408 G4double var4 = var2 - pDoppler2 * cosTheta;
409 G4double var = var4 * var4 - var3 + pDoppler2 * var3;
410 if (var > 0.) {
412 G4double scale = photonEnergy0 / var3;
413
415 photonE = (var4 - varSqrt) * scale;
416 }
417 else {
418 photonE = (var4 + varSqrt) * scale;
419 }
420 }
421 else {
422 photonE = -1.;
423 }
424 } while (iteration <= maxDopplerIterations && photonE > eMax);
425
426
427
428
429 if (iteration >= maxDopplerIterations) {
430 photonE = photonEoriginal;
431 bindingE = 0.;
432 }
433
434
436 photonDirection1.rotateUz(photonDirection0);
438
440
441 if (photonEnergy1 > 0.) {
443 }
444 else {
445
446 photonEnergy1 = 0.;
450 return;
451 }
452
453
454 G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
455
456
457 if (eKineticEnergy < 0.0) {
459 return;
460 }
461
462 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
463
464 G4double electronE = photonEnergy0 * (1. -
epsilon) + electron_mass_c2;
465 G4double electronP2 = electronE * electronE - electron_mass_c2 * electron_mass_c2;
468 if (electronP2 > 0.) {
469 cosThetaE = (eTotalEnergy + photonEnergy1) * (1. -
epsilon) / std::sqrt(electronP2);
470 sinThetaE = -1. * sqrt(1. - cosThetaE * cosThetaE);
471 }
472
473 G4double eDirX = sinThetaE * std::cos(phi);
474 G4double eDirY = sinThetaE * std::sin(phi);
476
478 eDirection.rotateUz(photonDirection0);
480 fvect->push_back(dp);
481
482
483 if (verboseLevel > 3) {
484 G4cout <<
"Started atomic de-excitation " << fAtomDeexcitation <<
G4endl;
485 }
486
487 if (nullptr != fAtomDeexcitation && iteration < maxDopplerIterations) {
490 size_t nbefore = fvect->size();
494 size_t nafter = fvect->size();
495 if (nafter > nbefore) {
496 for (size_t i = nbefore; i < nafter; ++i) {
497
498 if (bindingE >= ((*fvect)[i])->GetKineticEnergy()) {
499
500 bindingE -= ((*fvect)[i])->GetKineticEnergy();
501 }
502 else {
503
504
505 delete (*fvect)[i];
506 (*fvect)[i] = nullptr;
507 }
508 }
509 }
510 }
511 }
512 bindingE = std::max(bindingE, 0.0);
514}
G4double epsilon(G4double density, G4double temperature)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
const G4Material * GetMaterial() const
const G4String & GetName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4double BindingEnergy(G4int Z, G4int shellIndex) const
G4int SelectRandomShell(G4int Z) const
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)