59 G4cout <<
"Born ionisation model is constructed " <<
G4endl;
112 G4cout <<
"Calling G4DNABornIonisationModel2::Initialise()" <<
G4endl;
118 description <<
"You are trying to initialized G4DNABornIonisationModel2 "
122 description <<
"G4DNABornIonisationModel2 was already initiliased "
124 G4Exception(
"G4DNABornIonisationModel2::Initialise",
"bornIonInit",
131 char *path = getenv(
"G4LEDATA");
136 std::ostringstream fullFileName;
137 fullFileName << path;
139 if(particleName ==
"e-")
147 fullFileName <<
"/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
151 fullFileName <<
"/dna/sigmadiff_ionisation_e_born.dat";
154 else if(particleName ==
"proton")
162 fullFileName <<
"/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
166 fullFileName <<
"/dna/sigmadiff_ionisation_p_born.dat";
172 G4double scaleFactor = (1.e-22 / 3.343) *
m*
m;
178 std::ifstream diffCrossSection(fullFileName.str().c_str());
180 if (!diffCrossSection)
183 description <<
"Missing data file:" <<
G4endl << fullFileName.str() <<
G4endl;
184 G4Exception(
"G4DNABornIonisationModel2::Initialise",
"em0003",
194 for (
int j=0; j<5; j++)
204 while(!diffCrossSection.eof())
208 diffCrossSection>>tDummy>>eDummy;
212 for (
int j=0; j<5; j++)
214 diffCrossSection>>
tmp;
221 fProbaShellMap[j][tDummy].push_back(fDiffCrossSectionData[j][tDummy][eDummy]);
238 G4cout <<
"Born ionisation model is initialized " <<
G4endl
271 G4cout <<
"Calling CrossSectionPerVolume() of G4DNABornIonisationModel2"
291 G4double A = 1.39241700556072800000E-009 ;
292 G4double B = -8.52610412942622630000E-002 ;
293 sigma = sigma *
G4Exp(A*(ekin/
eV)+B);
300 G4cout <<
"__________________________________" <<
G4endl;
301 G4cout <<
"G4DNABornIonisationModel2 - XS INFO START" <<
G4endl;
303 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
304 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./
cm) << G4endl;
305 G4cout <<
"G4DNABornIonisationModel2 - XS INFO END" <<
G4endl;
308 return sigma*waterDensity;
322 G4cout <<
"Calling SampleSecondaries() of G4DNABornIonisationModel2"
332 G4double totalEnergy = k + particleMass;
333 G4double pSquare = k * (totalEnergy + particleMass);
334 G4double totalMomentum = std::sqrt(pSquare);
336 G4int ionizationShell = 0;
370 if (secondaryKinetic>0)
373 fvect->push_back(dp);
380 G4double finalPx = totalMomentum*primaryDirection.
x() - deltaTotalMomentum*deltaDirection.
x();
381 G4double finalPy = totalMomentum*primaryDirection.
y() - deltaTotalMomentum*deltaDirection.
y();
382 G4double finalPz = totalMomentum*primaryDirection.
z() - deltaTotalMomentum*deltaDirection.
z();
383 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
384 finalPx /= finalMomentum;
385 finalPy /= finalMomentum;
386 finalPz /= finalMomentum;
389 direction.
set(finalPx,finalPy,finalPz);
400 size_t secNumberInit = 0;
401 size_t secNumberFinal = 0;
407 if (k<bindingEnergy)
return;
410 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
417 secNumberInit = fvect->size();
419 secNumberFinal = fvect->size();
421 if(secNumberFinal > secNumberInit)
423 for (
size_t i=secNumberInit; i<secNumberFinal; ++i)
426 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
429 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
444 if(bindingEnergy < 0.0)
445 G4Exception(
"G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
485 G4double maximumEnergyTransfer = 0.;
487 maximumEnergyTransfer =
k;
505 G4double maxEnergy = maximumEnergyTransfer;
506 G4int nEnergySteps = 50;
509 G4double stpEnergy(std::pow(maxEnergy / value,
510 1. / static_cast<G4double>(nEnergySteps - 1)));
520 if (differentialCrossSection >= crossSectionMaximum)
521 crossSectionMaximum = differentialCrossSection;
526 G4double secondaryElectronKineticEnergy = 0.;
534 return secondaryElectronKineticEnergy;
540 G4double maximumKineticEnergyTransfer = 4.
552 if (differentialCrossSection >= crossSectionMaximum)
553 crossSectionMaximum = differentialCrossSection;
556 G4double secondaryElectronKineticEnergy = 0.;
559 secondaryElectronKineticEnergy =
G4UniformRand()* maximumKineticEnergyTransfer;
564 return secondaryElectronKineticEnergy;
618 G4int ionizationLevelIndex)
642 std::vector<G4double>::iterator
t2 = std::upper_bound(
fTdummyVec.begin(),
646 std::vector<G4double>::iterator
t1 = t2 - 1;
650 if (energyTransfer <=
fVecm[(*t1)].back()
651 && energyTransfer <=
fVecm[(*t2)].back())
653 std::vector<G4double>::iterator e12 = std::upper_bound(
fVecm[(*t1)].begin(),
656 std::vector<G4double>::iterator e11 = e12 - 1;
658 std::vector<G4double>::iterator e22 = std::upper_bound(
fVecm[(*t2)].begin(),
661 std::vector<G4double>::iterator e21 = e22 - 1;
677 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
710 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
713 G4double a = (std::log10(xs2) - std::log10(xs1))
714 / (std::log10(e2) - std::log10(e1));
715 G4double b = std::log10(xs2) - a * std::log10(e2);
717 value = (std::pow(10., sigma));
731 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 &&
fasterCode)
735 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
741 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) &&
fasterCode)
745 value = (d1 + (d2 -
d1) * (e - e1) / (e2 -
e1));
813 value += valuesBuffer[i];
824 if (valuesBuffer[i] > value)
826 delete[] valuesBuffer;
829 value -= valuesBuffer[i];
833 delete[] valuesBuffer;
846 G4double secondaryElectronKineticEnergy = 0.;
857 if (secondaryElectronKineticEnergy < 0.)
861 return secondaryElectronKineticEnergy;
868 G4int ionizationLevelIndex,
891 std::vector<G4double>::iterator
k2 = std::upper_bound(
fTdummyVec.begin(),
894 std::vector<G4double>::iterator
k1 = k2 - 1;
911 std::vector<G4double>::iterator prob12 =
912 std::upper_bound(
fProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
916 std::vector<G4double>::iterator prob11 = prob12 - 1;
918 std::vector<G4double>::iterator prob22 =
919 std::upper_bound(
fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
923 std::vector<G4double>::iterator prob21 = prob22 - 1;
927 valuePROB21 = *prob21;
928 valuePROB22 = *prob22;
929 valuePROB12 = *prob12;
930 valuePROB11 = *prob11;
937 nrjTransf11 =
fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
938 nrjTransf12 =
fNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
939 nrjTransf21 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
940 nrjTransf22 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
954 std::vector<G4double>::iterator prob22 =
955 std::upper_bound(
fProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
959 std::vector<G4double>::iterator prob21 = prob22 - 1;
963 valuePROB21 = *prob21;
964 valuePROB22 = *prob22;
968 nrjTransf21 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
969 nrjTransf22 =
fNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
996 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1000 if (nrjTransfProduct != 0.)