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 = 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)));
569 G4int step(nEnergySteps);
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;
TriDimensionMap diffCrossSectionData
virtual void Initialise(const G4ParticleDefinition *particle, const G4DataVector &=*(new G4DataVector()), G4ParticleChangeForGamma *fpChangeForGamme=nullptr)
Initialise Method called once at the beginning of the simulation. It is used to setup the list of the...
The G4DNAPTBAugerModel class Implement the PTB Auger model.
void SetHighELimit(const G4String &material, const G4String &particle, G4double lim)
SetHighEnergyLimit.
TriDimensionMap fEnergySecondaryData
static constexpr double MeV
virtual ~G4DNAPTBIonisationModel()
~G4DNAPTBIonisationModel Destructor
static constexpr double keV
void ComputeAugerEffect(std::vector< G4DynamicParticle * > *fvect, const G4String &materialNameIni, G4double bindingEnergy)
ComputeAugerEffect Main method to be called by the ionisation model.
const G4ThreeVector & GetMomentumDirection() const
const G4String & GetParticleName() const
G4DNAPTBIonisationStructure ptbStructure
G4double GetHighELimit(const G4String &material, const G4String &particle)
GetHighEnergyLimit.
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Hep3Vector & rotateUz(const Hep3Vector &)
void AddCrossSectionData(G4String materialName, G4String particleName, G4String fileCS, G4String fileDiffCS, G4double scaleFactor)
AddCrossSectionData Method used during the initialization of the model class to add a new material...
G4double GetPDGMass() const
std::map< G4String, std::map< G4String, std::vector< double > > > fTMapWithVec
void LoadCrossSectionData(const G4String &particleName)
LoadCrossSectionData Method to loop on all the registered materials in the model and load the corresp...
static constexpr double proton_mass_c2
static constexpr double m
G4double RandomizeEjectedElectronEnergyFromCumulated(G4ParticleDefinition *particleDefinition, G4double k, G4int shell, const G4String &materialName)
RandomizeEjectedElectronEnergyFromCumulated Uses the cumulated tables to find the energy of the eject...
G4double GetLowELimit(const G4String &material, const G4String &particle)
GetLowEnergyLimit.
G4ParticleDefinition * GetDefinition() const
VecMapWithShell fProbaShellMap
void SetLowELimit(const G4String &material, const G4String &particle, G4double lim)
SetLowEnergyLimit.
const XML_Char int const XML_Char * value
static G4Proton * ProtonDefinition()
static constexpr double electron_mass_c2
G4DNAPTBIonisationModel(const G4String &applyToMaterial="all", const G4ParticleDefinition *p=0, const G4String &nam="DNAPTBIonisationModel", const G4bool isAuger=true)
G4DNAPTBIonisationModel Constructor.
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
static constexpr double twopi
static constexpr double eV
G4double IonisationEnergy(G4int level, const G4String &materialName)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4String &materialName, const G4DynamicParticle *, G4ParticleChangeForGamma *particleChangeForGamma, G4double tmin, G4double tmax)
SampleSecondaries If the model is selected for the ModelInterface then SampleSecondaries will be call...
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell, const G4String &materialName)
G4int NumberOfLevels(const G4String &materialName)
static G4Electron * Electron()
G4GLOB_DLL std::ostream G4cerr
G4int verboseLevel
verbose level
static G4Electron * ElectronDefinition()
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4String &materialName, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
CrossSectionPerVolume Mandatory for every model the CrossSectionPerVolume method is in charge of retu...
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void ReadDiffCSFile(const G4String &materialName, const G4String &particleName, const G4String &file, const G4double scaleFactor)
ReadDiffCSFile Method to read the differential cross section files.
TableMapData * GetTableData()
GetTableData.
G4DNAPTBAugerModel * fDNAPTBAugerModel
PTB Auger model instanciated in the constructor and deleted in the destructor of the class...
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
LogLogInterpolate.
void RandomizeEjectedElectronDirection(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4double outgoingParticleEnergy, G4double &cosTheta, G4double &phi)
RandomizeEjectedElectronDirection Method to calculate the ejected electron direction.
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double xs11, G4double xs12, G4double xs21, G4double xs22, G4double t1, G4double t2, G4double t, G4double e)
QuadInterpolator.
G4double GetKineticEnergy() const
static constexpr double cm
G4GLOB_DLL std::ostream G4cout
virtual void Initialise()
Initialise Set the verbose value.
G4int RandomSelectShell(G4double k, const G4String &particle, const G4String &materialName)
RandomSelectShell Method to randomely select a shell from the data table uploaded. The size of the table (number of columns) is used to determine the total number of possible shells.
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double bindingEnergy(G4int A, G4int Z)
std::map< G4String, std::map< G4String, G4DNACrossSectionDataSet *, std::less< G4String > > > TableMapData
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, const G4String &materialName)