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)
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();
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;
653 G4double xEl = sinThetaE * std::cos(phiEl);
654 G4double yEl = sinThetaE * std::sin(phiEl);
657 eDirection.rotateUz(particleDirection0);
660 fvect->push_back(electron);
668 if (localEnergyDeposit < 0)
671 <<
"G4PenelopeIonisationModel::SampleSecondaries - Negative energy deposit"
673 localEnergyDeposit=0.;
679 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
680 G4cout <<
"Energy balance from G4PenelopeIonisation" <<
G4endl;
681 G4cout <<
"Incoming primary energy: " << kineticEnergy0/
keV <<
" keV" <<
G4endl;
682 G4cout <<
"-----------------------------------------------------------" <<
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;
691 localEnergyDeposit+energyInAuger)/
keV <<
693 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
699 localEnergyDeposit+energyInAuger-kineticEnergy0);
700 if (energyDiff > 0.05*
keV)
701 G4cout <<
"Warning from G4PenelopeIonisation: problem with energy conservation: " <<
703 " keV (final) vs. " <<
704 kineticEnergy0/
keV <<
" keV (initial)" << G4endl;
723 size_t numberOfOscillators = theTable->size();
734 for (
size_t i=0;i<numberOfOscillators-1;i++)
761 G4double rb = kineticEnergy + 2.0*electron_mass_c2;
779 if (resEne > cutEnergy && resEne < kineticEnergy)
781 cps = kineticEnergy*rb;
784 if (resEne > 1.0e-6*kineticEnergy)
786 G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
787 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
791 QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
792 QM *= (1.0-QM*0.5/electron_mass_c2);
796 XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
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)
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);
878 cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
880 G4cout <<
"SampleFinalStateElectron: sampled close collision " <<
G4endl;
891 G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
893 - (QS*0.5/electron_mass_c2));
894 G4double QTREV = Q*(Q+2.0*electron_mass_c2);
897 if (cosThetaPrimary > 1.)
898 cosThetaPrimary = 1.0;
901 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
905 G4cout <<
"SampleFinalStateElectron: sampled distant longitudinal collision " <<
G4endl;
915 G4cout <<
"SampleFinalStateElectron: sampled distant transverse collision " <<
G4endl;
936 size_t numberOfOscillators = theTable->size();
946 for (
size_t i=0;i<numberOfOscillators-1;i++)
966 G4double rb = kineticEnergy + 2.0*electron_mass_c2;
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;
993 if (resEne > cutEnergy && resEne < kineticEnergy)
995 cps = kineticEnergy*rb;
998 if (resEne > 1.0e-6*kineticEnergy)
1000 G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
1001 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
1005 QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
1006 QM *= (1.0-QM*0.5/electron_mass_c2);
1010 XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
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)
1065 G4bool loopAgain =
false;
1070 G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1075 G4double deltaE = rk*kineticEnergy;
1080 cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
1082 G4cout <<
"SampleFinalStatePositron: sampled close collision " <<
G4endl;
1092 G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
1094 - (QS*0.5/electron_mass_c2));
1095 G4double QTREV = Q*(Q+2.0*electron_mass_c2);
1098 if (cosThetaPrimary > 1.)
1099 cosThetaPrimary = 1.0;
1102 cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
1106 G4cout <<
"SampleFinalStatePositron: sampled distant longitudinal collision " <<
G4endl;
1117 G4cout <<
"SampleFinalStatePositron: sampled distant transverse collision " <<
G4endl;
G4double cosThetaSecondary
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4VAtomDeexcitation * fAtomDeexcitation
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.
CLHEP::Hep3Vector G4ThreeVector
void SampleFinalStateElectron(const G4Material *, G4double cutEnergy, G4double kineticEnergy)
G4double HighEnergyLimit() const
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *theParticle, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy=DBL_MAX)
G4bool IsPIXEActive() const
G4PenelopeOscillatorManager * oscManager
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
Returns the table of cross sections for the given particle, given material and given cut as a G4Penel...
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *)
G4double BindingEnergy() const
#define G4MUTEX_INITIALIZER
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double fIntrinsicLowEnergyLimit
void SetHighEnergyLimit(G4double)
G4GLOB_DLL std::ostream G4cout
virtual ~G4PenelopeIonisationModel()
double A(double temperature)
size_t GetTableSize() const
const G4ThreeVector & GetMomentumDirection() const
G4PenelopeIonisationXSHandler * theCrossSectionHandler
static const double twopi
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)
void SetParticle(const G4ParticleDefinition *)
G4double GetTotNbOfAtomsPerVolume() const
static G4ProductionCutsTable * GetProductionCutsTable()
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
static G4Positron * Positron()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4double fIntrinsicHighEnergyLimit
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 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.
G4ThreeVector G4ParticleMomentum
G4double bindingEnergy(G4int A, G4int Z)
G4PenelopeIonisationModel(const G4ParticleDefinition *p=0, const G4String &processName="PenIoni")
void SampleFinalStatePositron(const G4Material *, G4double cutEnergy, G4double kineticEnergy)
static G4AtomicTransitionManager * Instance()
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