60 G4cout <<
"Born ionisation model is constructed " <<
G4endl;
65 fAtomDeexcitation = 0;
67 fpMolWaterDensity = 0;
91 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator
pos;
92 for (pos = tableData.begin(); pos != tableData.end(); ++
pos)
110 if (verboseLevel > 3)
112 G4cout <<
"Calling G4DNABornIonisationModel1::Initialise()" <<
G4endl;
117 G4String fileElectron(
"dna/sigma_ionisation_e_born");
118 G4String fileProton(
"dna/sigma_ionisation_p_born");
126 G4double scaleFactor = (1.e-22 / 3.343) *
m*
m;
128 char *path = getenv(
"G4LEDATA");
148 std::ostringstream eFullFileName;
150 if (fasterCode) eFullFileName << path <<
"/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
151 if (!fasterCode) eFullFileName << path <<
"/dna/sigmadiff_ionisation_e_born.dat";
153 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
155 if (!eDiffCrossSection)
157 if (fasterCode)
G4Exception(
"G4DNABornIonisationModel1::Initialise",
"em0003",
158 FatalException,
"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat");
160 if (!fasterCode)
G4Exception(
"G4DNABornIonisationModel1::Initialise",
"em0003",
161 FatalException,
"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
175 for (
int j=0; j<5; j++)
177 eProbaShellMap[j].clear();
178 pProbaShellMap[j].clear();
180 eDiffCrossSectionData[j].clear();
181 pDiffCrossSectionData[j].clear();
183 eNrjTransfData[j].clear();
184 pNrjTransfData[j].clear();
188 eTdummyVec.push_back(0.);
189 while(!eDiffCrossSection.eof())
193 eDiffCrossSection>>tDummy>>eDummy;
194 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
197 for (
int j=0; j<5; j++)
199 eDiffCrossSection>> tmp;
201 eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
205 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
206 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
210 if (!eDiffCrossSection.eof() && !fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
212 if (!fasterCode) eVecm[tDummy].push_back(eDummy);
221 tableFile[
proton] = fileProton;
231 tableData[
proton] = tableP;
235 std::ostringstream pFullFileName;
237 if (fasterCode) pFullFileName << path <<
"/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
239 if (!fasterCode) pFullFileName << path <<
"/dna/sigmadiff_ionisation_p_born.dat";
241 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
243 if (!pDiffCrossSection)
245 if (fasterCode)
G4Exception(
"G4DNABornIonisationModel1::Initialise",
"em0003",
246 FatalException,
"Missing data file:/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat");
248 if (!fasterCode)
G4Exception(
"G4DNABornIonisationModel1::Initialise",
"em0003",
249 FatalException,
"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
252 pTdummyVec.push_back(0.);
253 while(!pDiffCrossSection.eof())
257 pDiffCrossSection>>tDummy>>eDummy;
258 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
259 for (
int j=0; j<5; j++)
261 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
265 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
266 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
270 if (!pDiffCrossSection.eof() && !fasterCode) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
272 if (!fasterCode) pVecm[tDummy].push_back(eDummy);
278 if (particle==electronDef)
284 if (particle==protonDef)
292 G4cout <<
"Born ionisation model is initialized " <<
G4endl
310 isInitialised =
true;
321 if (verboseLevel > 3)
323 G4cout <<
"Calling CrossSectionPerVolume() of G4DNABornIonisationModel1"
344 if(waterDensity!= 0.0)
349 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
350 pos1 = lowEnergyLimit.find(particleName);
351 if (pos1 != lowEnergyLimit.end())
353 lowLim = pos1->second;
356 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
357 pos2 = highEnergyLimit.find(particleName);
358 if (pos2 != highEnergyLimit.end())
360 highLim = pos2->second;
363 if (ekin >= lowLim && ekin < highLim)
365 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator
pos;
366 pos = tableData.find(particleName);
368 if (pos != tableData.end())
379 G4double A = 1.39241700556072800000E-009 ;
380 G4double B = -8.52610412942622630000E-002 ;
381 sigma = sigma *
G4Exp(A*(ekin/
eV)+B);
389 G4Exception(
"G4DNABornIonisationModel1::CrossSectionPerVolume",
"em0002",
394 if (verboseLevel > 2)
396 G4cout <<
"__________________________________" <<
G4endl;
397 G4cout <<
"G4DNABornIonisationModel1 - XS INFO START" <<
G4endl;
398 G4cout <<
"Kinetic energy(eV)=" << ekin/
eV <<
" particle : " << particleName <<
G4endl;
399 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
400 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
401 G4cout <<
"G4DNABornIonisationModel1 - XS INFO END" <<
G4endl;
405 return sigma*waterDensity;
417 if (verboseLevel > 3)
419 G4cout <<
"Calling SampleSecondaries() of G4DNABornIonisationModel1"
430 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
431 pos1 = lowEnergyLimit.find(particleName);
433 if (pos1 != lowEnergyLimit.end())
435 lowLim = pos1->second;
438 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
439 pos2 = highEnergyLimit.find(particleName);
441 if (pos2 != highEnergyLimit.end())
443 highLim = pos2->second;
446 if (k >= lowLim && k < highLim)
450 G4double totalEnergy = k + particleMass;
451 G4double pSquare = k * (totalEnergy + particleMass);
452 G4double totalMomentum = std::sqrt(pSquare);
454 G4int ionizationShell = 0;
456 if (!fasterCode) ionizationShell = RandomSelect(k,particleName);
467 ionizationShell = RandomSelect(k,particleName);
474 G4int secNumberInit = 0;
475 G4int secNumberFinal = 0;
481 if (k<bindingEnergy)
return;
485 if(fAtomDeexcitation)
489 if (ionizationShell <5 && ionizationShell >1)
493 else if (ionizationShell <2)
508 secNumberInit = fvect->size();
510 secNumberFinal = fvect->size();
515 if (fasterCode ==
false)
517 secondaryKinetic = RandomizeEjectedElectronEnergy(particle->
GetDefinition(),k,ionizationShell);
522 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->
GetDefinition(),k,ionizationShell);
535 G4double finalPx = totalMomentum*primaryDirection.
x() - deltaTotalMomentum*deltaDirection.
x();
536 G4double finalPy = totalMomentum*primaryDirection.
y() - deltaTotalMomentum*deltaDirection.
y();
537 G4double finalPz = totalMomentum*primaryDirection.
z() - deltaTotalMomentum*deltaDirection.
z();
538 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
539 finalPx /= finalMomentum;
540 finalPy /= finalMomentum;
541 finalPz /= finalMomentum;
544 direction.
set(finalPx,finalPy,finalPz);
552 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
554 for (
G4int j=secNumberInit; j < secNumberFinal; j++)
556 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
571 if (secondaryKinetic>0)
574 fvect->push_back(dp);
595 G4double maximumEnergyTransfer = 0.;
597 maximumEnergyTransfer = k;
615 G4double maxEnergy = maximumEnergyTransfer;
616 G4int nEnergySteps = 50;
620 1. / static_cast<G4double>(nEnergySteps - 1)));
621 G4int step(nEnergySteps);
630 if (differentialCrossSection >= crossSectionMaximum)
631 crossSectionMaximum = differentialCrossSection;
636 G4double secondaryElectronKineticEnergy = 0.;
644 return secondaryElectronKineticEnergy;
650 G4double maximumKineticEnergyTransfer = 4.
662 if (differentialCrossSection >= crossSectionMaximum)
663 crossSectionMaximum = differentialCrossSection;
666 G4double secondaryElectronKineticEnergy = 0.;
669 secondaryElectronKineticEnergy =
G4UniformRand()* maximumKineticEnergyTransfer;
674 return secondaryElectronKineticEnergy;
728 G4int ionizationLevelIndex)
732 if (energyTransfer >= waterStructure.
IonisationEnergy(ionizationLevelIndex))
750 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
754 std::vector<double>::iterator t1 = t2 - 1;
757 if (energyTransfer <= eVecm[(*t1)].back()
758 && energyTransfer <= eVecm[(*t2)].back())
760 std::vector<double>::iterator e12 =
761 std::upper_bound(eVecm[(*t1)].begin(),
764 std::vector<double>::iterator e11 = e12 - 1;
766 std::vector<double>::iterator e22 =
767 std::upper_bound(eVecm[(*t2)].begin(),
770 std::vector<double>::iterator e21 = e22 - 1;
779 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
780 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
781 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
782 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
791 std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),
794 std::vector<double>::iterator t1 = t2 - 1;
796 std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),
799 std::vector<double>::iterator e11 = e12 - 1;
801 std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),
804 std::vector<double>::iterator e21 = e22 - 1;
813 xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
814 xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
815 xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
816 xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
820 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
823 sigma = QuadInterpolator(valueE11,
853 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
856 G4double a = (std::log10(xs2) - std::log10(xs1))
857 / (std::log10(e2) - std::log10(e1));
858 G4double b = std::log10(xs2) - a * std::log10(e2);
859 G4double sigma = a * std::log10(e) + b;
860 value = (std::pow(10., sigma));
874 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
878 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
884 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
888 value = (d1 + (d2 -
d1) * (e - e1) / (e2 - e1));
920 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
921 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
936 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator
pos;
939 if (pos != tableData.end())
955 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator
pos;
956 pos = tableData.find(particle);
958 if (pos != tableData.end())
973 value += valuesBuffer[i];
984 if (valuesBuffer[i] > value)
986 delete[] valuesBuffer;
989 value -= valuesBuffer[i];
993 delete[] valuesBuffer;
998 G4Exception(
"G4DNABornIonisationModel1::RandomSelect",
1001 "Model not applicable to particle type.");
1015 G4double secondaryElectronKineticEnergy = 0.;
1027 if (secondaryElectronKineticEnergy < 0.)
1031 return secondaryElectronKineticEnergy;
1038 G4int ionizationLevelIndex,
1058 std::vector<double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),
1061 std::vector<double>::iterator k1 = k2 - 1;
1075 if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1076 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
1078 std::vector<double>::iterator prob12 =
1079 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1080 eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1083 std::vector<double>::iterator prob11 = prob12 - 1;
1085 std::vector<double>::iterator prob22 =
1086 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1087 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1090 std::vector<double>::iterator prob21 = prob22 - 1;
1094 valuePROB21 = *prob21;
1095 valuePROB22 = *prob22;
1096 valuePROB12 = *prob12;
1097 valuePROB11 = *prob11;
1104 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1105 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1106 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1107 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1119 if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
1121 std::vector<double>::iterator prob22 =
1122 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1123 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1126 std::vector<double>::iterator prob21 = prob22 - 1;
1130 valuePROB21 = *prob21;
1131 valuePROB22 = *prob22;
1135 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1136 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1138 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1146 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1166 std::vector<double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),
1170 std::vector<double>::iterator k1 = k2 - 1;
1185 if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1186 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1188 std::vector<double>::iterator prob12 =
1189 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1190 pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1193 std::vector<double>::iterator prob11 = prob12 - 1;
1195 std::vector<double>::iterator prob22 =
1196 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1197 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1200 std::vector<double>::iterator prob21 = prob22 - 1;
1204 valuePROB21 = *prob21;
1205 valuePROB22 = *prob22;
1206 valuePROB12 = *prob12;
1207 valuePROB11 = *prob11;
1214 nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1215 nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1216 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1217 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1230 if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1232 std::vector<double>::iterator prob22 =
1233 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1234 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1237 std::vector<double>::iterator prob21 = prob22 - 1;
1241 valuePROB21 = *prob21;
1242 valuePROB22 = *prob22;
1246 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1247 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1249 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1257 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1274 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1278 if (nrjTransfProduct != 0.)
1280 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()
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)
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
G4ParticleChangeForGamma * fParticleChangeForGamma
G4double TransferedEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
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
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual G4double GetPartialCrossSection(const G4Material *, G4int, const G4ParticleDefinition *, G4double)
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.
G4DNABornIonisationModel1(const G4ParticleDefinition *p=0, const G4String &nam="DNABornIonisationModel")
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
virtual ~G4DNABornIonisationModel1()
void SetAngularDistribution(G4VEmAngularDistribution *)
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
const G4Track * GetCurrentTrack() const
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static constexpr double MeV
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
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
static const G4double pos
G4ParticleChangeForGamma * GetParticleChangeForGamma()
const G4Material * GetMaterial() const