75 :
G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false),
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)
81 fIntrinsicLowEnergyLimit = 100.0*
eV;
82 fIntrinsicHighEnergyLimit = 100.0*
GeV;
112 if (theCrossSectionHandler)
113 delete theCrossSectionHandler;
122 if (verboseLevel > 3)
123 G4cout <<
"Calling G4PenelopeIonisationModel::Initialise()" <<
G4endl;
127 if (!fAtomDeexcitation)
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;
136 if (fAtomDeexcitation)
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;
166 SetParticle(particle);
175 nBins =
std::max(nBins,(
size_t)100);
178 if (theCrossSectionHandler)
180 delete theCrossSectionHandler;
181 theCrossSectionHandler = 0;
194 theCrossSectionHandler->
BuildXSTable(theMat,theCuts.at(i),particle,
199 if (verboseLevel > 2) {
200 G4cout <<
"Penelope Ionisation model v2008 is initialized " << G4endl
212 isInitialised =
true;
222 if (verboseLevel > 3)
223 G4cout <<
"Calling G4PenelopeIonisationModel::InitialiseLocal()" <<
G4endl;
236 theCrossSectionHandler = theModel->theCrossSectionHandler;
239 nBins = theModel->nBins;
242 verboseLevel = theModel->verboseLevel;
274 if (verboseLevel > 3)
275 G4cout <<
"Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" <<
G4endl;
284 if (!theCrossSectionHandler)
300 if (verboseLevel > 0)
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);
312 theCrossSectionHandler->
BuildXSTable(material,cutEnergy,theParticle);
328 if (verboseLevel > 3)
329 G4cout <<
"Material " << material->
GetName() <<
" has " << atPerMol <<
330 "atoms per molecule" <<
G4endl;
334 moleculeDensity = atomDensity/atPerMol;
336 G4double crossPerVolume = crossPerMolecule*moleculeDensity;
338 if (verboseLevel > 2)
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;
392 if (verboseLevel > 3)
393 G4cout <<
"Calling ComputeDEDX() of G4PenelopeIonisationModel" <<
G4endl;
397 if (!theCrossSectionHandler)
413 if (verboseLevel > 0)
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);
425 theCrossSectionHandler->
BuildXSTable(material,cutEnergy,theParticle);
444 moleculeDensity = atomDensity/atPerMol;
446 G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
448 if (verboseLevel > 2)
451 G4cout <<
"Stopping power < " << cutEnergy/
keV <<
" keV at " <<
452 kineticEnergy/
keV <<
" keV = " <<
453 sPowerPerVolume/(
keV/
mm) <<
" keV/mm" << G4endl;
455 return sPowerPerVolume;
463 return fIntrinsicLowEnergyLimit;
505 if (verboseLevel > 3)
506 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeIonisationModel" <<
G4endl;
511 if (kineticEnergy0 <= fIntrinsicLowEnergyLimit)
525 kineticEnergy1=kineticEnergy0;
528 cosThetaSecondary=1.0;
529 targetOscillator = -1;
532 SampleFinalStateElectron(material,cutE,kineticEnergy0);
534 SampleFinalStatePositron(material,cutE,kineticEnergy0);
539 G4Exception(
"G4PenelopeIonisationModel::SamplingSecondaries()",
543 if (energySecondary == 0)
return;
545 if (verboseLevel > 3)
547 G4cout <<
"G4PenelopeIonisationModel::SamplingSecondaries() for " <<
549 G4cout <<
"Final eKin = " << kineticEnergy1 <<
" keV" <<
G4endl;
550 G4cout <<
"Final cosTheta = " << cosThetaPrimary <<
G4endl;
551 G4cout <<
"Delta-ray eKin = " << energySecondary <<
" keV" <<
G4endl;
552 G4cout <<
"Delta-ray cosTheta = " << cosThetaSecondary <<
G4endl;
557 G4double sint = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
559 G4double dirx = sint * std::cos(phiPrimary);
560 G4double diry = sint * std::sin(phiPrimary);
564 electronDirection1.
rotateUz(particleDirection0);
566 if (kineticEnergy1 > 0)
576 G4double ionEnergyInPenelopeDatabase =
577 (*theTable)[targetOscillator]->GetIonisationEnergy();
581 G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
582 G4int Z = (
G4int) (*theTable)[targetOscillator]->GetParentZ();
591 if (Z > 0 && shFlag<30)
593 shell = transitionManager->
Shell(Z,shFlag-1);
601 energySecondary += ionEnergyInPenelopeDatabase-
bindingEnergy;
608 if (energySecondary < 0)
614 localEnergyDeposit += energySecondary;
615 energySecondary = 0.0;
623 if (fAtomDeexcitation && !fPIXEflag && shell)
628 size_t nBefore = fvect->size();
630 size_t nAfter = fvect->size();
632 if (nAfter > nBefore)
634 for (
size_t j=nBefore;j<nAfter;j++)
636 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
637 localEnergyDeposit -= itsEnergy;
639 energyInFluorescence += itsEnergy;
641 energyInAuger += itsEnergy;
648 if (energySecondary > cutE)
651 G4double sinThetaE = std::sqrt(1.-cosThetaSecondary*cosThetaSecondary);
653 G4double xEl = sinThetaE * std::cos(phiEl);
654 G4double yEl = sinThetaE * std::sin(phiEl);
657 eDirection.
rotateUz(particleDirection0);
659 eDirection,energySecondary) ;
660 fvect->push_back(electron);
664 localEnergyDeposit += energySecondary;
668 if (localEnergyDeposit < 0)
671 <<
"G4PenelopeIonisationModel::SampleSecondaries - Negative energy deposit"
673 localEnergyDeposit=0.;
677 if (verboseLevel > 1)
679 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
680 G4cout <<
"Energy balance from G4PenelopeIonisation" <<
G4endl;
681 G4cout <<
"Incoming primary energy: " << kineticEnergy0/
keV <<
" keV" <<
G4endl;
682 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
683 G4cout <<
"Outgoing primary energy: " << kineticEnergy1/
keV <<
" keV" <<
G4endl;
685 if (energyInFluorescence)
686 G4cout <<
"Fluorescence x-rays: " << energyInFluorescence/
keV <<
" keV" <<
G4endl;
688 G4cout <<
"Auger electrons: " << energyInAuger/
keV <<
" keV" <<
G4endl;
689 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
690 G4cout <<
"Total final state: " << (energySecondary+energyInFluorescence+kineticEnergy1+
691 localEnergyDeposit+energyInAuger)/
keV <<
693 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
696 if (verboseLevel > 0)
698 G4double energyDiff = std::fabs(energySecondary+energyInFluorescence+kineticEnergy1+
699 localEnergyDeposit+energyInAuger-kineticEnergy0);
700 if (energyDiff > 0.05*
keV)
701 G4cout <<
"Warning from G4PenelopeIonisation: problem with energy conservation: " <<
702 (energySecondary+energyInFluorescence+kineticEnergy1+localEnergyDeposit+energyInAuger)/
keV <<
703 " keV (final) vs. " <<
704 kineticEnergy0/
keV <<
" keV (initial)" << G4endl;
710 void G4PenelopeIonisationModel::SampleFinalStateElectron(
const G4Material* mat,
723 size_t numberOfOscillators = theTable->size();
731 targetOscillator = numberOfOscillators-1;
734 for (
size_t i=0;i<numberOfOscillators-1;i++)
746 targetOscillator = (
G4int) i;
752 if (verboseLevel > 3)
754 G4cout <<
"SampleFinalStateElectron: sampled oscillator #" << targetOscillator <<
"." <<
G4endl;
755 G4cout <<
"Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/
eV <<
757 G4cout <<
"Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/
eV <<
" eV "
768 G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
770 G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
771 G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
779 if (resEne > cutEnergy && resEne < kineticEnergy)
781 cps = kineticEnergy*rb;
784 if (resEne > 1.0e-6*kineticEnergy)
798 XHDT = XHDT0*invResEne;
817 G4double EE = kineticEnergy + ionEne;
826 XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
827 (1.0-amol)*std::log(rcl*rrl1))/EE;
834 if (XHTOT < 1.e-14*
barn)
836 kineticEnergy1=kineticEnergy;
839 cosThetaSecondary=1.0;
840 targetOscillator = numberOfOscillators-1;
862 rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));
864 rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
867 G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf);
874 kineticEnergy1 = kineticEnergy - deltaE;
875 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
877 energySecondary = deltaE - ionEne;
878 cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*
electron_mass_c2)));
879 if (verboseLevel > 3)
880 G4cout <<
"SampleFinalStateElectron: sampled close collision " <<
G4endl;
887 kineticEnergy1 = kineticEnergy - deltaE;
896 cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
897 if (cosThetaPrimary > 1.)
898 cosThetaPrimary = 1.0;
900 energySecondary = deltaE - ionEne;
901 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
902 if (cosThetaSecondary > 1.0)
903 cosThetaSecondary = 1.0;
904 if (verboseLevel > 3)
905 G4cout <<
"SampleFinalStateElectron: sampled distant longitudinal collision " <<
G4endl;
910 cosThetaPrimary = 1.0;
912 energySecondary = deltaE - ionEne;
913 cosThetaSecondary = 0.5;
914 if (verboseLevel > 3)
915 G4cout <<
"SampleFinalStateElectron: sampled distant transverse collision " <<
G4endl;
923 void G4PenelopeIonisationModel::SampleFinalStatePositron(
const G4Material* mat,
936 size_t numberOfOscillators = theTable->size();
944 targetOscillator = numberOfOscillators-1;
946 for (
size_t i=0;i<numberOfOscillators-1;i++)
951 targetOscillator = (
G4int) i;
956 if (verboseLevel > 3)
958 G4cout <<
"SampleFinalStatePositron: sampled oscillator #" << targetOscillator <<
"." <<
G4endl;
959 G4cout <<
"Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/
eV
961 G4cout <<
"Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/
eV
973 G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0);
975 G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
976 G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12;
981 G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
983 G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
984 G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
993 if (resEne > cutEnergy && resEne < kineticEnergy)
995 cps = kineticEnergy*rb;
998 if (resEne > 1.0e-6*kineticEnergy)
1012 XHDT = XHDT0*invResEne;
1037 XHC = ((1.0/rcl-1.0)+bha1*std::log(rcl)+bha2*rl1
1038 + (bha3/2.0)*(rcl*rcl-1.0)
1039 + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy;
1043 G4double XHTOT = XHC + XHDL + XHDT;
1046 if (XHTOT < 1.e-14*
barn)
1048 kineticEnergy1=kineticEnergy;
1049 cosThetaPrimary=1.0;
1050 energySecondary=0.0;
1051 cosThetaSecondary=1.0;
1052 targetOscillator = numberOfOscillators-1;
1065 G4bool loopAgain =
false;
1070 G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1075 G4double deltaE = rk*kineticEnergy;
1076 kineticEnergy1 = kineticEnergy - deltaE;
1077 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
1079 energySecondary = deltaE - ionEne;
1080 cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*
electron_mass_c2)));
1081 if (verboseLevel > 3)
1082 G4cout <<
"SampleFinalStatePositron: sampled close collision " <<
G4endl;
1089 kineticEnergy1 = kineticEnergy - deltaE;
1097 cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
1098 if (cosThetaPrimary > 1.)
1099 cosThetaPrimary = 1.0;
1101 energySecondary = deltaE - ionEne;
1102 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
1103 if (cosThetaSecondary > 1.0)
1104 cosThetaSecondary = 1.0;
1105 if (verboseLevel > 3)
1106 G4cout <<
"SampleFinalStatePositron: sampled distant longitudinal collision " <<
G4endl;
1111 cosThetaPrimary = 1.0;
1113 energySecondary = deltaE - ionEne;
1114 cosThetaSecondary = 0.5;
1116 if (verboseLevel > 3)
1117 G4cout <<
"SampleFinalStatePositron: sampled distant transverse collision " <<
G4endl;
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4PenelopeOscillatorTable * GetOscillatorTableIonisation(const G4Material *)
G4double LowEnergyLimit() const
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
static constexpr double mm
G4double GetSoftStoppingPower(G4double energy) const
Returns the total stopping power due to soft collisions.
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
void SetVerboseLevel(G4int vl)
Setter for the verbosity level.
G4double HighEnergyLimit() const
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *theParticle, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy=DBL_MAX)
G4bool IsPIXEActive() const
const G4String & GetName() const
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
static G4Electron * Definition()
G4double GetHardCrossSection(G4double energy) const
Returns hard cross section at the given energy.
G4ParticleDefinition * GetDefinition() const
const G4PenelopeCrossSection * GetCrossSectionTableForCouple(const G4ParticleDefinition *, const G4Material *, const G4double cut) const
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *)
G4double BindingEnergy() const
#define G4MUTEX_INITIALIZER
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
void SetHighEnergyLimit(G4double)
G4GLOB_DLL std::ostream G4cout
virtual ~G4PenelopeIonisationModel()
double A(double temperature)
size_t GetTableSize() const
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
G4double GetTotalCrossSection(G4double energy) const
Returns total cross section at the given energy.
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static constexpr double eV
G4ParticleChangeForLoss * fParticleChange
static G4PenelopeOscillatorManager * GetOscillatorManager()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double GetTotNbOfAtomsPerVolume() const
static G4ProductionCutsTable * GetProductionCutsTable()
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
static G4Positron * Positron()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
static G4EmParameters * Instance()
G4double GetDensityCorrection(const G4Material *, const G4double energy) const
Returns the density coeection for the material at the given energy.
static constexpr double GeV
static G4Electron * Electron()
static constexpr double pi
G4VAtomDeexcitation * AtomDeexcitation()
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
void BuildXSTable(const G4Material *, G4double cut, const G4ParticleDefinition *, G4bool isMaster=true)
This can be inkoved only by the master.
G4double GetNormalizedShellCrossSection(size_t shellID, G4double energy) const
Returns the hard cross section for the given shell (normalized to 1)
const G4ParticleDefinition * fParticle
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static constexpr double barn
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
G4double bindingEnergy(G4int A, G4int Z)
G4PenelopeIonisationModel(const G4ParticleDefinition *p=0, const G4String &processName="PenIoni")
static G4AtomicTransitionManager * Instance()
static constexpr double keV
const G4String & PIXEElectronCrossSectionModel()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
static G4Gamma * Definition()
const G4Material * GetMaterial() const