284{
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
301
302 if (verboseLevel > 3) {
303 G4cout <<
"G4LivermoreComptonModel::SampleSecondaries() E(MeV)= "
306 }
307
308
309
311 return ;
312
313 G4double e0m = photonEnergy0 / electron_mass_c2 ;
315
316
319
321
322 G4double epsilon0Local = 1. / (1. + 2. * e0m);
323 G4double epsilon0Sq = epsilon0Local * epsilon0Local;
325 G4double alpha2 = 0.5 * (1. - epsilon0Sq);
326
327 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
328
329
335
336 if (verboseLevel > 3) {
337 G4cout <<
"Started loop to sample gamma energy" <<
G4endl;
338 }
339
340 do {
342 {
345 }
346 else
347 {
348 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) *
G4UniformRand();
349 epsilon = std::sqrt(epsilonSq);
350 }
351
353 sinT2 = oneCosT * (2. - oneCosT);
354 G4double x = std::sqrt(oneCosT/2.) * cm / wlPhoton;
355 G4double scatteringFunction = ComputeScatteringFunction(x, Z);
356 gReject = (1. -
epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction;
357
359
361 G4double sinTheta = std::sqrt (sinT2);
363 G4double dirx = sinTheta * std::cos(phi);
364 G4double diry = sinTheta * std::sin(phi);
366
367
368
369
370
371
372
373 static G4int maxDopplerIterations = 1000;
379
381
382 if (verboseLevel > 3) {
384 }
385
386 do
387 {
388 ++iteration;
389
392
393 if (verboseLevel > 3) {
394 G4cout <<
"Shell ID= " << shellIdx
395 <<
" Ebind(keV)= " << bindingE/keV <<
G4endl;
396 }
397
398 eMax = photonEnergy0 - bindingE;
399
400
401
403 if (verboseLevel > 3) {
405 }
406
407 G4double pDoppler = pSample * fine_structure_const;
408 G4double pDoppler2 = pDoppler * pDoppler;
410 G4double var3 = var2*var2 - pDoppler2;
411 G4double var4 = var2 - pDoppler2 * cosTheta;
412 G4double var = var4*var4 - var3 + pDoppler2 * var3;
413 if (var > 0.)
414 {
416 G4double scale = photonEnergy0 / var3;
417
418 if (
G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; }
419 else { photonE = (var4 + varSqrt) * scale; }
420 }
421 else
422 {
423 photonE = -1.;
424 }
425 } while (iteration <= maxDopplerIterations && photonE > eMax);
426
427
428
429
430
431
432
433 if (iteration >= maxDopplerIterations)
434 {
435 photonE = photonEoriginal;
436 bindingE = 0.;
437 }
438
439
440
442 photonDirection1.rotateUz(photonDirection0);
444
446
447 if (photonEnergy1 > 0.) {
449
450 } else {
451
452 photonEnergy1 = 0.;
456 return;
457 }
458
459
460 G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
461
462
463 if(eKineticEnergy < 0.0) {
465 return;
466 }
467
468 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
469
470 G4double electronE = photonEnergy0 * (1. -
epsilon) + electron_mass_c2;
472 electronE*electronE - electron_mass_c2*electron_mass_c2;
475 if (electronP2 > 0.)
476 {
477 cosThetaE = (eTotalEnergy + photonEnergy1 )*
478 (1. -
epsilon) / std::sqrt(electronP2);
479 sinThetaE = -1. * sqrt(1. - cosThetaE * cosThetaE);
480 }
481
482 G4double eDirX = sinThetaE * std::cos(phi);
483 G4double eDirY = sinThetaE * std::sin(phi);
485
487 eDirection.rotateUz(photonDirection0);
489 eDirection,eKineticEnergy) ;
490 fvect->push_back(dp);
491
492
493
494
495 if (verboseLevel > 3) {
496 G4cout <<
"Started atomic de-excitation " << fAtomDeexcitation <<
G4endl;
497 }
498
499 if(fAtomDeexcitation && iteration < maxDopplerIterations) {
502 size_t nbefore = fvect->size();
506 size_t nafter = fvect->size();
507 if(nafter > nbefore) {
508 for (size_t i=nbefore; i<nafter; ++i) {
509
510 if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
511 {
512
513 bindingE -= ((*fvect)[i])->GetKineticEnergy();
514 }
515 else
516 {
517
518
519 delete (*fvect)[i];
520 (*fvect)[i]=0;
521 }
522 }
523 }
524 }
525 }
526
527 if(bindingE < 0.0)
528 G4Exception(
"G4LivermoreComptonModel::SampleSecondaries()",
530
532}
double epsilon(double density, double temperature)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
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 G4String & GetName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
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)