51 slaterEffectiveCharge[0] = 0.;
52 slaterEffectiveCharge[1] = 0.;
53 slaterEffectiveCharge[2] = 0.;
58 lowEnergyLimitForZ1 = 0 * eV;
59 lowEnergyLimitForZ2 = 0 * eV;
60 lowEnergyLimitOfModelForZ1 = 100 * eV;
61 lowEnergyLimitOfModelForZ2 = 1 * keV;
62 killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
63 killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
75 G4cout <<
"Rudd ionisation model is constructed " <<
G4endl;
83 fAtomDeexcitation = 0;
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");
150 tableFile[proton] = fileProton;
152 lowEnergyLimit[proton] = lowEnergyLimitForZ1;
153 highEnergyLimit[proton] = 500. * keV;
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;
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"
307 particleDefinition != instance->
GetIon(
"hydrogen")
309 particleDefinition != instance->
GetIon(
"alpha++")
311 particleDefinition != instance->
GetIon(
"alpha+")
313 particleDefinition != instance->
GetIon(
"helium")
321 || particleDefinition == instance->
GetIon(
"hydrogen")
324 lowLim = lowEnergyLimitOfModelForZ1;
326 if ( particleDefinition == instance->
GetIon(
"alpha++")
327 || particleDefinition == instance->
GetIon(
"alpha+")
328 || particleDefinition == instance->
GetIon(
"helium")
331 lowLim = lowEnergyLimitOfModelForZ2;
351 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
352 pos2 = highEnergyLimit.find(particleName);
354 if (pos2 != highEnergyLimit.end())
356 highLim = pos2->second;
363 if (k < lowLim) k = lowLim;
367 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
368 pos = tableData.find(particleName);
370 if (pos != tableData.end())
380 G4Exception(
"G4DNARuddIonisationModel::CrossSectionPerVolume",
"em0002",
386 if (verboseLevel > 2)
388 G4cout <<
"__________________________________" <<
G4endl;
389 G4cout <<
"G4DNARuddIonisationModel - XS INFO START" <<
G4endl;
391 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/cm/cm <<
G4endl;
392 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) <<
G4endl;
395 G4cout <<
"G4DNARuddIonisationModel - XS INFO END" <<
G4endl;
398 return sigma*waterDensity;
410 if (verboseLevel > 3)
412 G4cout <<
"Calling SampleSecondaries() of G4DNARuddIonisationModel"
426 lowLim = killBelowEnergyForZ1;
433 lowLim = killBelowEnergyForZ2;
450 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
451 pos2 = highEnergyLimit.find(particleName);
453 if (pos2 != highEnergyLimit.end())
455 highLim = pos2->second;
458 if (k >= lowLim && k <= highLim)
469 G4int ionizationShell = RandomSelect(k,particleName);
475 if (k<bindingEnergy)
return;
482 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
490 fvect->push_back(dp);
516 size_t secNumberInit = 0;
517 size_t secNumberFinal = 0;
519 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
522 if(fAtomDeexcitation && ionizationShell == 4)
526 secNumberInit = fvect->size();
528 secNumberFinal = fvect->size();
530 if(secNumberFinal > secNumberInit)
532 for (
size_t i=secNumberInit; i<secNumberFinal; ++i)
535 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
538 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
553 if(bindingEnergy < 0.0)
554 G4Exception(
"G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
606 G4double maximumKineticEnergyTransfer = 0.;
612 || particleDefinition == instance->
GetIon(
"hydrogen"))
614 maximumKineticEnergyTransfer = 4. * (electron_mass_c2 / proton_mass_c2) * k;
617 else if (particleDefinition == instance->
GetIon(
"helium")
618 || particleDefinition == instance->
GetIon(
"alpha+")
619 || particleDefinition == instance->
GetIon(
"alpha++"))
621 maximumKineticEnergyTransfer = 4. * (0.511 / 3728) * k;
631 DifferentialCrossSection(particleDefinition, k, value, shell);
632 if (differentialCrossSection >= crossSectionMaximum)
633 crossSectionMaximum = differentialCrossSection;
640 secElecKinetic =
G4UniformRand()* maximumKineticEnergyTransfer;
641 }
while(
G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
646 return (secElecKinetic);
694 G4int ionizationLevelIndex)
712 const G4int j = ionizationLevelIndex;
727 const G4double Bj[5] = { 12.60 * eV, 14.70 * eV, 18.40 * eV, 32.20 * eV, 540
761 const G4double Gj[5] = { 0.99, 1.11, 1.11, 0.52, 1. };
771 G4double w = wBig / Bj[ionizationLevelIndex];
780 G4bool isProtonOrHydrogen =
false;
784 || particleDefinition == instance->
GetIon(
"hydrogen"))
786 isProtonOrHydrogen =
true;
787 tau = (electron_mass_c2 / proton_mass_c2) * k;
790 else if (particleDefinition == instance->
GetIon(
"helium")
791 || particleDefinition == instance->
GetIon(
"alpha+")
792 || particleDefinition == instance->
GetIon(
"alpha++"))
795 tau = (0.511 / 3728.) * k;
799 * std::pow((Ry / Bj[ionizationLevelIndex]), 2);
801 S = 4. *
pi * Bohr_radius * Bohr_radius *
n
805 G4double v2 = tau / Bj[ionizationLevelIndex];
810 G4double wc = 4. * v2 - 2. * v - (Ry / (4. * Bj[ionizationLevelIndex]));
812 wc = 4. * v2 - 2. * v
815 G4double L1 = (
C1 * std::pow(v, (D1))) / (1. + E1 * std::pow(v, (D1 + 4.)));
817 G4double H1 = (A1 * std::log(1. + v2)) / (v2 + (B1 / v2));
818 G4double H2 = (A2 / v2) + (B2 / (v2 * v2));
821 G4double F2 = (L2 * H2) / (L2 + H2);
824 CorrectionFactor(particleDefinition, k) * Gj[j]
825 * (
S / Bj[ionizationLevelIndex])
827 / (std::pow((1. + w), 3)
828 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
831 sigma = CorrectionFactor(particleDefinition, k) * Gj[j]
834 / (std::pow((1. + w), 3)
835 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
837 if ((particleDefinition == instance->
GetIon(
"hydrogen"))
838 && (ionizationLevelIndex == 4))
843 / (std::pow((1. + w), 3)
844 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
851 if (isProtonOrHydrogen)
856 if (particleDefinition == instance->
GetIon(
"alpha++"))
858 slaterEffectiveCharge[0] = 0.;
859 slaterEffectiveCharge[1] = 0.;
860 slaterEffectiveCharge[2] = 0.;
861 sCoefficient[0] = 0.;
862 sCoefficient[1] = 0.;
863 sCoefficient[2] = 0.;
866 else if (particleDefinition == instance->
GetIon(
"alpha+"))
868 slaterEffectiveCharge[0] = 2.0;
870 slaterEffectiveCharge[1] = 2.0;
871 slaterEffectiveCharge[2] = 2.0;
873 sCoefficient[0] = 0.7;
874 sCoefficient[1] = 0.15;
875 sCoefficient[2] = 0.15;
878 else if (particleDefinition == instance->
GetIon(
"helium"))
880 slaterEffectiveCharge[0] = 1.7;
881 slaterEffectiveCharge[1] = 1.15;
882 slaterEffectiveCharge[2] = 1.15;
883 sCoefficient[0] = 0.5;
884 sCoefficient[1] = 0.25;
885 sCoefficient[2] = 0.25;
894 sigma = Gj[j] * (
S / Bj[ionizationLevelIndex])
896 / (std::pow((1. + w), 3)
897 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
903 / (std::pow((1. + w), 3)
904 * (1. +
G4Exp(alphaConst * (w - wc) / v))));
909 zEff -= (sCoefficient[0]
910 * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.)
912 * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.)
914 * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.));
916 return zEff * zEff * sigma;
932 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
933 G4double value = 1. -
G4Exp(-2 * r) * ((2. * r + 2.) * r + 1.);
948 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
950 -
G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
966 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
969 * ((((2. / 3. * r + 4. / 3.) * r + 2.) * r + 2.) * r + 1.);
984 G4double tElectron = 0.511 / 3728. * t;
987 G4double value = std::sqrt(2. * tElectron / H) / (energyTransferred / H)
988 * (slaterEffectiveChg / shellNumber);
1004 }
else if (particleDefinition == instance->
GetIon(
"hydrogen"))
1006 G4double value = (std::log10(k / eV) - 4.2) / 0.5;
1008 return ((0.6 / (1 +
G4Exp(value))) + 0.9);
1029 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator
pos;
1030 pos = tableData.find(particle);
1032 if (pos != tableData.end())
1048 value += valuesBuffer[i];
1059 if (valuesBuffer[i] > value)
1061 delete[] valuesBuffer;
1064 value -= valuesBuffer[i];
1068 delete[] valuesBuffer;
1073 G4Exception(
"G4DNARuddIonisationModel::RandomSelect",
1076 "Model not applicable to particle type.");
1084G4double G4DNARuddIonisationModel::PartialCrossSection(
const G4Track& track)
1096 std::map<G4String, G4double, std::less<G4String> >::iterator pos1;
1097 pos1 = lowEnergyLimit.find(particleName);
1099 if (pos1 != lowEnergyLimit.end())
1101 lowLim = pos1->second;
1104 std::map<G4String, G4double, std::less<G4String> >::iterator pos2;
1105 pos2 = highEnergyLimit.find(particleName);
1107 if (pos2 != highEnergyLimit.end())
1109 highLim = pos2->second;
1112 if (k >= lowLim && k <= highLim)
1114 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator
pos;
1115 pos = tableData.find(particleName);
1117 if (pos != tableData.end())
1126 G4Exception(
"G4DNARuddIonisationModel::PartialCrossSection",
1129 "Model not applicable to particle type.");
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4GLOB_DLL std::ostream G4cout
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual size_t NumberOfComponents(void) const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
virtual G4bool LoadData(const G4String &argFileName)
static G4DNAGenericIonsManager * Instance(void)
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
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual ~G4DNARuddIonisationModel()
G4DNARuddIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationModel")
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
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
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
const G4Track * GetCurrentTrack() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4int GetLeptonNumber() const
G4double GetPDGCharge() const
const G4String & GetParticleName() 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)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)