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;
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")
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))
980 else interpolatedvalue1 = xs11;
984 else interpolatedvalue2 = xs21;
987 if(interpolatedvalue1!=interpolatedvalue2) value =
LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
988 else value = interpolatedvalue1;
static const G4double e1[44]
static const G4double e2[44]
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
static constexpr double m
static constexpr double twopi
static constexpr double keV
static constexpr double eV
static constexpr double MeV
static constexpr double cm
G4GLOB_DLL std::ostream G4cerr
G4GLOB_DLL std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
The G4DNAPTBAugerModel class Implement the PTB Auger model.
void ComputeAugerEffect(std::vector< G4DynamicParticle * > *fvect, const G4String &materialNameIni, G4double bindingEnergy)
ComputeAugerEffect Main method to be called by the ionisation model.
virtual void Initialise()
Initialise Set the verbose value.
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...
VecMapWithShell fProbaShellMap
void ReadDiffCSFile(const G4String &materialName, const G4String &particleName, const G4String &file, const G4double scaleFactor)
ReadDiffCSFile Method to read the differential cross section files.
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...
std::map< G4String, std::map< G4String, std::vector< double > > > fTMapWithVec
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, const G4String &materialName)
TriDimensionMap diffCrossSectionData
G4DNAPTBIonisationModel(const G4String &applyToMaterial="all", const G4ParticleDefinition *p=0, const G4String &nam="DNAPTBIonisationModel", const G4bool isAuger=true)
G4DNAPTBIonisationModel Constructor.
TriDimensionMap fEnergySecondaryData
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell, const G4String &materialName)
G4double RandomizeEjectedElectronEnergyFromCumulated(G4ParticleDefinition *particleDefinition, G4double k, G4int shell, const G4String &materialName)
RandomizeEjectedElectronEnergyFromCumulated Uses the cumulated tables to find the energy of the eject...
G4DNAPTBAugerModel * fDNAPTBAugerModel
PTB Auger model instanciated in the constructor and deleted in the destructor of the class.
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...
G4int verboseLevel
verbose level
void RandomizeEjectedElectronDirection(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4double outgoingParticleEnergy, G4double &cosTheta, G4double &phi)
RandomizeEjectedElectronDirection Method to calculate the ejected electron direction.
virtual ~G4DNAPTBIonisationModel()
~G4DNAPTBIonisationModel Destructor
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.
G4DNAPTBIonisationStructure ptbStructure
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
LogLogInterpolate.
G4double IonisationEnergy(G4int level, const G4String &materialName)
G4int NumberOfLevels(const G4String &materialName)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * ElectronDefinition()
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double GetPDGMass() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
G4int RandomSelectShell(G4double k, const G4String &particle, const G4String &materialName)
RandomSelectShell Method to randomely select a shell from the data table uploaded....
TableMapData * GetTableData()
GetTableData.
void SetHighELimit(const G4String &material, const G4String &particle, G4double lim)
SetHighEnergyLimit.
std::map< G4String, std::map< G4String, G4DNACrossSectionDataSet *, std::less< G4String > > > TableMapData
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....
void SetLowELimit(const G4String &material, const G4String &particle, G4double lim)
SetLowEnergyLimit.
G4double GetHighELimit(const G4String &material, const G4String &particle)
GetHighEnergyLimit.
void LoadCrossSectionData(const G4String &particleName)
LoadCrossSectionData Method to loop on all the registered materials in the model and load the corresp...
G4double GetLowELimit(const G4String &material, const G4String &particle)
GetLowEnergyLimit.
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double bindingEnergy(G4int A, G4int Z)