60 G4cout <<
"Born ionisation model is constructed " <<
G4endl;
65 fAtomDeexcitation = 0;
67 fpMolWaterDensity = 0;
109 if (verboseLevel > 3)
111 G4cout <<
"Calling G4DNABornIonisationModel2::Initialise()" <<
G4endl;
114 if(fParticleDef != 0 && particle != fParticleDef)
117 description <<
"You are trying to initialized G4DNABornIonisationModel2 "
121 description <<
"G4DNABornIonisationModel2 was already initiliased "
123 G4Exception(
"G4DNABornIonisationModel2::Initialise",
"bornIonInit",
127 fParticleDef = particle;
130 char *path = getenv(
"G4LEDATA");
135 std::ostringstream fullFileName;
136 fullFileName << path;
138 if(particleName ==
"e-")
140 fTableFile =
"dna/sigma_ionisation_e_born";
141 fLowEnergyLimit = 11.*
eV;
142 fHighEnergyLimit = 1.*
MeV;
146 fullFileName <<
"/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
150 fullFileName <<
"/dna/sigmadiff_ionisation_e_born.dat";
153 else if(particleName ==
"proton")
155 fTableFile =
"dna/sigma_ionisation_p_born";
156 fLowEnergyLimit = 500. *
keV;
157 fHighEnergyLimit = 100. *
MeV;
161 fullFileName <<
"/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
165 fullFileName <<
"/dna/sigmadiff_ionisation_p_born.dat";
170 G4double scaleFactor = (1.e-22 / 3.343) *
m*
m;
176 std::ifstream diffCrossSection(fullFileName.str().c_str());
178 if (!diffCrossSection)
181 description <<
"Missing data file:" <<
G4endl << fullFileName.str() <<
G4endl;
182 G4Exception(
"G4DNABornIonisationModel2::Initialise",
"em0003",
194 for (
int j=0; j<5; j++)
196 fProbaShellMap[j].clear();
197 fDiffCrossSectionData[j].clear();
198 fNrjTransfData[j].clear();
202 fTdummyVec.push_back(0.);
203 while(!diffCrossSection.eof())
207 diffCrossSection>>tDummy>>eDummy;
208 if (tDummy != fTdummyVec.back()) fTdummyVec.push_back(tDummy);
211 for (
int j=0; j<5; j++)
213 diffCrossSection>> tmp;
215 fDiffCrossSectionData[j][tDummy][eDummy] = tmp;
219 fNrjTransfData[j][tDummy][fDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
220 fProbaShellMap[j][tDummy].push_back(fDiffCrossSectionData[j][tDummy][eDummy]);
224 if (!diffCrossSection.eof() && !fasterCode) fDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
226 if (!fasterCode) fVecm[tDummy].push_back(eDummy);
237 G4cout <<
"Born ionisation model is initialized " <<
G4endl
255 isInitialised =
true;
266 if (verboseLevel > 3)
268 G4cout <<
"Calling CrossSectionPerVolume() of G4DNABornIonisationModel2"
272 if (particleDefinition != fParticleDef)
return 0;
280 if(waterDensity!= 0.0)
282 if (ekin >= fLowEnergyLimit && ekin < fHighEnergyLimit)
290 G4double A = 1.39241700556072800000E-009 ;
291 G4double B = -8.52610412942622630000E-002 ;
292 sigma = sigma *
G4Exp(A*(ekin/
eV)+B);
297 if (verboseLevel > 2)
299 G4cout <<
"__________________________________" <<
G4endl;
300 G4cout <<
"G4DNABornIonisationModel2 - XS INFO START" <<
G4endl;
302 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
303 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
304 G4cout <<
"G4DNABornIonisationModel2 - XS INFO END" <<
G4endl;
308 return sigma*waterDensity;
320 if (verboseLevel > 3)
322 G4cout <<
"Calling SampleSecondaries() of G4DNABornIonisationModel2"
328 if (k >= fLowEnergyLimit && k < fHighEnergyLimit)
332 G4double totalEnergy = k + particleMass;
333 G4double pSquare = k * (totalEnergy + particleMass);
334 G4double totalMomentum = std::sqrt(pSquare);
336 G4int ionizationShell = 0;
338 if (!fasterCode) ionizationShell = RandomSelect(k);
349 ionizationShell = RandomSelect(k);
356 G4int secNumberInit = 0;
357 G4int secNumberFinal = 0;
363 if (k<bindingEnergy)
return;
367 if(fAtomDeexcitation)
371 if (ionizationShell <5 && ionizationShell >1)
375 else if (ionizationShell <2)
390 secNumberInit = fvect->size();
392 secNumberFinal = fvect->size();
397 if (fasterCode ==
false)
399 secondaryKinetic = RandomizeEjectedElectronEnergy(particle->
GetDefinition(),k,ionizationShell);
404 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->
GetDefinition(),k,ionizationShell);
417 G4double finalPx = totalMomentum*primaryDirection.
x() - deltaTotalMomentum*deltaDirection.
x();
418 G4double finalPy = totalMomentum*primaryDirection.
y() - deltaTotalMomentum*deltaDirection.
y();
419 G4double finalPz = totalMomentum*primaryDirection.
z() - deltaTotalMomentum*deltaDirection.
z();
420 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
421 finalPx /= finalMomentum;
422 finalPy /= finalMomentum;
423 finalPz /= finalMomentum;
426 direction.
set(finalPx,finalPy,finalPz);
434 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
436 for (
G4int j=secNumberInit; j < secNumberFinal; j++)
438 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
453 if (secondaryKinetic>0)
456 fvect->push_back(dp);
477 G4double maximumEnergyTransfer = 0.;
479 maximumEnergyTransfer = k;
497 G4double maxEnergy = maximumEnergyTransfer;
498 G4int nEnergySteps = 50;
502 1. / static_cast<G4double>(nEnergySteps - 1)));
503 G4int step(nEnergySteps);
512 if (differentialCrossSection >= crossSectionMaximum)
513 crossSectionMaximum = differentialCrossSection;
518 G4double secondaryElectronKineticEnergy = 0.;
526 return secondaryElectronKineticEnergy;
532 G4double maximumKineticEnergyTransfer = 4.
544 if (differentialCrossSection >= crossSectionMaximum)
545 crossSectionMaximum = differentialCrossSection;
548 G4double secondaryElectronKineticEnergy = 0.;
551 secondaryElectronKineticEnergy =
G4UniformRand()* maximumKineticEnergyTransfer;
556 return secondaryElectronKineticEnergy;
610 G4int ionizationLevelIndex)
614 if (energyTransfer >= waterStructure.
IonisationEnergy(ionizationLevelIndex))
630 std::vector<double>::iterator t2 = std::upper_bound(fTdummyVec.begin(),
634 std::vector<double>::iterator t1 = t2 - 1;
637 if (energyTransfer <= fVecm[(*t1)].back()
638 && energyTransfer <= fVecm[(*t2)].back())
640 std::vector<double>::iterator e12 = std::upper_bound(fVecm[(*t1)].begin(),
643 std::vector<double>::iterator e11 = e12 - 1;
645 std::vector<double>::iterator e22 = std::upper_bound(fVecm[(*t2)].begin(),
648 std::vector<double>::iterator e21 = e22 - 1;
657 xs11 = fDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
658 xs12 = fDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
659 xs21 = fDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
660 xs22 = fDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
664 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
667 sigma = QuadInterpolator(valueE11,
697 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
700 G4double a = (std::log10(xs2) - std::log10(xs1))
701 / (std::log10(e2) - std::log10(e1));
702 G4double b = std::log10(xs2) - a * std::log10(e2);
703 G4double sigma = a * std::log10(e) + b;
704 value = (std::pow(10., sigma));
718 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
722 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
728 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
732 value = (d1 + (d2 -
d1) * (e - e1) / (e2 - e1));
764 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
765 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
800 value += valuesBuffer[i];
811 if (valuesBuffer[i] > value)
813 delete[] valuesBuffer;
816 value -= valuesBuffer[i];
820 delete[] valuesBuffer;
833 G4double secondaryElectronKineticEnergy = 0.;
845 if (secondaryElectronKineticEnergy < 0.)
849 return secondaryElectronKineticEnergy;
856 G4int ionizationLevelIndex,
875 std::vector<double>::iterator k2 = std::upper_bound(fTdummyVec.begin(),
878 std::vector<double>::iterator k1 = k2 - 1;
892 if (random <= fProbaShellMap[ionizationLevelIndex][(*k1)].back()
893 && random <= fProbaShellMap[ionizationLevelIndex][(*k2)].back())
895 std::vector<double>::iterator prob12 =
896 std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
897 fProbaShellMap[ionizationLevelIndex][(*k1)].end(),
900 std::vector<double>::iterator prob11 = prob12 - 1;
902 std::vector<double>::iterator prob22 =
903 std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
904 fProbaShellMap[ionizationLevelIndex][(*k2)].end(),
907 std::vector<double>::iterator prob21 = prob22 - 1;
911 valuePROB21 = *prob21;
912 valuePROB22 = *prob22;
913 valuePROB12 = *prob12;
914 valuePROB11 = *prob11;
921 nrjTransf11 = fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
922 nrjTransf12 = fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
923 nrjTransf21 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
924 nrjTransf22 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
936 if (random > fProbaShellMap[ionizationLevelIndex][(*k1)].back())
938 std::vector<double>::iterator prob22 =
939 std::upper_bound(fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
940 fProbaShellMap[ionizationLevelIndex][(*k2)].end(),
943 std::vector<double>::iterator prob21 = prob22 - 1;
947 valuePROB21 = *prob21;
948 valuePROB22 = *prob22;
952 nrjTransf21 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
953 nrjTransf22 = fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
955 G4double interpolatedvalue2 = Interpolate(valuePROB21,
963 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
980 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
984 if (nrjTransfProduct != 0.)
986 nrj = QuadInterpolator(valuePROB11,
void set(double x, double y, double z)
static G4Electron * ElectronDefinition()
G4double LowEnergyLimit() const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
std::ostringstream G4ExceptionDescription
static constexpr double proton_mass_c2
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
static G4Proton * ProtonDefinition()
G4VEmAngularDistribution * GetAngularDistribution()
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
double B(double temperature)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4DNABornIonisationModel2(const G4ParticleDefinition *p=0, const G4String &nam="DNABornIonisationModel")
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double electron_mass_c2
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
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
static constexpr double m
const XML_Char int const XML_Char * value
const G4ThreeVector & GetMomentumDirection() const
static constexpr double cm
static constexpr double eV
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual ~G4DNABornIonisationModel2()
virtual size_t NumberOfComponents(void) const
G4double IonisationEnergy(G4int level)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &=*(new G4DataVector()))
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double TransferedEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
virtual G4double GetPartialCrossSection(const G4Material *, G4int, const G4ParticleDefinition *, G4double)
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
void SetAngularDistribution(G4VEmAngularDistribution *)
const G4Track * GetCurrentTrack() const
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static constexpr double MeV
G4VAtomDeexcitation * AtomDeexcitation()
void SetLowEnergyLimit(G4double)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
G4double bindingEnergy(G4int A, G4int Z)
static constexpr double keV
G4ParticleChangeForGamma * GetParticleChangeForGamma()
const G4Material * GetMaterial() const