85 G4cout <<
"MicroElec inelastic model is constructed " <<
G4endl;
90 fAtomDeexcitation =
nullptr;
105 for (
auto & pos : tableData)
123 if (verboseLevel > 3)
124 G4cout <<
"Calling G4MicroElecInelasticModel::Initialise()" <<
G4endl;
128 G4String fileElectron(
"microelec/sigma_inelastic_e_Si");
129 G4String fileProton(
"microelec/sigma_inelastic_p_Si");
136 G4double scaleFactor = 1e-18 * cm *cm;
141 tableFile[electron] = fileElectron;
142 lowEnergyLimit[electron] = 16.7 * eV;
143 highEnergyLimit[electron] = 100.0 * MeV;
149 tableData[electron] = tableE;
153 std::ostringstream eFullFileName;
155 if (fasterCode) eFullFileName << path <<
"/microelec/sigmadiff_cumulated_inelastic_e_Si.dat";
156 else eFullFileName << path <<
"/microelec/sigmadiff_inelastic_e_Si.dat";
158 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
160 if (!eDiffCrossSection)
162 if (fasterCode)
G4Exception(
"G4MicroElecInelasticModel::Initialise",
"em0003",
163 FatalException,
"Missing data file:/microelec/sigmadiff_cumulated_inelastic_e_Si.dat");
165 else G4Exception(
"G4MicroElecInelasticModel::Initialise",
"em0003",
166 FatalException,
"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
177 for (
int j=0; j<6; j++)
179 eProbaShellMap[j].clear();
180 pProbaShellMap[j].clear();
182 eDiffCrossSectionData[j].clear();
183 pDiffCrossSectionData[j].clear();
185 eNrjTransfData[j].clear();
186 pNrjTransfData[j].clear();
190 eTdummyVec.push_back(0.);
191 while(!eDiffCrossSection.eof())
195 eDiffCrossSection>>tDummy>>eDummy;
196 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
199 for (
int j=0; j<6; j++)
201 eDiffCrossSection>> tmp;
203 eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
207 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
208 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
213 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
214 eVecm[tDummy].push_back(eDummy);
222 tableFile[proton] = fileProton;
223 lowEnergyLimit[proton] = 50. * keV;
224 highEnergyLimit[proton] = 10. * GeV;
229 tableData[proton] = tableP;
232 std::ostringstream pFullFileName;
234 if (fasterCode) pFullFileName << path <<
"/microelec/sigmadiff_cumulated_inelastic_p_Si.dat";
235 else pFullFileName << path <<
"/microelec/sigmadiff_inelastic_p_Si.dat";
237 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
239 if (!pDiffCrossSection)
241 if (fasterCode)
G4Exception(
"G4MicroElecInelasticModel::Initialise",
"em0003",
242 FatalException,
"Missing data file:/microelec/sigmadiff_cumulated_inelastic_p_Si.dat");
244 else G4Exception(
"G4MicroElecInelasticModel::Initialise",
"em0003",
245 FatalException,
"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
248 pTdummyVec.push_back(0.);
249 while(!pDiffCrossSection.eof())
253 pDiffCrossSection>>tDummy>>eDummy;
254 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
255 for (
int j=0; j<6; j++)
257 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
261 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
262 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
267 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
268 pVecm[tDummy].push_back(eDummy);
273 if (particle==electronDef)
279 if (particle==protonDef)
287 G4cout <<
"MicroElec Inelastic model is initialized " <<
G4endl
292 <<
" with mass (amu) " << particle->
GetPDGMass()/proton_mass_c2
300 if (isInitialised) {
return; }
302 isInitialised =
true;
313 if (verboseLevel > 3)
314 G4cout <<
"Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" <<
G4endl;
329 if (Mion_c2 > proton_mass_c2)
335 if (verboseLevel > 3)
337 <<
"Particle : " << nameLocal <<
", mass : " << Mion_c2/proton_mass_c2 <<
"*mp, charge " << Zeff
338 <<
", Ekin (eV) = " << ekin/eV <<
G4endl ;
340 ekin *= proton_mass_c2/Mion_c2 ;
341 nameLocal =
"proton" ;
343 if (verboseLevel > 3)
345 <<
"Particle : " << nameLocal <<
", Ekin (eV) = " << ekin/eV <<
G4endl ;
350 auto pos1 = lowEnergyLimit.find(nameLocal);
351 if (pos1 != lowEnergyLimit.end())
353 lowLim = pos1->second;
356 auto pos2 = highEnergyLimit.find(nameLocal);
357 if (pos2 != highEnergyLimit.end())
359 highLim = pos2->second;
362 if (ekin >= lowLim && ekin < highLim)
364 auto pos = tableData.find(nameLocal);
365 if (pos != tableData.end())
368 if (table !=
nullptr)
375 G4Exception(
"G4MicroElecInelasticModel::CrossSectionPerVolume",
"em0002",
388 if (verboseLevel > 3)
390 G4cout <<
"---> Kinetic energy (eV)=" << ekin/eV <<
G4endl;
391 G4cout <<
" - Cross section per Si atom (cm^2)=" << sigma*Zeff2/cm2 <<
G4endl;
392 G4cout <<
" - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./cm) <<
G4endl;
396 return sigma*density*Zeff2;
408 if (verboseLevel > 3)
409 G4cout <<
"Calling SampleSecondaries() of G4MicroElecInelasticModel" <<
G4endl;
419 G4String nameLocal2 = particleName ;
422 if (particleMass > proton_mass_c2)
424 k *= proton_mass_c2/particleMass ;
426 nameLocal2 =
"proton" ;
429 auto pos1 = lowEnergyLimit.find(nameLocal2);
430 if (pos1 != lowEnergyLimit.end())
432 lowLim = pos1->second;
435 auto pos2 = highEnergyLimit.find(nameLocal2);
436 if (pos2 != highEnergyLimit.end())
438 highLim = pos2->second;
441 if (k >= lowLim && k < highLim)
444 G4double totalEnergy = ekin + particleMass;
445 G4double pSquare = ekin * (totalEnergy + particleMass);
446 G4double totalMomentum = std::sqrt(pSquare);
450 Shell = RandomSelect(k,nameLocal2);
460 if (verboseLevel > 3)
462 G4cout <<
"---> Kinetic energy (eV)=" << k/eV <<
G4endl ;
463 G4cout <<
"Shell: " << Shell <<
", energy: " << bindingEnergy/eV <<
G4endl;
467 std::size_t secNumberInit = 0;
468 std::size_t secNumberFinal = 0;
471 if (k<bindingEnergy)
return;
475 if(fAtomDeexcitation && Shell > 2) {
489 secNumberInit = fvect->size();
491 secNumberFinal = fvect->size();
498 secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef,k,Shell);
502 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(PartDef,k,Shell);
505 if (verboseLevel > 3)
508 G4cout <<
"Shell: " << Shell <<
" Kin. energy (eV)=" << k/eV
509 <<
" Sec. energy (eV)=" << secondaryKinetic/eV <<
G4endl;
519 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
520 G4double finalPx = totalMomentum*primaryDirection.
x() - deltaTotalMomentum*deltaDirection.
x();
521 G4double finalPy = totalMomentum*primaryDirection.
y() - deltaTotalMomentum*deltaDirection.
y();
522 G4double finalPz = totalMomentum*primaryDirection.
z() - deltaTotalMomentum*deltaDirection.
z();
523 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
524 finalPx /= finalMomentum;
525 finalPy /= finalMomentum;
526 finalPz /= finalMomentum;
529 direction.
set(finalPx,finalPy,finalPz);
538 for (std::size_t j=secNumberInit; j < secNumberFinal; ++j) {
539 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
544 if (secondaryKinetic>0)
547 fvect->push_back(dp);
560 if ((k+SiStructure.
Energy(shell))/2. > k) maximumEnergyTransfer=k;
561 else maximumEnergyTransfer = (k+SiStructure.
Energy(shell))/2.;
566 G4double maxEnergy = maximumEnergyTransfer;
567 G4int nEnergySteps = 100;
570 G4double stpEnergy(std::pow(maxEnergy/value, 1./
static_cast<G4double>(nEnergySteps-1)));
571 G4int step(nEnergySteps);
576 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
580 G4double secondaryElectronKineticEnergy=0.;
583 secondaryElectronKineticEnergy =
G4UniformRand() * (maximumEnergyTransfer-SiStructure.
Energy(shell));
587 return secondaryElectronKineticEnergy;
592 G4double maximumEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
596 G4double maxEnergy = maximumEnergyTransfer;
597 G4int nEnergySteps = 100;
600 G4double stpEnergy(std::pow(maxEnergy/value, 1./
static_cast<G4double>(nEnergySteps-1)));
601 G4int step(nEnergySteps);
606 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
610 G4double secondaryElectronKineticEnergy = 0.;
613 secondaryElectronKineticEnergy =
G4UniformRand() * (maximumEnergyTransfer-SiStructure.
Energy(shell));
617 return secondaryElectronKineticEnergy;
666 if (energyTransfer >= SiStructure.
Energy(LevelIndex))
683 auto t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
686 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
688 auto e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
690 auto e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
700 xs11 = eDiffCrossSectionData[LevelIndex][valueT1][valueE11];
701 xs12 = eDiffCrossSectionData[LevelIndex][valueT1][valueE12];
702 xs21 = eDiffCrossSectionData[LevelIndex][valueT2][valueE21];
703 xs22 = eDiffCrossSectionData[LevelIndex][valueT2][valueE22];
710 auto t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
712 if (energyTransfer <= pVecm[(*t1)].back() && energyTransfer <= pVecm[(*t2)].back() )
714 auto e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
717 auto e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
727 xs11 = pDiffCrossSectionData[LevelIndex][valueT1][valueE11];
728 xs12 = pDiffCrossSectionData[LevelIndex][valueT1][valueE12];
729 xs21 = pDiffCrossSectionData[LevelIndex][valueT2][valueE21];
730 xs22 = pDiffCrossSectionData[LevelIndex][valueT2][valueE22];
734 sigma = QuadInterpolator( valueE11, valueE12,
755 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
758 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
759 G4double b = std::log10(xs2) - a*std::log10(e2);
760 G4double sigma = a*std::log10(e) + b;
761 value = (std::pow(10.,sigma));
766 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
770 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
775 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0))
779 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
794 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
795 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
796 G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
806 auto pos = tableData.find(particle);
807 if (pos != tableData.cend())
811 if (table !=
nullptr)
822 value += valuesBuffer[i];
833 if (valuesBuffer[i] > value)
835 delete[] valuesBuffer;
838 value -= valuesBuffer[i];
841 if (valuesBuffer)
delete[] valuesBuffer;
860 G4double secondaryElectronKineticEnergy = 0.;
866 - SiStructure.
Energy(shell);
868 if (secondaryElectronKineticEnergy < 0.)
871 return secondaryElectronKineticEnergy;
878 G4int ionizationLevelIndex,
895 G4double maximumEnergyTransfer1 = 0;
896 G4double maximumEnergyTransfer2 = 0;
897 G4double maximumEnergyTransferP = 4.* (electron_mass_c2 / proton_mass_c2) * k;
898 G4double bindingEnergy = SiStructure.
Energy(ionizationLevelIndex)*1e6;
903 auto k2 = std::upper_bound(eTdummyVec.begin(),
920 if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
921 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
924 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
925 eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
927 auto prob11 = prob12 - 1;
929 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
930 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
932 auto prob21 = prob22 - 1;
936 valuePROB21 = *prob21;
937 valuePROB22 = *prob22;
938 valuePROB12 = *prob12;
939 valuePROB11 = *prob11;
942 if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
943 else nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
946 if ((valueK1+bindingEnergy)/2. > valueK1) maximumEnergyTransfer1=valueK1;
947 else maximumEnergyTransfer1 = (valueK1+bindingEnergy)/2.;
949 nrjTransf12 = maximumEnergyTransfer1;
951 else nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
953 if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
954 else nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
957 if ((valueK2+bindingEnergy)/2. > valueK2) maximumEnergyTransfer2=valueK2;
958 else maximumEnergyTransfer2 = (valueK2+bindingEnergy)/2.;
960 nrjTransf22 = maximumEnergyTransfer2;
963 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
966 if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
969 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
970 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
973 auto prob21 = prob22 - 1;
977 valuePROB21 = *prob21;
978 valuePROB22 = *prob22;
980 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
981 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
983 G4double interpolatedvalue2 = Interpolate(valuePROB21,
990 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
998 auto k2 = std::upper_bound(pTdummyVec.begin(),
1017 if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1018 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1021 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1022 pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1025 auto prob11 = prob12 - 1;
1028 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1029 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1032 auto prob21 = prob22 - 1;
1036 valuePROB21 = *prob21;
1037 valuePROB22 = *prob22;
1038 valuePROB12 = *prob12;
1039 valuePROB11 = *prob11;
1042 if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
1043 else nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1044 if(valuePROB12 == 1) nrjTransf12 = maximumEnergyTransferP;
1045 else nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1046 if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
1047 else nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1048 if(valuePROB22 == 1) nrjTransf22 = maximumEnergyTransferP;
1049 else nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1054 if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1057 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1058 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1061 auto prob21 = prob22 - 1;
1065 valuePROB21 = *prob21;
1066 valuePROB22 = *prob22;
1068 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1069 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1071 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1078 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1085 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1088 if (nrjTransfProduct != 0.)
1090 nrj = QuadInterpolator(valuePROB11,
const char * G4FindDataDir(const char *)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4GLOB_DLL std::ostream G4cout
void set(double x, double y, double z)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * ElectronDefinition()
static G4Electron * Electron()
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
const G4Material * GetBaseMaterial() const
G4double GetTotNbOfAtomsPerVolume() const
G4bool LoadData(const G4String &argFileName) override
G4double FindValue(G4double e, G4int componentId=0) const override
size_t NumberOfComponents(void) const override
const G4VEMDataSet * GetComponent(G4int componentId) const override
G4double TransferedEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4ParticleChangeForGamma * fParticleChangeForGamma
G4MicroElecInelasticModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="MicroElecInelasticModel")
virtual ~G4MicroElecInelasticModel()
G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax) override
G4double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
G4double Energy(G4int level)
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
static G4NistManager * Instance()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4double GetPDGMass() const
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
void SetHighEnergyLimit(G4double)
G4VEmAngularDistribution * GetAngularDistribution()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void SetLowEnergyLimit(G4double)
void SetDeexcitationFlag(G4bool val)
void SetAngularDistribution(G4VEmAngularDistribution *)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)