53 G4cout <<
"PTB ionisation model is constructed " <<
G4endl;
83 G4cout <<
"Calling G4DNAPTBIonisationModel::Initialise()" <<
G4endl;
86 G4double scaleFactorBorn = (1.e-22 / 3.343) *
m*
m;
95 if(particle == electronDef)
103 "dna/sigma_ionisation_e-_PTB_THF",
104 "dna/sigmadiff_cumulated_ionisation_e-_PTB_THF",
111 "dna/sigma_ionisation_e-_PTB_PY",
112 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PY",
119 "dna/sigma_ionisation_e-_PTB_PU",
120 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU",
127 "dna/sigma_ionisation_e-_PTB_TMP",
128 "dna/sigmadiff_cumulated_ionisation_e-_PTB_TMP",
135 "dna/sigma_ionisation_e_born",
136 "dna/sigmadiff_ionisation_e_born",
145 "dna/sigma_ionisation_e-_PTB_THF",
146 "dna/sigmadiff_cumulated_ionisation_e-_PTB_THF",
153 "dna/sigma_ionisation_e-_PTB_PY",
154 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PY",
161 "dna/sigma_ionisation_e-_PTB_PY",
162 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PY",
169 "dna/sigma_ionisation_e-_PTB_PU",
170 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU",
177 "dna/sigma_ionisation_e-_PTB_PU",
178 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU",
185 "dna/sigma_ionisation_e-_PTB_TMP",
186 "dna/sigmadiff_cumulated_ionisation_e-_PTB_TMP",
192 else if (particle == protonDef)
200 "dna/sigma_ionisation_p_HKS_THF",
201 "dna/sigmadiff_cumulated_ionisation_p_PTB_THF",
209 "dna/sigma_ionisation_p_HKS_PY",
210 "dna/sigmadiff_cumulated_ionisation_p_PTB_PY",
227 "dna/sigma_ionisation_p_HKS_TMP",
228 "dna/sigmadiff_cumulated_ionisation_p_PTB_TMP",
258 G4cout <<
"Calling CrossSectionPerVolume() of G4DNAPTBIonisationModel" <<
G4endl;
271 if (ekin >= lowLim && ekin < highLim)
277 sigma = (*tableData)[materialName][particleName]->FindValue(ekin);
281 G4cout <<
"__________________________________" <<
G4endl;
282 G4cout <<
"°°° G4DNAPTBIonisationModel - XS INFO START" <<
G4endl;
283 G4cout <<
"°°° Kinetic energy(eV)=" << ekin/
eV <<
" particle : " << particleName <<
G4endl;
284 G4cout <<
"°°° Cross section per "<< materialName <<
" molecule (cm^2)=" << sigma/
cm/
cm <<
G4endl;
285 G4cout <<
"°°° G4DNAPTBIonisationModel - XS INFO END" <<
G4endl;
304 G4cout <<
"Calling SampleSecondaries() of G4DNAPTBIonisationModel" <<
G4endl;
317 if (k >= lowLim && k < highLim)
321 G4double totalEnergy = k + particleMass;
322 G4double pSquare = k * (totalEnergy + particleMass);
323 G4double totalMomentum = std::sqrt(pSquare);
334 if(materialName!=
"G4_WATER")
344 if(secondaryKinetic<=0)
346 G4cout<<
"Fatal error *************************************** "<<secondaryKinetic/
eV<<
G4endl;
358 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
359 G4double dirX = sinTheta*std::cos(phi);
360 G4double dirY = sinTheta*std::sin(phi);
363 deltaDirection.
rotateUz(primaryDirection);
374 G4double finalPx = totalMomentum*primaryDirection.
x() - deltaTotalMomentum*deltaDirection.
x();
375 G4double finalPy = totalMomentum*primaryDirection.
y() - deltaTotalMomentum*deltaDirection.
y();
376 G4double finalPz = totalMomentum*primaryDirection.
z() - deltaTotalMomentum*deltaDirection.
z();
377 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
378 finalPx /= finalMomentum;
379 finalPy /= finalMomentum;
380 finalPz /= finalMomentum;
385 G4cout<<
"Fatal error ****************************"<<
G4endl;
397 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
399 if(scatteredEnergy<=0)
401 G4cout<<
"Fatal error ****************************"<<
G4endl;
419 fvect->push_back(dp);
440 char *path = std::getenv(
"G4LEDATA");
444 G4Exception(
"G4DNAPTBIonisationModel::ReadAllDiffCSFiles",
"em0006",
450 std::ostringstream fullFileName;
451 fullFileName << path <<
"/"<< file<<
".dat";
454 std::ifstream diffCrossSection (fullFileName.str().c_str());
456 std::stringstream endPath;
457 if (!diffCrossSection)
459 endPath <<
"Missing data file: "<<
file;
460 G4Exception(
"G4DNAPTBIonisationModel::Initialise",
"em0003",
471 while(std::getline(diffCrossSection, line))
475 std::istringstream testIss(line);
485 else if(line.empty())
494 std::istringstream iss(line);
514 if(materialName!=
"G4_WATER")
518 fEnergySecondaryData[materialName][particleName][shell][
T][diffCrossSectionData[materialName][particleName][shell][
T][
E] ]=
E;
522 fProbaShellMap[materialName][particleName][shell][
T].push_back(diffCrossSectionData[materialName][particleName][shell][T][E]);
526 diffCrossSectionData[materialName][particleName][shell][
T][
E]*=scaleFactor;
565 G4double maxEnergy = maximumEnergyTransfer;
566 G4int nEnergySteps = 50;
568 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
574 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
581 G4double secondaryElectronKineticEnergy=0.;
590 return secondaryElectronKineticEnergy;
616 if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
619 G4double secondaryElectronKineticEnergy = 0.;
622 secondaryElectronKineticEnergy =
G4UniformRand() * maximumKineticEnergyTransfer;
626 return secondaryElectronKineticEnergy;
644 else if (secKinetic <= 200.*
eV)
652 cosTheta = std::sqrt(1.-sin2O);
665 if (secKinetic>100*
eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
676 G4int ionizationLevelIndex,
683 G4double kSE (energyTransfer-shellEnergy);
685 if (energyTransfer >= shellEnergy)
702 std::vector<double>::iterator
t2 = std::upper_bound(
fTMapWithVec[materialName][particleName].begin(),
fTMapWithVec[materialName][particleName].end(), k);
703 std::vector<double>::iterator
t1 = t2-1;
708 std::vector<double>::iterator e12 = std::upper_bound(
fEMapWithVector[materialName][particleName][(*t1)].begin(),
fEMapWithVector[materialName][particleName][(*t1)].end(), kSE);
709 std::vector<double>::iterator e11 = e12-1;
711 std::vector<double>::iterator e22 = std::upper_bound(
fEMapWithVector[materialName][particleName][(*t2)].begin(),
fEMapWithVector[materialName][particleName][(*t2)].end(), kSE);
712 std::vector<double>::iterator e21 = e22-1;
721 xs11 =
diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT1][valueE11];
722 xs12 =
diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT1][valueE12];
723 xs21 =
diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT2][valueE21];
724 xs22 =
diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT2][valueE22];
731 std::vector<double>::iterator
t2 = std::upper_bound(
fTMapWithVec[materialName][particleName].begin(),
fTMapWithVec[materialName][particleName].end(), k);
732 std::vector<double>::iterator
t1 = t2-1;
734 std::vector<double>::iterator e12 = std::upper_bound(
fEMapWithVector[materialName][particleName][(*t1)].begin(),
fEMapWithVector[materialName][particleName][(*t1)].end(), kSE);
735 std::vector<double>::iterator e11 = e12-1;
737 std::vector<double>::iterator e22 = std::upper_bound(
fEMapWithVector[materialName][particleName][(*t2)].begin(),
fEMapWithVector[materialName][particleName][(*t2)].end(), kSE);
738 std::vector<double>::iterator e21 = e22-1;
747 xs11 =
diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT1][valueE11];
748 xs12 =
diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT1][valueE12];
749 xs21 =
diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT2][valueE21];
750 xs22 =
diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT2][valueE22];
753 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
799 G4double ejectedElectronEnergy = 0.;
829 std::vector<double>::iterator
k2 = std::upper_bound(
fTMapWithVec[materialName][particleName].begin(),
fTMapWithVec[materialName][particleName].end(), k);
830 std::vector<double>::iterator
k1 = k2-1;
840 G4cerr<<
"**************** Fatal error ******************"<<
G4endl;
841 G4cerr<<
"G4DNAPTBIonisationModel::RandomizeEjectedElectronEnergyFromCumulated"<<
G4endl;
842 G4cerr<<
"You have *k1 > *k2 with k1 "<<*k1<<
" and k2 "<<*k2<<
G4endl;
843 G4cerr<<
"This may be because the energy of the incident particle is to high for the data table."<<
G4endl;
853 std::vector<double>::iterator cumulCS12 = std::upper_bound(
fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k1)].begin(),
854 fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k1)].end(), random);
855 std::vector<double>::iterator cumulCS11 = cumulCS12-1;
857 std::vector<double>::iterator cumulCS22 = std::upper_bound(
fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k2)].begin(),
858 fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k2)].end(), random);
859 std::vector<double>::iterator cumulCS21 = cumulCS22-1;
864 valueCumulCS11 = *cumulCS11;
865 valueCumulCS12 = *cumulCS12;
866 valueCumulCS21 = *cumulCS21;
867 valueCumulCS22 = *cumulCS22;
883 if(cumulCS12==
fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k1)].end())
888 secElecE21 =
fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK2][valueCumulCS21];
889 secElecE22 =
fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK2][valueCumulCS22];
897 secElecE11 =
fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK1][valueCumulCS11];
898 secElecE12 =
fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK1][valueCumulCS12];
899 secElecE21 =
fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK2][valueCumulCS21];
900 secElecE22 =
fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK2][valueCumulCS22];
904 valueCumulCS21, valueCumulCS22,
905 secElecE11, secElecE12,
906 secElecE21, secElecE22,
915 if(k-ejectedElectronEnergy-
bindingEnergy<=0 || ejectedElectronEnergy<=0)
924 G4cout<<
"surrounding k values: valueK1 valueK2\n"<<valueK1<<
" "<<valueK2<<
G4endl;
925 G4cout<<
"surrounding E values: secElecE11 secElecE12 secElecE21 secElecE22\n"
926 <<secElecE11<<
" "<<secElecE12<<
" "<<secElecE21<<
" "<<secElecE22<<
" "<<
G4endl;
927 G4cout<<
"surrounding cumulCS values: valueCumulCS11 valueCumulCS12 valueCumulCS21 valueCumulCS22\n"
928 <<valueCumulCS11<<
" "<<valueCumulCS12<<
" "<<valueCumulCS21<<
" "<<valueCumulCS22<<
" "<<
G4endl;
934 return ejectedElectronEnergy*
eV;
949 if ((e2-e1)!=0 && xs1 !=0 && xs2 !=0)
953 value = std::pow(10.,(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)) );
959 if ((e2-e1)!=0 && (xs1 ==0 || xs2 ==0))
963 value = (d1 + (d2 -
d1)*(e - e1)/ (e2 -
e1));
980 else interpolatedvalue1 = xs11;
984 else interpolatedvalue2 = xs21;
987 if(interpolatedvalue1!=interpolatedvalue2) value =
LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
988 else value = interpolatedvalue1;