53 slaterEffectiveCharge[0] = 0.;
54 slaterEffectiveCharge[1] = 0.;
55 slaterEffectiveCharge[2] = 0.;
60 lowEnergyLimitForZ1 = 0 * eV;
61 lowEnergyLimitForZ2 = 0 * eV;
62 lowEnergyLimitOfModelForZ1 = 100 * eV;
63 lowEnergyLimitOfModelForZ2 = 1 * keV;
64 killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
65 killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
77 G4cout <<
"Rudd ionisation model is constructed " <<
G4endl;
97 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
98 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
116 if (verboseLevel > 3)
118 G4cout <<
"Calling G4DNARuddIonisationModel::Initialise()" <<
G4endl;
123 G4String fileProton(
"dna/sigma_ionisation_p_rudd");
124 G4String fileHydrogen(
"dna/sigma_ionisation_h_rudd");
125 G4String fileAlphaPlusPlus(
"dna/sigma_ionisation_alphaplusplus_rudd");
126 G4String fileAlphaPlus(
"dna/sigma_ionisation_alphaplus_rudd");
127 G4String fileHelium(
"dna/sigma_ionisation_he_rudd");
132 hydrogenDef = instance->
GetIon(
"hydrogen");
134 alphaPlusDef = instance->
GetIon(
"alpha+");
135 heliumDef = instance->
GetIon(
"helium");
150 tableFile[proton] = fileProton;
152 lowEnergyLimit[proton] = lowEnergyLimitForZ1;
153 highEnergyLimit[proton] = 500. * keV;
160 tableProton->LoadData(fileProton);
161 tableData[proton] = tableProton;
166 tableFile[hydrogen] = fileHydrogen;
168 lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
169 highEnergyLimit[hydrogen] = 100. * MeV;
176 tableHydrogen->LoadData(fileHydrogen);
178 tableData[hydrogen] = tableHydrogen;
183 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
185 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
186 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
193 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
195 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
200 tableFile[alphaPlus] = fileAlphaPlus;
202 lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
203 highEnergyLimit[alphaPlus] = 400. * MeV;
210 tableAlphaPlus->LoadData(fileAlphaPlus);
211 tableData[alphaPlus] = tableAlphaPlus;
216 tableFile[helium] = fileHelium;
218 lowEnergyLimit[helium] = lowEnergyLimitForZ2;
219 highEnergyLimit[helium] = 400. * MeV;
226 tableHelium->LoadData(fileHelium);
227 tableData[helium] = tableHelium;
231 if (particle==protonDef)
237 if (particle==hydrogenDef)
243 if (particle==heliumDef)
249 if (particle==alphaPlusDef)
255 if (particle==alphaPlusPlusDef)
263 G4cout <<
"Rudd ionisation model is initialized " <<
G4endl
281 isInitialised =
true;
293 if (verboseLevel > 3)
295 G4cout <<
"Calling CrossSectionPerVolume() of G4DNARuddIonisationModel"
302 particleDefinition != protonDef
304 particleDefinition != hydrogenDef
306 particleDefinition != alphaPlusPlusDef
308 particleDefinition != alphaPlusDef
310 particleDefinition != heliumDef
317 if ( particleDefinition == protonDef
318 || particleDefinition == hydrogenDef
321 lowLim = lowEnergyLimitOfModelForZ1;
323 if ( particleDefinition == alphaPlusPlusDef
324 || particleDefinition == alphaPlusDef
325 || particleDefinition == heliumDef
328 lowLim = lowEnergyLimitOfModelForZ2;
348 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
349 pos2 = highEnergyLimit.find(particleName);
351 if (pos2 != highEnergyLimit.end())
353 highLim = pos2->second;
360 if (k < lowLim) k = lowLim;
364 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
365 pos = tableData.find(particleName);
367 if (pos != tableData.end())
370 if (table !=
nullptr)
377 G4Exception(
"G4DNARuddIonisationModel::CrossSectionPerVolume",
"em0002",
383 if (verboseLevel > 2)
385 G4cout <<
"__________________________________" <<
G4endl;
386 G4cout <<
"G4DNARuddIonisationModel - XS INFO START" <<
G4endl;
388 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/cm/cm <<
G4endl;
389 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) <<
G4endl;
392 G4cout <<
"G4DNARuddIonisationModel - XS INFO END" <<
G4endl;
395 return sigma*waterDensity;
407 if (verboseLevel > 3)
409 G4cout <<
"Calling SampleSecondaries() of G4DNARuddIonisationModel"
420 lowLim = killBelowEnergyForZ1;
427 lowLim = killBelowEnergyForZ2;
444 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
445 pos2 = highEnergyLimit.find(particleName);
447 if (pos2 != highEnergyLimit.end())
449 highLim = pos2->second;
452 if (k >= lowLim && k <= highLim)
463 G4int ionizationShell = RandomSelect(k,particleName);
469 if (k<bindingEnergy)
return;
476 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
484 fvect->push_back(dp);
510 size_t secNumberInit = 0;
511 size_t secNumberFinal = 0;
513 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
516 if((fAtomDeexcitation !=
nullptr) && ionizationShell == 4)
520 secNumberInit = fvect->size();
522 secNumberFinal = fvect->size();
524 if(secNumberFinal > secNumberInit)
526 for (
size_t i=secNumberInit; i<secNumberFinal; ++i)
529 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
532 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
547 if(bindingEnergy < 0.0)
548 G4Exception(
"G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
600 G4double maximumKineticEnergyTransfer = 0.;
602 if (particleDefinition == protonDef
603 || particleDefinition == hydrogenDef)
605 maximumKineticEnergyTransfer = 4. * (electron_mass_c2 / proton_mass_c2) * k;
608 else if (particleDefinition == heliumDef
609 || particleDefinition == alphaPlusDef
610 || particleDefinition == alphaPlusPlusDef)
612 maximumKineticEnergyTransfer = 4. * (0.511 / 3728) * k;
622 DifferentialCrossSection(particleDefinition, k, value, shell);
623 if (differentialCrossSection >= crossSectionMaximum)
624 crossSectionMaximum = differentialCrossSection;
631 secElecKinetic =
G4UniformRand()* maximumKineticEnergyTransfer;
632 }
while(
G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
637 return (secElecKinetic);
685 G4int ionizationLevelIndex)
703 const G4int j = ionizationLevelIndex;
718 const G4double Bj[5] = { 12.60 * eV, 14.70 * eV, 18.40 * eV, 32.20 * eV, 540
752 const G4double Gj[5] = { 0.99, 1.11, 1.11, 0.52, 1. };
759 G4double w = wBig / Bj[ionizationLevelIndex];
768 G4bool isProtonOrHydrogen =
false;
771 if (particleDefinition == protonDef
772 || particleDefinition == hydrogenDef)
774 isProtonOrHydrogen =
true;
775 tau = (electron_mass_c2 / proton_mass_c2) * k;
778 else if (particleDefinition == heliumDef
779 || particleDefinition == alphaPlusDef
780 || particleDefinition == alphaPlusPlusDef)
783 tau = (0.511 / 3728.) * k;
786 G4double S = 4. * pi * Bohr_radius * Bohr_radius *
n
787 * gpow->
powN((Ry / Bj[ionizationLevelIndex]), 2);
789 S = 4. * pi * Bohr_radius * Bohr_radius *
n
793 G4double v2 = tau / Bj[ionizationLevelIndex];
798 G4double wc = 4. * v2 - 2. * v - (Ry / (4. * Bj[ionizationLevelIndex]));
800 wc = 4. * v2 - 2. * v
806 G4double H2 = (A2 / v2) + (B2 / (v2 * v2));
809 G4double F2 = (L2 * H2) / (L2 + H2);
812 CorrectionFactor(particleDefinition, k) * Gj[j]
813 * (
S / Bj[ionizationLevelIndex])
815 / (gpow->
powN((1. + w), 3)
816 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
819 sigma = CorrectionFactor(particleDefinition, k) * Gj[j]
822 / (gpow->
powN((1. + w), 3)
823 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
825 if ((particleDefinition == hydrogenDef)
826 && (ionizationLevelIndex == 4))
831 / (gpow->
powN((1. + w), 3)
832 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
839 if (isProtonOrHydrogen)
844 if (particleDefinition == alphaPlusPlusDef)
846 slaterEffectiveCharge[0] = 0.;
847 slaterEffectiveCharge[1] = 0.;
848 slaterEffectiveCharge[2] = 0.;
849 sCoefficient[0] = 0.;
850 sCoefficient[1] = 0.;
851 sCoefficient[2] = 0.;
854 else if (particleDefinition == alphaPlusDef)
856 slaterEffectiveCharge[0] = 2.0;
858 slaterEffectiveCharge[1] = 2.0;
859 slaterEffectiveCharge[2] = 2.0;
861 sCoefficient[0] = 0.7;
862 sCoefficient[1] = 0.15;
863 sCoefficient[2] = 0.15;
866 else if (particleDefinition == heliumDef)
868 slaterEffectiveCharge[0] = 1.7;
869 slaterEffectiveCharge[1] = 1.15;
870 slaterEffectiveCharge[2] = 1.15;
871 sCoefficient[0] = 0.5;
872 sCoefficient[1] = 0.25;
873 sCoefficient[2] = 0.25;
882 sigma = Gj[j] * (
S / Bj[ionizationLevelIndex])
884 / (gpow->
powN((1. + w), 3)
885 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
891 / (gpow->
powN((1. + w), 3)
892 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
897 zEff -= (sCoefficient[0]
898 * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.)
900 * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.)
902 * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.));
904 return zEff * zEff * sigma;
920 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
921 G4double value = 1. -
G4Exp(-2 * r) * ((2. * r + 2.) * r + 1.);
936 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
938 -
G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
954 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
957 * ((((2. / 3. * r + 4. / 3.) * r + 2.) * r + 2.) * r + 1.);
972 G4double tElectron = 0.511 / 3728. * t;
975 G4double value = std::sqrt(2. * tElectron / H) / (energyTransferred / H)
976 * (slaterEffectiveChg / shellNumber);
991 if (particleDefinition == hydrogenDef)
995 return ((0.6 / (1 +
G4Exp(value))) + 0.9);
1014 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator
pos;
1015 pos = tableData.find(particle);
1017 if (pos != tableData.end())
1021 if (table !=
nullptr)
1033 value += valuesBuffer[i];
1044 if (valuesBuffer[i] > value)
1046 delete[] valuesBuffer;
1049 value -= valuesBuffer[i];
1053 delete[] valuesBuffer;
1058 G4Exception(
"G4DNARuddIonisationModel::RandomSelect",
1061 "Model not applicable to particle type.");
1069G4double G4DNARuddIonisationModel::PartialCrossSection(
const G4Track& track)
1081 std::map<G4String, G4double, std::less<G4String> >::iterator pos1;
1082 pos1 = lowEnergyLimit.find(particleName);
1084 if (pos1 != lowEnergyLimit.end())
1086 lowLim = pos1->second;
1089 std::map<G4String, G4double, std::less<G4String> >::iterator pos2;
1090 pos2 = highEnergyLimit.find(particleName);
1092 if (pos2 != highEnergyLimit.end())
1094 highLim = pos2->second;
1097 if (k >= lowLim && k <= highLim)
1099 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator
pos;
1100 pos = tableData.find(particleName);
1102 if (pos != tableData.end())
1105 if (table !=
nullptr)
1111 G4Exception(
"G4DNARuddIonisationModel::PartialCrossSection",
1114 "Model not applicable to particle type.");
G4double S(G4double temp)
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)
G4GLOB_DLL std::ostream G4cout
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
size_t NumberOfComponents() const override
const G4VEMDataSet * GetComponent(G4int componentId) const override
G4double FindValue(G4double e, G4int componentId=0) const override
static G4DNAGenericIonsManager * Instance()
G4ParticleDefinition * GetIon(const G4String &name)
const std::vector< G4double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
static G4DNAMolecularMaterial * Instance()
G4ParticleChangeForGamma * fParticleChangeForGamma
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax) override
~G4DNARuddIonisationModel() override
G4DNARuddIonisationModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="DNARuddIonisationModel")
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4double IonisationEnergy(G4int level)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
std::size_t GetIndex() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4int GetLeptonNumber() const
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Pow * GetInstance()
G4double logZ(G4int Z) const
G4double powN(G4double x, G4int n) const
G4double powA(G4double A, G4double y) const
static G4Proton * ProtonDefinition()
static G4Proton * Proton()
const G4DynamicParticle * GetDynamicParticle() const
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
void SetHighEnergyLimit(G4double)
G4VEmAngularDistribution * GetAngularDistribution()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void SetLowEnergyLimit(G4double)
void SetDeexcitationFlag(G4bool val)
void SetAngularDistribution(G4VEmAngularDistribution *)
void ProposeTrackStatus(G4TrackStatus status)
const G4Track * GetCurrentTrack() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)