281{
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
298
299 if (verboseLevel > 3) {
300 G4cout <<
"G4LivermoreComptonModel::SampleSecondaries() E(MeV)= "
303 }
304
305
306
308 return ;
309
310 G4double e0m = photonEnergy0 / electron_mass_c2 ;
312
313
316
318
319 G4double epsilon0Local = 1. / (1. + 2. * e0m);
320 G4double epsilon0Sq = epsilon0Local * epsilon0Local;
323
324 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
325
326
332
333 if (verboseLevel > 3) {
334 G4cout <<
"Started loop to sample gamma energy" <<
G4endl;
335 }
336
337 do {
339 {
342 }
343 else
344 {
345 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) *
G4UniformRand();
346 epsilon = std::sqrt(epsilonSq);
347 }
348
350 sinT2 = oneCosT * (2. - oneCosT);
351 G4double x = std::sqrt(oneCosT/2.) * cm / wlPhoton;
352 G4double scatteringFunction = ComputeScatteringFunction(x,
Z);
353 gReject = (1. -
epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
354
356
358 G4double sinTheta = std::sqrt (sinT2);
360 G4double dirx = sinTheta * std::cos(phi);
361 G4double diry = sinTheta * std::sin(phi);
363
364
365
366
367
368
369
370 static G4int maxDopplerIterations = 1000;
376
378
379 if (verboseLevel > 3) {
381 }
382
383 do
384 {
385 ++iteration;
386
389
390 if (verboseLevel > 3) {
391 G4cout <<
"Shell ID= " << shellIdx
392 <<
" 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.)
411 {
413 G4double scale = photonEnergy0 / var3;
414
415 if (
G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
416 else { photonE = (var4 + varSqrt) * scale; }
417 }
418 else
419 {
420 photonE = -1.;
421 }
422 } while (iteration <= maxDopplerIterations && photonE > eMax);
423
424
425
426
427 if (iteration >= maxDopplerIterations)
428 {
429 photonE = photonEoriginal;
430 bindingE = 0.;
431 }
432
433
435 photonDirection1.rotateUz(photonDirection0);
437
439
440 if (photonEnergy1 > 0.) {
442
443 } else {
444
445 photonEnergy1 = 0.;
449 return;
450 }
451
452
453 G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
454
455
456 if(eKineticEnergy < 0.0) {
458 return;
459 }
460
461 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
462
463 G4double electronE = photonEnergy0 * (1. -
epsilon) + electron_mass_c2;
465 electronE*electronE - electron_mass_c2*electron_mass_c2;
468 if (electronP2 > 0.)
469 {
470 cosThetaE = (eTotalEnergy + photonEnergy1 )*
471 (1. -
epsilon) / std::sqrt(electronP2);
472 sinThetaE = -1. * sqrt(1. - cosThetaE * cosThetaE);
473 }
474
475 G4double eDirX = sinThetaE * std::cos(phi);
476 G4double eDirY = sinThetaE * std::sin(phi);
478
480 eDirection.rotateUz(photonDirection0);
482 eDirection,eKineticEnergy) ;
483 fvect->push_back(dp);
484
485
486 if (verboseLevel > 3) {
487 G4cout <<
"Started atomic de-excitation " << fAtomDeexcitation <<
G4endl;
488 }
489
490 if(nullptr != fAtomDeexcitation && iteration < maxDopplerIterations) {
493 size_t nbefore = fvect->size();
497 size_t nafter = fvect->size();
498 if(nafter > nbefore) {
499 for (size_t i=nbefore; i<nafter; ++i) {
500
501 if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
502 {
503
504 bindingE -= ((*fvect)[i])->GetKineticEnergy();
505 }
506 else
507 {
508
509
510 delete (*fvect)[i];
511 (*fvect)[i]=0;
512 }
513 }
514 }
515 }
516 }
517 bindingE = std::max(bindingE, 0.0);
519}
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)