76 fAtomDeexcitation(0),fPIXEflag(
false),kineticEnergy1(0.*
eV),
77 cosThetaPrimary(1.0),energySecondary(0.*
eV),
78 cosThetaSecondary(0.0),targetOscillator(-1),
79 theCrossSectionHandler(0),fLocalTable(
false)
123 G4cout <<
"Calling G4PenelopeIonisationModel::Initialise()" <<
G4endl;
130 G4cout <<
"WARNING from G4PenelopeIonisationModel " <<
G4endl;
131 G4cout <<
"Atomic de-excitation module is not instantiated, so there will not be ";
133 G4cout <<
"Please make sure this is intended" <<
G4endl;
145 G4cout <<
"======================================================================" <<
G4endl;
146 G4cout <<
"The G4PenelopeIonisationModel is being used with the PIXE flag ON." <<
G4endl;
147 G4cout <<
"Atomic de-excitation will be produced statistically by the PIXE " <<
G4endl;
148 G4cout <<
"interface by using the shell cross section --> " << theModel <<
G4endl;
149 G4cout <<
"The built-in model procedure for atomic de-excitation is disabled. " <<
G4endl;
150 G4cout <<
"*Please be sure this is intended*, or disable PIXE by" <<
G4endl;
160 G4cout <<
"======================================================================" <<
G4endl;
200 G4cout <<
"Penelope Ionisation model v2008 is initialized " << G4endl
223 G4cout <<
"Calling G4PenelopeIonisationModel::InitialiseLocal()" <<
G4endl;
275 G4cout <<
"Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" <<
G4endl;
304 ed <<
"Unable to retrieve the cross section table for " << theParticle->
GetParticleName() <<
305 " in " << material->
GetName() <<
", cut = " << cutEnergy/
keV <<
" keV " <<
G4endl;
306 ed <<
"This can happen only in Unit Tests or via G4EmCalculator" <<
G4endl;
307 G4Exception(
"G4PenelopeIonisationModel::CrossSectionPerVolume()",
311 G4AutoLock lock(&PenelopeIonisationModelMutex);
329 G4cout <<
"Material " << material->
GetName() <<
" has " << atPerMol <<
330 "atoms per molecule" <<
G4endl;
334 moleculeDensity = atomDensity/atPerMol;
336 G4double crossPerVolume = crossPerMolecule*moleculeDensity;
341 G4cout <<
"Mean free path for delta emission > " << cutEnergy/
keV <<
" keV at " <<
342 energy/
keV <<
" keV = " << (1./crossPerVolume)/
mm <<
" mm" << G4endl;
345 G4cout <<
"Total free path for ionisation (no threshold) at " <<
346 energy/
keV <<
" keV = " << (1./totalCross)/
mm <<
" mm" << G4endl;
348 return crossPerVolume;
363 G4cout <<
"*** G4PenelopeIonisationModel -- WARNING ***" <<
G4endl;
364 G4cout <<
"Penelope Ionisation model v2008 does not calculate cross section _per atom_ " <<
G4endl;
365 G4cout <<
"so the result is always zero. For physics values, please invoke " <<
G4endl;
366 G4cout <<
"GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" <<
G4endl;
393 G4cout <<
"Calling ComputeDEDX() of G4PenelopeIonisationModel" <<
G4endl;
417 ed <<
"Unable to retrieve the cross section table for " << theParticle->
GetParticleName() <<
418 " in " << material->
GetName() <<
", cut = " << cutEnergy/
keV <<
" keV " <<
G4endl;
419 ed <<
"This can happen only in Unit Tests or via G4EmCalculator" <<
G4endl;
420 G4Exception(
"G4PenelopeIonisationModel::ComputeDEDXPerVolume()",
424 G4AutoLock lock(&PenelopeIonisationModelMutex);
444 moleculeDensity = atomDensity/atPerMol;
446 G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
451 G4cout <<
"Stopping power < " << cutEnergy/
keV <<
" keV at " <<
452 kineticEnergy/
keV <<
" keV = " <<
453 sPowerPerVolume/(
keV/
mm) <<
" keV/mm" << G4endl;
455 return sPowerPerVolume;
506 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeIonisationModel" <<
G4endl;
539 G4Exception(
"G4PenelopeIonisationModel::SamplingSecondaries()",
547 G4cout <<
"G4PenelopeIonisationModel::SamplingSecondaries() for " <<
559 G4double dirx = sint * std::cos(phiPrimary);
560 G4double diry = sint * std::sin(phiPrimary);
564 electronDirection1.
rotateUz(particleDirection0);
576 G4double ionEnergyInPenelopeDatabase =
591 if (Z > 0 && shFlag<30)
593 shell = transitionManager->
Shell(Z,shFlag-1);
628 size_t nBefore = fvect->size();
630 size_t nAfter = fvect->size();
634 for (
size_t j=nBefore;j<nAfter;j++)
636 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
637 if (itsEnergy < localEnergyDeposit)
639 localEnergyDeposit -= itsEnergy;
641 energyInFluorescence += itsEnergy;
643 energyInAuger += itsEnergy;
648 (*fvect)[j] =
nullptr;
662 G4double xEl = sinThetaE * std::cos(phiEl);
663 G4double yEl = sinThetaE * std::sin(phiEl);
666 eDirection.
rotateUz(particleDirection0);
669 fvect->push_back(electron);
677 if (localEnergyDeposit < 0)
679 G4Exception(
"G4PenelopeIonisationModel::SampleSecondaries()",
680 "em2099",
JustWarning,
"WARNING: Negative local energy deposit");
681 localEnergyDeposit=0.;
687 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
688 G4cout <<
"Energy balance from G4PenelopeIonisation" <<
G4endl;
689 G4cout <<
"Incoming primary energy: " << kineticEnergy0/
keV <<
" keV" <<
G4endl;
690 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
693 if (energyInFluorescence)
694 G4cout <<
"Fluorescence x-rays: " << energyInFluorescence/
keV <<
" keV" <<
G4endl;
696 G4cout <<
"Auger electrons: " << energyInAuger/
keV <<
" keV" <<
G4endl;
697 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
699 localEnergyDeposit+energyInAuger)/
keV <<
701 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
707 localEnergyDeposit+energyInAuger-kineticEnergy0);
708 if (energyDiff > 0.05*
keV)
709 G4cout <<
"Warning from G4PenelopeIonisation: problem with energy conservation: " <<
711 " keV (final) vs. " <<
712 kineticEnergy0/
keV <<
" keV (initial)" << G4endl;
731 size_t numberOfOscillators = theTable->size();
742 for (
size_t i=0;i<numberOfOscillators-1;i++)
787 if (resEne > cutEnergy && resEne < kineticEnergy)
789 cps = kineticEnergy*
rb;
792 if (resEne > 1.0
e-6*kineticEnergy)
806 XHDT = XHDT0*invResEne;
825 G4double EE = kineticEnergy + ionEne;
834 XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
835 (1.0-amol)*std::log(rcl*rrl1))/EE;
842 if (XHTOT < 1.
e-14*
barn)
870 rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));
872 rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
888 G4cout <<
"SampleFinalStateElectron: sampled close collision " <<
G4endl;
905 if (cosThetaPrimary > 1.)
906 cosThetaPrimary = 1.0;
909 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
913 G4cout <<
"SampleFinalStateElectron: sampled distant longitudinal collision " <<
G4endl;
923 G4cout <<
"SampleFinalStateElectron: sampled distant transverse collision " <<
G4endl;
944 size_t numberOfOscillators = theTable->size();
954 for (
size_t i=0;i<numberOfOscillators-1;i++)
981 G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0);
983 G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
984 G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12;
1001 if (resEne > cutEnergy && resEne < kineticEnergy)
1003 cps = kineticEnergy*
rb;
1004 cp = std::sqrt(cps);
1006 if (resEne > 1.0
e-6*kineticEnergy)
1020 XHDT = XHDT0*invResEne;
1045 XHC = ((1.0/rcl-1.0)+bha1*std::log(rcl)+bha2*rl1
1046 + (bha3/2.0)*(rcl*rcl-1.0)
1047 + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy;
1051 G4double XHTOT = XHC + XHDL + XHDT;
1054 if (XHTOT < 1.
e-14*
barn)
1073 G4bool loopAgain =
false;
1078 G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1083 G4double deltaE = rk*kineticEnergy;
1090 G4cout <<
"SampleFinalStatePositron: sampled close collision " <<
G4endl;
1106 if (cosThetaPrimary > 1.)
1107 cosThetaPrimary = 1.0;
1110 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
1114 G4cout <<
"SampleFinalStatePositron: sampled distant longitudinal collision " <<
G4endl;
1125 G4cout <<
"SampleFinalStatePositron: sampled distant transverse collision " <<
G4endl;