86 G4cout <<
"Rudd ionisation model is constructed " <<
G4endl;
108 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
127 G4cout <<
"Calling G4DNARuddIonisationExtendedModel::Initialise()" <<
G4endl;
131 G4String fileProton(
"dna/sigma_ionisation_p_rudd");
132 G4String fileHydrogen(
"dna/sigma_ionisation_h_rudd");
133 G4String fileAlphaPlusPlus(
"dna/sigma_ionisation_alphaplusplus_rudd");
134 G4String fileAlphaPlus(
"dna/sigma_ionisation_alphaplus_rudd");
135 G4String fileHelium(
"dna/sigma_ionisation_he_rudd");
136 G4String fileLithium(
"dna/sigma_ionisation_li_rudd");
137 G4String fileBeryllium(
"dna/sigma_ionisation_be_rudd");
138 G4String fileBoron(
"dna/sigma_ionisation_b_rudd");
139 G4String fileCarbon(
"dna/sigma_ionisation_c_rudd");
140 G4String fileNitrogen(
"dna/sigma_ionisation_n_rudd");
141 G4String fileOxygen(
"dna/sigma_ionisation_o_rudd");
142 G4String fileSilicon(
"dna/sigma_ionisation_si_rudd");
143 G4String fileIron(
"dna/sigma_ionisation_fe_rudd");
215 tableHydrogen->
LoadData(fileHydrogen);
222 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
232 tableAlphaPlusPlus->
LoadData(fileAlphaPlusPlus);
234 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
249 tableAlphaPlus->
LoadData(fileAlphaPlus);
285 tableLithium->
LoadData(fileLithium);
305 tableBeryllium->
LoadData(fileBeryllium);
385 tableNitrogen->
LoadData(fileNitrogen);
404 tableSilicon->
LoadData(fileSilicon);
435 if (particle==protonDef)
441 if (particle==hydrogenDef)
447 if (particle==heliumDef)
453 if (particle==alphaPlusDef)
459 if (particle==alphaPlusPlusDef)
465 if (particle==lithiumDef)
471 if (particle==berylliumDef)
477 if (particle==boronDef)
483 if (particle==carbonDef)
489 if (particle==nitrogenDef)
495 if (particle==oxygenDef)
501 if (particle==siliconDef)
507 if (particle==ironDef)
517 G4cout <<
"Rudd ionisation model is initialized " << G4endl
550 G4cout <<
"Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" <<
G4endl;
560 particleDefinition != instance->
GetIon(
"hydrogen")
562 particleDefinition != instance->
GetIon(
"alpha++")
564 particleDefinition != instance->
GetIon(
"alpha+")
566 particleDefinition != instance->
GetIon(
"helium")
599 || particleDefinition == instance->
GetIon(
"hydrogen")
604 else if ( particleDefinition == instance->
GetIon(
"alpha++")
605 || particleDefinition == instance->
GetIon(
"alpha+")
606 || particleDefinition == instance->
GetIon(
"helium")
619 if(waterDensity!= 0.0)
624 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
629 highLim = pos2->second;
637 if (k < lowLim) k = lowLim;
641 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
654 G4Exception(
"G4DNARuddIonisationExtendedModel::CrossSectionPerVolume",
"em0002",
662 G4cout <<
"__________________________________" <<
G4endl;
663 G4cout <<
"G4DNARuddIonisationExtendedModel - XS INFO START" <<
G4endl;
665 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
666 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
668 G4cout <<
"G4DNARuddIonisationExtendedModel - XS INFO END" <<
G4endl;
674 return sigma*waterDensity;
693 G4cout <<
"Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" <<
G4endl;
739 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
744 if (k >= lowLim && k < highLim)
761 G4int secNumberInit = 0;
762 G4int secNumberFinal = 0;
767 if (k<bindingEnergy)
return;
775 if (ionizationShell <5 && ionizationShell >1)
779 else if (ionizationShell <2)
795 secNumberInit = fvect->size();
797 secNumberFinal = fvect->size();
840 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
842 for (
G4int j=secNumberInit; j < secNumberFinal; j++) {
844 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
860 fvect->push_back(dp);
902 value_sampling =
RejectionFunction(particleDefinition, k, proposed_energy, shell);
904 }
while(random1 > value_sampling);
906 return(proposed_energy);
957 G4int ionizationLevelIndex)
959 const G4int j=ionizationLevelIndex;
962 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
979 Bj_energy = Bj[ionizationLevelIndex];
982 G4double energyTransfer = proposed_ws + Bj_energy;
983 proposed_ws/=Bj_energy;
988 tau = (electron_mass_c2 / particleDefinition->
GetPDGMass()) * k;
994 if((tau/
MeV)<5.447761194e-2)
996 v2 = tau / Bj_energy;
997 beta2 = 2.*tau / electron_mass_c2;
1002 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
1003 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
1007 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
1008 G4double rejection_term = 1.+std::exp(alphaConst*(proposed_ws - wc) / v);
1009 rejection_term = (1./rejection_term)*
CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
1015 || particleDefinition == instance->
GetIon(
"hydrogen")
1018 return(rejection_term);
1025 G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.));
1026 G4double Zeffion = Z*(1.-std::exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
1027 rejection_term*=Zeffion*Zeffion;
1030 else if (particleDefinition == instance->
GetIon(
"alpha++") )
1041 else if (particleDefinition == instance->
GetIon(
"alpha+") )
1054 else if (particleDefinition == instance->
GetIon(
"helium") )
1079 rejection_term*= zEff * zEff;
1082 return (rejection_term);
1091 G4int ionizationLevelIndex)
1094 const G4int j=ionizationLevelIndex;
1139 Bj_energy = Bj[ionizationLevelIndex];
1144 tau = (electron_mass_c2 / particle->
GetPDGMass()) * k;
1150 if((tau/
MeV)<5.447761194e-2)
1152 v2 = tau / Bj_energy;
1153 beta2 = 2.*tau / electron_mass_c2;
1158 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
1159 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
1164 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
1166 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
1167 G4double H2 = (A2/v2) + (B2/(v2*v2));
1177 maximumEnergy = 4.* (electron_mass_c2 / particle->
GetPDGMass()) * k;
1182 G4double gamma = 1./sqrt(1.-beta2);
1183 maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
1184 (1.+2.*gamma*(electron_mass_c2/particle->
GetPDGMass())+pow(electron_mass_c2/particle->
GetPDGMass(), 2.) );
1191 G4double wmax = maximumEnergy/Bj_energy;
1192 G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
1195 G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
1196 proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
1198 proposed_ws/= ( F1*c + F2*c - 2.*randVal );
1199 proposed_ws*=Bj_energy;
1201 return(proposed_ws);
1215 G4double r =
R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1216 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
1231 G4double r =
R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1232 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
1248 G4double r =
R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1249 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
1264 G4double tElectron = 0.511/3728. * t;
1267 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
1280 if (particleDefinition == instance->
GetIon(
"hydrogen") && shell < 4)
1284 return((0.6/(1+std::exp(value))) + 0.9);
1301 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
1321 value += valuesBuffer[i];
1332 if (valuesBuffer[i] > value)
1334 delete[] valuesBuffer;
1337 value -= valuesBuffer[i];
1340 if (valuesBuffer)
delete[] valuesBuffer;
1346 G4Exception(
"G4DNARuddIonisationExtendedModel::RandomSelect",
"em0002",
1367 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1372 lowLim = pos1->second;
1375 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1380 highLim = pos2->second;
1383 if (k >= lowLim && k <= highLim)
1385 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
1398 G4Exception(
"G4DNARuddIonisationExtendedModel::PartialCrossSection",
"em0002",
G4VAtomDeexcitation * fAtomDeexcitation
G4DNARuddIonisationExtendedModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationExtendedModel")
G4double LowEnergyLimit() const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
G4double slaterEffectiveCharge[3]
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
const G4DynamicParticle * GetDynamicParticle() const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
virtual ~G4DNARuddIonisationExtendedModel()
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *particleDefinition, G4double incomingParticleEnergy, G4int shell)
static G4Proton * ProtonDefinition()
G4VEmAngularDistribution * GetAngularDistribution()
std::map< G4double, G4double > lowEnergyLimitOfModelForA
const std::vector< G4double > * fpWaterDensity
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double S_1s(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
G4int GetAtomicNumber() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4ParticleChangeForGamma * fParticleChangeForGamma
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
G4double RejectionFunction(G4ParticleDefinition *particle, G4double k, G4double proposed_ws, G4int ionizationLevelIndex)
G4double Sum(G4double energy, const G4String &particle)
const G4ThreeVector & GetMomentumDirection() const
G4double S_2p(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
static G4DNAGenericIonsManager * Instance(void)
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4int GetAtomicMass() const
virtual size_t NumberOfComponents(void) const
static G4IonTable * GetIonTable()
G4double IonisationEnergy(G4int level)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double R(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
G4double PartialCrossSection(const G4Track &track)
G4double ProposedSampledEnergy(G4ParticleDefinition *particle, G4double k, G4int ionizationLevelIndex)
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Method used by DNA physics model to create a water molecule.
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
void SetAngularDistribution(G4VEmAngularDistribution *)
const G4double x[NPOINTSGL]
G4DNAWaterIonisationStructure waterStructure
const G4Track * GetCurrentTrack() const
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static MCTruthManager * instance
G4VAtomDeexcitation * AtomDeexcitation()
G4double S_2s(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
std::map< G4double, G4double > killBelowEnergyForA
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
static const double eplus
G4double GetPDGCharge() const
G4double CorrectionFactor(G4ParticleDefinition *particleDefinition, G4double k, G4int shell)
std::map< G4double, G4double > lowEnergyLimitForA
G4ThreeVector G4ParticleMomentum
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4double bindingEnergy(G4int A, G4int Z)
G4int GetLeptonNumber() const
static const G4double pos
G4int RandomSelect(G4double energy, const G4String &particle)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
const G4Material * GetMaterial() const
G4ParticleDefinition * GetIon(const G4String &name)