54 slaterEffectiveCharge[0]=0.;
55 slaterEffectiveCharge[1]=0.;
56 slaterEffectiveCharge[2]=0.;
61 lowEnergyLimitForA[1] = 0 * eV;
62 lowEnergyLimitForA[2] = 0 * eV;
63 lowEnergyLimitForA[3] = 0 * eV;
64 lowEnergyLimitOfModelForA[1] = 100 * eV;
65 lowEnergyLimitOfModelForA[4] = 1 * keV;
66 lowEnergyLimitOfModelForA[5] = 0.5 * MeV;
67 killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1];
68 killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4];
69 killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5];
82 G4cout <<
"Rudd ionisation model is constructed " <<
G4endl;
87 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)
117 G4cout <<
"Calling G4DNARuddIonisationExtendedModel::Initialise()" <<
G4endl;
121 G4String fileProton(
"dna/sigma_ionisation_p_rudd");
122 G4String fileHydrogen(
"dna/sigma_ionisation_h_rudd");
123 G4String fileAlphaPlusPlus(
"dna/sigma_ionisation_alphaplusplus_rudd");
124 G4String fileAlphaPlus(
"dna/sigma_ionisation_alphaplus_rudd");
125 G4String fileHelium(
"dna/sigma_ionisation_he_rudd");
126 G4String fileCarbon(
"dna/sigma_ionisation_c_rudd");
127 G4String fileNitrogen(
"dna/sigma_ionisation_n_rudd");
128 G4String fileOxygen(
"dna/sigma_ionisation_o_rudd");
129 G4String fileIron(
"dna/sigma_ionisation_fe_rudd");
158 tableFile[proton] = fileProton;
159 lowEnergyLimit[proton] = lowEnergyLimitForA[1];
160 highEnergyLimit[proton] = 500. * keV;
168 tableData[proton] = tableProton;
173 tableFile[hydrogen] = fileHydrogen;
175 lowEnergyLimit[hydrogen] = lowEnergyLimitForA[1];
176 highEnergyLimit[hydrogen] = 100. * MeV;
183 tableHydrogen->
LoadData(fileHydrogen);
185 tableData[hydrogen] = tableHydrogen;
190 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
192 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForA[4];
193 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
200 tableAlphaPlusPlus->
LoadData(fileAlphaPlusPlus);
202 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
207 tableFile[alphaPlus] = fileAlphaPlus;
209 lowEnergyLimit[alphaPlus] = lowEnergyLimitForA[4];
210 highEnergyLimit[alphaPlus] = 400. * MeV;
217 tableAlphaPlus->
LoadData(fileAlphaPlus);
218 tableData[alphaPlus] = tableAlphaPlus;
223 tableFile[helium] = fileHelium;
225 lowEnergyLimit[helium] = lowEnergyLimitForA[4];
226 highEnergyLimit[helium] = 400. * MeV;
234 tableData[helium] = tableHelium;
239 tableFile[carbon] = fileCarbon;
241 lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->
GetAtomicMass();
242 highEnergyLimit[carbon] = 1e6* particle->
GetAtomicMass() * MeV;
250 tableData[carbon] = tableCarbon;
255 tableFile[oxygen] = fileOxygen;
257 lowEnergyLimit[oxygen] = lowEnergyLimitForA[5]* particle->
GetAtomicMass();
258 highEnergyLimit[oxygen] = 1e6* particle->
GetAtomicMass()* MeV;
266 tableData[oxygen] = tableOxygen;
271 tableFile[nitrogen] = fileNitrogen;
273 lowEnergyLimit[nitrogen] = lowEnergyLimitForA[5]* particle->
GetAtomicMass();
274 highEnergyLimit[nitrogen] = 1e6* particle->
GetAtomicMass()* MeV;
281 tableNitrogen->
LoadData(fileNitrogen);
282 tableData[nitrogen] = tableNitrogen;
287 tableFile[iron] = fileIron;
289 lowEnergyLimit[iron] = lowEnergyLimitForA[5]* particle->
GetAtomicMass();
298 tableData[iron] = tableIron;
352 G4cout <<
"Rudd ionisation model is initialized " <<
G4endl
367 if (isInitialised) {
return; }
369 isInitialised =
true;
380 if (verboseLevel > 3)
381 G4cout <<
"Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" <<
G4endl;
391 particleDefinition != instance->
GetIon(
"hydrogen")
393 particleDefinition != instance->
GetIon(
"alpha++")
395 particleDefinition != instance->
GetIon(
"alpha+")
397 particleDefinition != instance->
GetIon(
"helium")
399 particleDefinition != instance->
GetIon(
"carbon")
401 particleDefinition != instance->
GetIon(
"nitrogen")
403 particleDefinition != instance->
GetIon(
"oxygen")
405 particleDefinition != instance->
GetIon(
"iron")
413 || particleDefinition == instance->
GetIon(
"hydrogen")
416 lowLim = lowEnergyLimitOfModelForA[1];
418 else if ( particleDefinition == instance->
GetIon(
"alpha++")
419 || particleDefinition == instance->
GetIon(
"alpha+")
420 || particleDefinition == instance->
GetIon(
"helium")
423 lowLim = lowEnergyLimitOfModelForA[4];
425 else lowLim = lowEnergyLimitOfModelForA[5];
433 if(waterDensity!= 0.0)
438 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
439 pos2 = highEnergyLimit.find(particleName);
441 if (pos2 != highEnergyLimit.end())
443 highLim = pos2->second;
451 if (k < lowLim) k = lowLim;
455 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
456 pos = tableData.find(particleName);
458 if (pos != tableData.end())
468 G4Exception(
"G4DNARuddIonisationExtendedModel::CrossSectionPerVolume",
"em0002",
474 if (verboseLevel > 2)
476 G4cout <<
"__________________________________" <<
G4endl;
477 G4cout <<
"°°° G4DNARuddIonisationExtendedModel - XS INFO START" <<
G4endl;
479 G4cout <<
"°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm <<
G4endl;
480 G4cout <<
"°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) <<
G4endl;
482 G4cout <<
"°°° G4DNARuddIonisationExtendedModel - XS INFO END" <<
G4endl;
488 return sigma*waterDensity;
501 if (verboseLevel > 3)
502 G4cout <<
"Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" <<
G4endl;
543 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
544 pos2 = highEnergyLimit.find(particleName);
546 if (pos2 != highEnergyLimit.end())highLim = pos2->second;
548 if (k >= lowLim && k < highLim)
559 G4int ionizationShell = RandomSelect(k,particleName);
565 G4int secNumberInit = 0;
566 G4int secNumberFinal = 0;
570 if(fAtomDeexcitation) {
574 if (ionizationShell <5 && ionizationShell >1)
578 else if (ionizationShell <2)
595 secNumberInit = fvect->size();
597 secNumberFinal = fvect->size();
601 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
605 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi, ionizationShell);
607 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
608 G4double dirX = sinTheta*std::cos(phi);
609 G4double dirY = sinTheta*std::sin(phi);
612 deltaDirection.
rotateUz(primaryDirection);
633 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
635 for (
G4int j=secNumberInit; j < secNumberFinal; j++) {
637 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
645 fvect->push_back(dp);
678 proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell);
682 for(
G4double en=0.; en<20.; en+=1.)
if(RejectionFunction(particleDefinition, k, en, shell) > max1)
683 max1=RejectionFunction(particleDefinition, k, en, shell);
687 value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
689 }
while(random1 > value_sampling);
691 return(proposed_energy);
697void G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronDirection(
G4ParticleDefinition* particleDefinition,
705 G4double maximumEnergyTransfer = 0.;
731 if( (k/MeV)/(particleDefinition->
GetPDGMass()/MeV) <= 0.1 )
733 maximumEnergyTransfer= 4.* (electron_mass_c2 / particleDefinition->
GetPDGMass()) * k;
739 G4double en_per_nucleon = k/approx_nuc_number;
740 G4double beta2 = 1. - 1./pow( (1.+(en_per_nucleon/electron_mass_c2)*(electron_mass_c2/proton_mass_c2)), 2.);
742 maximumEnergyTransfer = 2.*electron_mass_c2*(gamma*gamma-1.)/(1.+2.*gamma*(electron_mass_c2/particleDefinition->
GetPDGMass())+pow(electron_mass_c2/particleDefinition->
GetPDGMass(), 2.) );
746 maxSecKinetic = maximumEnergyTransfer-waterStructure.
IonisationEnergy(shell);
750 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
760 G4int ionizationLevelIndex)
762 const G4int j=ionizationLevelIndex;
765 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
770 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
782 Bj_energy = Bj[ionizationLevelIndex];
785 G4double energyTransfer = proposed_ws + Bj_energy;
786 proposed_ws/=Bj_energy;
791 tau = (electron_mass_c2 / particleDefinition->
GetPDGMass()) * k;
797 if((tau/MeV)<5.447761194e-2)
799 v2 = tau / Bj_energy;
800 beta2 = 2.*tau / electron_mass_c2;
805 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
806 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
810 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
811 G4double rejection_term = 1.+std::exp(alphaConst*(proposed_ws - wc) / v);
812 rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
818 || particleDefinition == instance->
GetIon(
"hydrogen")
821 return(rejection_term);
827 G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.));
828 G4double Zeffion = Z*(1.-std::exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
829 rejection_term*=Zeffion*Zeffion;
832 else if (particleDefinition == instance->
GetIon(
"alpha++") )
835 slaterEffectiveCharge[0]=0.;
836 slaterEffectiveCharge[1]=0.;
837 slaterEffectiveCharge[2]=0.;
843 else if (particleDefinition == instance->
GetIon(
"alpha+") )
846 slaterEffectiveCharge[0]=2.0;
848 slaterEffectiveCharge[1]=2.0;
849 slaterEffectiveCharge[2]=2.0;
852 sCoefficient[1]=0.15;
853 sCoefficient[2]=0.15;
856 else if (particleDefinition == instance->
GetIon(
"helium") )
859 slaterEffectiveCharge[0]=1.7;
860 slaterEffectiveCharge[1]=1.15;
861 slaterEffectiveCharge[2]=1.15;
863 sCoefficient[1]=0.25;
864 sCoefficient[2]=0.25;
876 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
877 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
878 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
880 rejection_term*= zEff * zEff;
883 return (rejection_term);
892 G4int ionizationLevelIndex)
895 const G4int j=ionizationLevelIndex;
904 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
940 Bj_energy = Bj[ionizationLevelIndex];
945 tau = (electron_mass_c2 / particle->
GetPDGMass()) * k;
951 if((tau/MeV)<5.447761194e-2)
953 v2 = tau / Bj_energy;
954 beta2 = 2.*tau / electron_mass_c2;
959 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
960 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
965 G4double L1 = (
C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
967 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
968 G4double H2 = (A2/v2) + (B2/(v2*v2));
976 if( (k/MeV)/(particle->
GetPDGMass()/MeV) <= 0.1 )
978 maximumEnergy = 4.* (electron_mass_c2 / particle->
GetPDGMass()) * k;
984 maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
985 (1.+2.*gamma*(electron_mass_c2/particle->
GetPDGMass())+pow(electron_mass_c2/particle->
GetPDGMass(), 2.) );
992 G4double wmax = maximumEnergy/Bj_energy;
993 G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
996 G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
997 proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
999 proposed_ws/= ( F1*c + F2*c - 2.*randVal );
1000 proposed_ws*=Bj_energy;
1002 return(proposed_ws);
1016 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1017 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
1032 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1033 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
1049 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1050 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
1065 G4double tElectron = 0.511/3728. * t;
1068 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
1081 if (particleDefinition == instance->
GetIon(
"hydrogen") && shell < 4)
1083 G4double value = (std::log10(k/eV)-4.2)/0.5;
1085 return((0.6/(1+std::exp(value))) + 0.9);
1102 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1103 pos = tableData.find(particle);
1105 if (pos != tableData.end())
1122 value += valuesBuffer[i];
1133 if (valuesBuffer[i] > value)
1135 delete[] valuesBuffer;
1138 value -= valuesBuffer[i];
1141 if (valuesBuffer)
delete[] valuesBuffer;
1147 G4Exception(
"G4DNARuddIonisationExtendedModel::RandomSelect",
"em0002",
1156G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(
const G4Track& track )
1168 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1169 pos1 = lowEnergyLimit.find(particleName);
1171 if (pos1 != lowEnergyLimit.end())
1173 lowLim = pos1->second;
1176 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1177 pos2 = highEnergyLimit.find(particleName);
1179 if (pos2 != highEnergyLimit.end())
1181 highLim = pos2->second;
1184 if (k >= lowLim && k <= highLim)
1186 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1187 pos = tableData.find(particleName);
1189 if (pos != tableData.end())
1199 G4Exception(
"G4DNARuddIonisationExtendedModel::PartialCrossSection",
"em0002",
G4DLLIMPORT std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
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)
static G4DNAMolecularMaterial * Instance()
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual ~G4DNARuddIonisationExtendedModel()
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4DNARuddIonisationExtendedModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationExtendedModel")
G4double IonisationEnergy(G4int level)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
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 GetAtomicNumber() const
G4double GetPDGMass() const
G4int GetAtomicMass() const
G4int GetLeptonNumber() const
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
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
void SetHighEnergyLimit(G4double)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void SetLowEnergyLimit(G4double)
void SetDeexcitationFlag(G4bool val)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)