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;
80 fAtomDeexcitation = 0;
90 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
91 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
110 if (verboseLevel > 3)
111 G4cout <<
"Calling G4DNARuddIonisationModel::Initialise()" <<
G4endl;
115 G4String fileProton(
"dna/sigma_ionisation_p_rudd");
116 G4String fileHydrogen(
"dna/sigma_ionisation_h_rudd");
117 G4String fileAlphaPlusPlus(
"dna/sigma_ionisation_alphaplusplus_rudd");
118 G4String fileAlphaPlus(
"dna/sigma_ionisation_alphaplus_rudd");
119 G4String fileHelium(
"dna/sigma_ionisation_he_rudd");
142 tableFile[proton] = fileProton;
144 lowEnergyLimit[proton] = lowEnergyLimitForZ1;
145 highEnergyLimit[proton] = 500. * keV;
153 tableData[proton] = tableProton;
158 tableFile[hydrogen] = fileHydrogen;
160 lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
161 highEnergyLimit[hydrogen] = 100. * MeV;
168 tableHydrogen->
LoadData(fileHydrogen);
170 tableData[hydrogen] = tableHydrogen;
175 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
177 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
178 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
185 tableAlphaPlusPlus->
LoadData(fileAlphaPlusPlus);
187 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
192 tableFile[alphaPlus] = fileAlphaPlus;
194 lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
195 highEnergyLimit[alphaPlus] = 400. * MeV;
202 tableAlphaPlus->
LoadData(fileAlphaPlus);
203 tableData[alphaPlus] = tableAlphaPlus;
208 tableFile[helium] = fileHelium;
210 lowEnergyLimit[helium] = lowEnergyLimitForZ2;
211 highEnergyLimit[helium] = 400. * MeV;
219 tableData[helium] = tableHelium;
223 if (particle==protonDef)
229 if (particle==hydrogenDef)
235 if (particle==heliumDef)
241 if (particle==alphaPlusDef)
247 if (particle==alphaPlusPlusDef)
255 G4cout <<
"Rudd ionisation model is initialized " <<
G4endl
270 if (isInitialised) {
return; }
272 isInitialised =
true;
284 if (verboseLevel > 3)
285 G4cout <<
"Calling CrossSectionPerVolume() of G4DNARuddIonisationModel" <<
G4endl;
295 particleDefinition != instance->
GetIon(
"hydrogen")
297 particleDefinition != instance->
GetIon(
"alpha++")
299 particleDefinition != instance->
GetIon(
"alpha+")
301 particleDefinition != instance->
GetIon(
"helium")
309 || particleDefinition == instance->
GetIon(
"hydrogen")
312 lowLim = lowEnergyLimitOfModelForZ1;
314 if ( particleDefinition == instance->
GetIon(
"alpha++")
315 || particleDefinition == instance->
GetIon(
"alpha+")
316 || particleDefinition == instance->
GetIon(
"helium")
319 lowLim = lowEnergyLimitOfModelForZ2;
326 if(waterDensity!= 0.0)
342 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
343 pos2 = highEnergyLimit.find(particleName);
345 if (pos2 != highEnergyLimit.end())
347 highLim = pos2->second;
354 if (k < lowLim) k = lowLim;
358 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
359 pos = tableData.find(particleName);
361 if (pos != tableData.end())
371 G4Exception(
"G4DNARuddIonisationModel::CrossSectionPerVolume",
"em0002",
377 if (verboseLevel > 2)
379 G4cout <<
"__________________________________" <<
G4endl;
380 G4cout <<
"°°° G4DNARuddIonisationModel - XS INFO START" <<
G4endl;
382 G4cout <<
"°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm <<
G4endl;
383 G4cout <<
"°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) <<
G4endl;
385 G4cout <<
"°°° G4DNARuddIonisationModel - XS INFO END" <<
G4endl;
391 if (verboseLevel > 2)
393 G4cout <<
"Warning : RuddIonisationModel: WATER DENSITY IS NULL" <<
G4endl;
397 return sigma*waterDensity;
410 if (verboseLevel > 3)
411 G4cout <<
"Calling SampleSecondaries() of G4DNARuddIonisationModel" <<
G4endl;
423 lowLim = killBelowEnergyForZ1;
430 lowLim = killBelowEnergyForZ2;
447 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
448 pos2 = highEnergyLimit.find(particleName);
450 if (pos2 != highEnergyLimit.end())
452 highLim = pos2->second;
455 if (k >= lowLim && k <= highLim)
466 G4int ionizationShell = RandomSelect(k,particleName);
472 G4int secNumberInit = 0;
473 G4int secNumberFinal = 0;
478 if(fAtomDeexcitation) {
482 if (ionizationShell <5 && ionizationShell >1)
486 else if (ionizationShell <2)
503 secNumberInit = fvect->size();
505 secNumberFinal = fvect->size();
509 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
513 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi);
515 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
516 G4double dirX = sinTheta*std::cos(phi);
517 G4double dirY = sinTheta*std::sin(phi);
520 deltaDirection.
rotateUz(primaryDirection);
541 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
543 for (
G4int j=secNumberInit; j < secNumberFinal; j++) {
545 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
559 fvect->push_back(dp);
584 G4double maximumKineticEnergyTransfer = 0.;
590 || particleDefinition == instance->
GetIon(
"hydrogen"))
592 maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
595 else if (particleDefinition == instance->
GetIon(
"helium")
596 || particleDefinition == instance->
GetIon(
"alpha+")
597 || particleDefinition == instance->
GetIon(
"alpha++"))
599 maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
606 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
607 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
615 secElecKinetic =
G4UniformRand() * maximumKineticEnergyTransfer;
616 }
while(
G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
621 return(secElecKinetic);
627void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(
G4ParticleDefinition* particleDefinition,
639 || particleDefinition == instance->
GetIon(
"hydrogen"))
641 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
644 else if (particleDefinition == instance->
GetIon(
"helium")
645 || particleDefinition == instance->
GetIon(
"alpha+")
646 || particleDefinition == instance->
GetIon(
"alpha++"))
648 maxSecKinetic = 4.* (0.511 / 3728) * k;
657 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
668 G4int ionizationLevelIndex)
686 const G4int j=ionizationLevelIndex;
701 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
735 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
741 if (wBig<0)
return 0.;
743 G4double w = wBig / Bj[ionizationLevelIndex];
751 G4bool isProtonOrHydrogen =
false;
755 || particleDefinition == instance->
GetIon(
"hydrogen"))
757 isProtonOrHydrogen =
true;
758 tau = (electron_mass_c2/proton_mass_c2) * k ;
761 else if ( particleDefinition == instance->
GetIon(
"helium")
762 || particleDefinition == instance->
GetIon(
"alpha+")
763 || particleDefinition == instance->
GetIon(
"alpha++"))
766 tau = (0.511/3728.) * k ;
769 G4double S = 4.*
pi * Bohr_radius*Bohr_radius *
n * std::pow((Ry/Bj[ionizationLevelIndex]),2);
770 if (j==4) S = 4.*
pi * Bohr_radius*Bohr_radius *
n * std::pow((Ry/waterStructure.
IonisationEnergy(ionizationLevelIndex)),2);
772 G4double v2 = tau / Bj[ionizationLevelIndex];
776 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj[ionizationLevelIndex]));
777 if (j==4) wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.
IonisationEnergy(ionizationLevelIndex)));
779 G4double L1 = (
C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
781 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
782 G4double H2 = (A2/v2) + (B2/(v2*v2));
787 G4double sigma = CorrectionFactor(particleDefinition, k)
788 * Gj[j] * (S/Bj[ionizationLevelIndex])
789 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
791 if (j==4) sigma = CorrectionFactor(particleDefinition, k)
793 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
795 if ( (particleDefinition == instance->
GetIon(
"hydrogen")) && (ionizationLevelIndex==4))
799 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
804 if(isProtonOrHydrogen)
809 if (particleDefinition == instance->
GetIon(
"alpha++") )
811 slaterEffectiveCharge[0]=0.;
812 slaterEffectiveCharge[1]=0.;
813 slaterEffectiveCharge[2]=0.;
819 else if (particleDefinition == instance->
GetIon(
"alpha+") )
821 slaterEffectiveCharge[0]=2.0;
823 slaterEffectiveCharge[1]=2.0;
824 slaterEffectiveCharge[2]=2.0;
827 sCoefficient[1]=0.15;
828 sCoefficient[2]=0.15;
831 else if (particleDefinition == instance->
GetIon(
"helium") )
833 slaterEffectiveCharge[0]=1.7;
834 slaterEffectiveCharge[1]=1.15;
835 slaterEffectiveCharge[2]=1.15;
837 sCoefficient[1]=0.25;
838 sCoefficient[2]=0.25;
847 sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
849 if (j==4) sigma = Gj[j] * (S/waterStructure.
IonisationEnergy(ionizationLevelIndex))
850 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
854 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
855 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
856 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
858 return zEff * zEff * sigma ;
874 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
875 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
890 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
891 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
907 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
908 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
923 G4double tElectron = 0.511/3728. * t;
926 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
943 if (particleDefinition == instance->
GetIon(
"hydrogen"))
945 G4double value = (std::log10(k/eV)-4.2)/0.5;
947 return((0.6/(1+std::exp(value))) + 0.9);
968 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
969 pos = tableData.find(particle);
971 if (pos != tableData.end())
987 value += valuesBuffer[i];
999 if (valuesBuffer[i] > value)
1001 delete[] valuesBuffer;
1004 value -= valuesBuffer[i];
1007 if (valuesBuffer)
delete[] valuesBuffer;
1013 G4Exception(
"G4DNARuddIonisationModel::RandomSelect",
"em0002",
1022G4double G4DNARuddIonisationModel::PartialCrossSection(
const G4Track& track )
1034 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1035 pos1 = lowEnergyLimit.find(particleName);
1037 if (pos1 != lowEnergyLimit.end())
1039 lowLim = pos1->second;
1042 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1043 pos2 = highEnergyLimit.find(particleName);
1045 if (pos2 != highEnergyLimit.end())
1047 highLim = pos2->second;
1050 if (k >= lowLim && k <= highLim)
1052 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1053 pos = tableData.find(particleName);
1055 if (pos != tableData.end())
1065 G4Exception(
"G4DNARuddIonisationModel::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
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()
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
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)