72 :
G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false),
73 fAtomDeexcitation(0),kineticEnergy1(0.*
eV),
74 cosThetaPrimary(1.0),energySecondary(0.*
eV),
75 cosThetaSecondary(0.0),targetOscillator(-1),
76 theCrossSectionHandler(0),fLocalTable(false)
78 fIntrinsicLowEnergyLimit = 100.0*
eV;
79 fIntrinsicHighEnergyLimit = 100.0*
GeV;
109 if (theCrossSectionHandler)
110 delete theCrossSectionHandler;
119 if (verboseLevel > 3)
120 G4cout <<
"Calling G4PenelopeIonisationModel::Initialise()" <<
G4endl;
124 if (!fAtomDeexcitation)
127 G4cout <<
"WARNING from G4PenelopeIonisationModel " <<
G4endl;
128 G4cout <<
"Atomic de-excitation module is not instantiated, so there will not be ";
130 G4cout <<
"Please make sure this is intended" <<
G4endl;
132 SetParticle(particle);
141 nBins =
std::max(nBins,(
size_t)100);
144 if (theCrossSectionHandler)
146 delete theCrossSectionHandler;
147 theCrossSectionHandler = 0;
160 theCrossSectionHandler->
BuildXSTable(theMat,theCuts.at(i),particle,
165 if (verboseLevel > 2) {
166 G4cout <<
"Penelope Ionisation model v2008 is initialized " << G4endl
178 isInitialised =
true;
188 if (verboseLevel > 3)
189 G4cout <<
"Calling G4PenelopeIonisationModel::InitialiseLocal()" <<
G4endl;
202 theCrossSectionHandler = theModel->theCrossSectionHandler;
205 nBins = theModel->nBins;
208 verboseLevel = theModel->verboseLevel;
240 if (verboseLevel > 3)
241 G4cout <<
"Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" <<
G4endl;
250 if (!theCrossSectionHandler)
266 if (verboseLevel > 0)
270 ed <<
"Unable to retrieve the cross section table for " << theParticle->
GetParticleName() <<
271 " in " << material->
GetName() <<
", cut = " << cutEnergy/
keV <<
" keV " <<
G4endl;
272 ed <<
"This can happen only in Unit Tests or via G4EmCalculator" <<
G4endl;
273 G4Exception(
"G4PenelopeIonisationModel::CrossSectionPerVolume()",
277 G4AutoLock lock(&PenelopeIonisationModelMutex);
278 theCrossSectionHandler->
BuildXSTable(material,cutEnergy,theParticle);
294 if (verboseLevel > 3)
295 G4cout <<
"Material " << material->
GetName() <<
" has " << atPerMol <<
296 "atoms per molecule" <<
G4endl;
300 moleculeDensity = atomDensity/atPerMol;
302 G4double crossPerVolume = crossPerMolecule*moleculeDensity;
304 if (verboseLevel > 2)
307 G4cout <<
"Mean free path for delta emission > " << cutEnergy/
keV <<
" keV at " <<
308 energy/
keV <<
" keV = " << (1./crossPerVolume)/
mm <<
" mm" << G4endl;
311 G4cout <<
"Total free path for ionisation (no threshold) at " <<
312 energy/
keV <<
" keV = " << (1./totalCross)/
mm <<
" mm" << G4endl;
314 return crossPerVolume;
329 G4cout <<
"*** G4PenelopeIonisationModel -- WARNING ***" <<
G4endl;
330 G4cout <<
"Penelope Ionisation model v2008 does not calculate cross section _per atom_ " <<
G4endl;
331 G4cout <<
"so the result is always zero. For physics values, please invoke " <<
G4endl;
332 G4cout <<
"GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" <<
G4endl;
358 if (verboseLevel > 3)
359 G4cout <<
"Calling ComputeDEDX() of G4PenelopeIonisationModel" <<
G4endl;
363 if (!theCrossSectionHandler)
379 if (verboseLevel > 0)
383 ed <<
"Unable to retrieve the cross section table for " << theParticle->
GetParticleName() <<
384 " in " << material->
GetName() <<
", cut = " << cutEnergy/
keV <<
" keV " <<
G4endl;
385 ed <<
"This can happen only in Unit Tests or via G4EmCalculator" <<
G4endl;
386 G4Exception(
"G4PenelopeIonisationModel::ComputeDEDXPerVolume()",
390 G4AutoLock lock(&PenelopeIonisationModelMutex);
391 theCrossSectionHandler->
BuildXSTable(material,cutEnergy,theParticle);
410 moleculeDensity = atomDensity/atPerMol;
412 G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
414 if (verboseLevel > 2)
417 G4cout <<
"Stopping power < " << cutEnergy/
keV <<
" keV at " <<
418 kineticEnergy/
keV <<
" keV = " <<
419 sPowerPerVolume/(
keV/
mm) <<
" keV/mm" << G4endl;
421 return sPowerPerVolume;
429 return fIntrinsicLowEnergyLimit;
471 if (verboseLevel > 3)
472 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeIonisationModel" <<
G4endl;
477 if (kineticEnergy0 <= fIntrinsicLowEnergyLimit)
491 kineticEnergy1=kineticEnergy0;
494 cosThetaSecondary=1.0;
495 targetOscillator = -1;
498 SampleFinalStateElectron(material,cutE,kineticEnergy0);
500 SampleFinalStatePositron(material,cutE,kineticEnergy0);
505 G4Exception(
"G4PenelopeIonisationModel::SamplingSecondaries()",
509 if (energySecondary == 0)
return;
511 if (verboseLevel > 3)
513 G4cout <<
"G4PenelopeIonisationModel::SamplingSecondaries() for " <<
515 G4cout <<
"Final eKin = " << kineticEnergy1 <<
" keV" <<
G4endl;
516 G4cout <<
"Final cosTheta = " << cosThetaPrimary <<
G4endl;
517 G4cout <<
"Delta-ray eKin = " << energySecondary <<
" keV" <<
G4endl;
518 G4cout <<
"Delta-ray cosTheta = " << cosThetaSecondary <<
G4endl;
523 G4double sint = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
525 G4double dirx = sint * std::cos(phiPrimary);
526 G4double diry = sint * std::sin(phiPrimary);
530 electronDirection1.
rotateUz(particleDirection0);
532 if (kineticEnergy1 > 0)
542 G4double ionEnergyInPenelopeDatabase =
543 (*theTable)[targetOscillator]->GetIonisationEnergy();
547 G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
548 G4int Z = (
G4int) (*theTable)[targetOscillator]->GetParentZ();
557 if (Z > 0 && shFlag<30)
559 shell = transitionManager->
Shell(Z,shFlag-1);
567 energySecondary += ionEnergyInPenelopeDatabase-
bindingEnergy;
574 if (energySecondary < 0)
580 localEnergyDeposit += energySecondary;
581 energySecondary = 0.0;
586 if (fAtomDeexcitation && shell)
591 size_t nBefore = fvect->size();
593 size_t nAfter = fvect->size();
595 if (nAfter > nBefore)
597 for (
size_t j=nBefore;j<nAfter;j++)
599 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
600 localEnergyDeposit -= itsEnergy;
602 energyInFluorescence += itsEnergy;
604 energyInAuger += itsEnergy;
611 if (energySecondary > cutE)
614 G4double sinThetaE = std::sqrt(1.-cosThetaSecondary*cosThetaSecondary);
616 G4double xEl = sinThetaE * std::cos(phiEl);
617 G4double yEl = sinThetaE * std::sin(phiEl);
620 eDirection.
rotateUz(particleDirection0);
622 eDirection,energySecondary) ;
623 fvect->push_back(electron);
627 localEnergyDeposit += energySecondary;
631 if (localEnergyDeposit < 0)
634 <<
"G4PenelopeIonisationModel::SampleSecondaries - Negative energy deposit"
636 localEnergyDeposit=0.;
640 if (verboseLevel > 1)
642 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
643 G4cout <<
"Energy balance from G4PenelopeIonisation" <<
G4endl;
644 G4cout <<
"Incoming primary energy: " << kineticEnergy0/
keV <<
" keV" <<
G4endl;
645 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
646 G4cout <<
"Outgoing primary energy: " << kineticEnergy1/
keV <<
" keV" <<
G4endl;
648 if (energyInFluorescence)
649 G4cout <<
"Fluorescence x-rays: " << energyInFluorescence/
keV <<
" keV" <<
G4endl;
651 G4cout <<
"Auger electrons: " << energyInAuger/
keV <<
" keV" <<
G4endl;
652 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
653 G4cout <<
"Total final state: " << (energySecondary+energyInFluorescence+kineticEnergy1+
654 localEnergyDeposit+energyInAuger)/
keV <<
656 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
659 if (verboseLevel > 0)
661 G4double energyDiff = std::fabs(energySecondary+energyInFluorescence+kineticEnergy1+
662 localEnergyDeposit+energyInAuger-kineticEnergy0);
663 if (energyDiff > 0.05*
keV)
664 G4cout <<
"Warning from G4PenelopeIonisation: problem with energy conservation: " <<
665 (energySecondary+energyInFluorescence+kineticEnergy1+localEnergyDeposit+energyInAuger)/
keV <<
666 " keV (final) vs. " <<
667 kineticEnergy0/
keV <<
" keV (initial)" << G4endl;
673 void G4PenelopeIonisationModel::SampleFinalStateElectron(
const G4Material* mat,
686 size_t numberOfOscillators = theTable->size();
694 targetOscillator = numberOfOscillators-1;
697 for (
size_t i=0;i<numberOfOscillators-1;i++)
709 targetOscillator = (
G4int) i;
715 if (verboseLevel > 3)
717 G4cout <<
"SampleFinalStateElectron: sampled oscillator #" << targetOscillator <<
"." <<
G4endl;
718 G4cout <<
"Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/
eV <<
720 G4cout <<
"Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/
eV <<
" eV "
731 G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
733 G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
734 G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
742 if (resEne > cutEnergy && resEne < kineticEnergy)
744 cps = kineticEnergy*rb;
747 if (resEne > 1.0e-6*kineticEnergy)
761 XHDT = XHDT0*invResEne;
780 G4double EE = kineticEnergy + ionEne;
789 XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
790 (1.0-amol)*std::log(rcl*rrl1))/EE;
797 if (XHTOT < 1.e-14*
barn)
799 kineticEnergy1=kineticEnergy;
802 cosThetaSecondary=1.0;
803 targetOscillator = numberOfOscillators-1;
825 rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));
827 rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
830 G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf);
837 kineticEnergy1 = kineticEnergy - deltaE;
838 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
840 energySecondary = deltaE - ionEne;
841 cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*
electron_mass_c2)));
842 if (verboseLevel > 3)
843 G4cout <<
"SampleFinalStateElectron: sampled close collision " <<
G4endl;
850 kineticEnergy1 = kineticEnergy - deltaE;
859 cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
860 if (cosThetaPrimary > 1.)
861 cosThetaPrimary = 1.0;
863 energySecondary = deltaE - ionEne;
864 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
865 if (cosThetaSecondary > 1.0)
866 cosThetaSecondary = 1.0;
867 if (verboseLevel > 3)
868 G4cout <<
"SampleFinalStateElectron: sampled distant longitudinal collision " <<
G4endl;
873 cosThetaPrimary = 1.0;
875 energySecondary = deltaE - ionEne;
876 cosThetaSecondary = 0.5;
877 if (verboseLevel > 3)
878 G4cout <<
"SampleFinalStateElectron: sampled distant transverse collision " <<
G4endl;
886 void G4PenelopeIonisationModel::SampleFinalStatePositron(
const G4Material* mat,
899 size_t numberOfOscillators = theTable->size();
907 targetOscillator = numberOfOscillators-1;
909 for (
size_t i=0;i<numberOfOscillators-1;i++)
914 targetOscillator = (
G4int) i;
919 if (verboseLevel > 3)
921 G4cout <<
"SampleFinalStatePositron: sampled oscillator #" << targetOscillator <<
"." <<
G4endl;
922 G4cout <<
"Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/
eV
924 G4cout <<
"Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/
eV
936 G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0);
938 G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
939 G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12;
944 G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
946 G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
947 G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
956 if (resEne > cutEnergy && resEne < kineticEnergy)
958 cps = kineticEnergy*rb;
961 if (resEne > 1.0e-6*kineticEnergy)
975 XHDT = XHDT0*invResEne;
1000 XHC = ((1.0/rcl-1.0)+bha1*std::log(rcl)+bha2*rl1
1001 + (bha3/2.0)*(rcl*rcl-1.0)
1002 + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy;
1006 G4double XHTOT = XHC + XHDL + XHDT;
1009 if (XHTOT < 1.e-14*
barn)
1011 kineticEnergy1=kineticEnergy;
1012 cosThetaPrimary=1.0;
1013 energySecondary=0.0;
1014 cosThetaSecondary=1.0;
1015 targetOscillator = numberOfOscillators-1;
1028 G4bool loopAgain =
false;
1033 G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1038 G4double deltaE = rk*kineticEnergy;
1039 kineticEnergy1 = kineticEnergy - deltaE;
1040 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
1042 energySecondary = deltaE - ionEne;
1043 cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*
electron_mass_c2)));
1044 if (verboseLevel > 3)
1045 G4cout <<
"SampleFinalStatePositron: sampled close collision " <<
G4endl;
1052 kineticEnergy1 = kineticEnergy - deltaE;
1060 cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
1061 if (cosThetaPrimary > 1.)
1062 cosThetaPrimary = 1.0;
1064 energySecondary = deltaE - ionEne;
1065 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
1066 if (cosThetaSecondary > 1.0)
1067 cosThetaSecondary = 1.0;
1068 if (verboseLevel > 3)
1069 G4cout <<
"SampleFinalStatePositron: sampled distant longitudinal collision " <<
G4endl;
1074 cosThetaPrimary = 1.0;
1076 energySecondary = deltaE - ionEne;
1077 cosThetaSecondary = 0.5;
1079 if (verboseLevel > 3)
1080 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()
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)
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)
void SetHighEnergyLimit(G4double)
double precision function energy(A, Z)
G4GLOB_DLL std::ostream G4cout
virtual ~G4PenelopeIonisationModel()
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)
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
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
G4double GetDensityCorrection(const G4Material *, const G4double energy) const
Returns the density coeection for the material at the given energy.
static G4Electron * Electron()
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 &)
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()
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