78 G4cout <<
"MicroElec inelastic model is constructed " <<
G4endl;
99 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
120 G4cout <<
"Calling G4MicroElecInelasticModel::Initialise()" <<
G4endl;
124 G4String fileElectron(
"microelec/sigma_inelastic_e_Si");
125 G4String fileProton(
"microelec/sigma_inelastic_p_Si");
135 char *path = getenv(
"G4LEDATA");
154 std::ostringstream eFullFileName;
156 if (
fasterCode) eFullFileName << path <<
"/microelec/sigmadiff_cumulated_inelastic_e_Si.dat";
157 else eFullFileName << path <<
"/microelec/sigmadiff_inelastic_e_Si.dat";
159 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
161 if (!eDiffCrossSection)
164 FatalException,
"Missing data file:/microelec/sigmadiff_cumulated_inelastic_e_Si.dat");
166 else G4Exception(
"G4MicroElecInelasticModel::Initialise",
"em0003",
167 FatalException,
"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
182 for (
int j=0; j<6; j++)
198 while(!eDiffCrossSection.eof())
202 eDiffCrossSection>>tDummy>>eDummy;
206 for (
int j=0; j<6; j++)
208 eDiffCrossSection>>
tmp;
215 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
221 eVecm[tDummy].push_back(eDummy);
246 std::ostringstream pFullFileName;
248 if (
fasterCode) pFullFileName << path <<
"/microelec/sigmadiff_cumulated_inelastic_p_Si.dat";
249 else pFullFileName << path <<
"/microelec/sigmadiff_inelastic_p_Si.dat";
251 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
253 if (!pDiffCrossSection)
256 FatalException,
"Missing data file:/microelec/sigmadiff_cumulated_inelastic_p_Si.dat");
258 else G4Exception(
"G4MicroElecInelasticModel::Initialise",
"em0003",
259 FatalException,
"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
263 while(!pDiffCrossSection.eof())
267 pDiffCrossSection>>tDummy>>eDummy;
269 for (
int j=0; j<6; j++)
275 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
276 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
281 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
282 pVecm[tDummy].push_back(eDummy);
288 if (particle==electronDef)
294 if (particle==protonDef)
302 G4cout <<
"MicroElec Inelastic model is initialized " << G4endl
331 G4cout <<
"Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" <<
G4endl;
364 G4cout <<
"Before scaling : " << G4endl
365 <<
"Particle : " << nameLocal <<
", mass : " << Mion_c2/
proton_mass_c2 <<
"*mp, charge " << Zeff
366 <<
", Ekin (eV) = " << ekin/
eV <<
G4endl ;
369 nameLocal =
"proton" ;
372 G4cout <<
"After scaling : " << G4endl
373 <<
"Particle : " << nameLocal <<
", Ekin (eV) = " << ekin/
eV <<
G4endl ;
379 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
383 lowLim = pos1->second;
386 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
390 highLim = pos2->second;
393 if (ekin >= lowLim && ekin < highLim)
395 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
408 G4Exception(
"G4MicroElecInelasticModel::CrossSectionPerVolume",
"em0002",
FatalException,
"Model not applicable to particle type.");
423 G4cout <<
" - Cross section per Si atom (cm^2)=" << sigma*Zeff2/
cm2 <<
G4endl;
424 G4cout <<
" - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./
cm) << G4endl;
428 return sigma*density*Zeff2;
443 G4cout <<
"Calling SampleSecondaries() of G4MicroElecInelasticModel" <<
G4endl;
453 G4String nameLocal2 = particleName ;
460 nameLocal2 =
"proton" ;
463 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
468 lowLim = pos1->second;
471 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
476 highLim = pos2->second;
479 if (k >= lowLim && k < highLim)
482 G4double totalEnergy = ekin + particleMass;
483 G4double pSquare = ekin * (totalEnergy + particleMass);
484 G4double totalMomentum = std::sqrt(pSquare);
507 G4cout <<
"Shell: " << Shell <<
", energy: " << bindingEnergy/
eV <<
G4endl;
512 G4int secNumberInit = 0;
513 G4int secNumberFinal = 0;
516 if (k<bindingEnergy)
return;
534 secNumberInit = fvect->size();
536 secNumberFinal = fvect->size();
554 G4cout <<
"Shell: " << Shell <<
" Kin. energy (eV)=" << k/
eV
555 <<
" Sec. energy (eV)=" << secondaryKinetic/
eV <<
G4endl;
567 G4double finalPx = totalMomentum*primaryDirection.
x() - deltaTotalMomentum*deltaDirection.
x();
568 G4double finalPy = totalMomentum*primaryDirection.
y() - deltaTotalMomentum*deltaDirection.
y();
569 G4double finalPz = totalMomentum*primaryDirection.
z() - deltaTotalMomentum*deltaDirection.
z();
570 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
571 finalPx /= finalMomentum;
572 finalPy /= finalMomentum;
573 finalPz /= finalMomentum;
576 direction.
set(finalPx,finalPy,finalPz);
584 for (
G4int j=secNumberInit; j < secNumberFinal; j++) {
585 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
590 if (secondaryKinetic>0)
593 fvect->push_back(dp);
613 G4double maxEnergy = maximumEnergyTransfer;
614 G4int nEnergySteps = 100;
617 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
618 G4int step(nEnergySteps);
623 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
628 G4double secondaryElectronKineticEnergy=0.;
635 return secondaryElectronKineticEnergy;
645 G4double maxEnergy = maximumEnergyTransfer;
646 G4int nEnergySteps = 100;
649 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
650 G4int step(nEnergySteps);
655 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
659 G4double secondaryElectronKineticEnergy = 0.;
666 return secondaryElectronKineticEnergy;
734 std::vector<double>::iterator
t1 = t2-1;
736 if (energyTransfer <=
eVecm[(*t1)].back() && energyTransfer <=
eVecm[(*t2)].back() )
738 std::vector<double>::iterator e12 = std::upper_bound(
eVecm[(*t1)].begin(),
eVecm[(*t1)].end(), energyTransfer);
739 std::vector<double>::iterator e11 = e12-1;
741 std::vector<double>::iterator e22 = std::upper_bound(
eVecm[(*t2)].begin(),
eVecm[(*t2)].end(), energyTransfer);
742 std::vector<double>::iterator e21 = e22-1;
763 std::vector<double>::iterator
t1 = t2-1;
764 if (energyTransfer <=
pVecm[(*t1)].back() && energyTransfer <=
pVecm[(*t2)].back() )
766 std::vector<double>::iterator e12 = std::upper_bound(
pVecm[(*t1)].begin(),
pVecm[(*t1)].end(), energyTransfer);
767 std::vector<double>::iterator e11 = e12-1;
769 std::vector<double>::iterator e22 = std::upper_bound(
pVecm[(*t2)].begin(),
pVecm[(*t2)].end(), energyTransfer);
770 std::vector<double>::iterator e21 = e22-1;
813 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
816 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
817 G4double b = std::log10(xs2) - a*std::log10(e2);
818 G4double sigma = a*std::log10(e) + b;
819 value = (std::pow(10.,sigma));
824 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 &&
fasterCode)
828 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
834 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0))
838 value = (d1 + (d2 -
d1) * (e - e1) / (e2 - e1));
866 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator
pos;
884 value += valuesBuffer[i];
895 if (valuesBuffer[i] > value)
897 delete[] valuesBuffer;
900 value -= valuesBuffer[i];
903 if (valuesBuffer)
delete[] valuesBuffer;
922 G4double secondaryElectronKineticEnergy = 0.;
932 if (secondaryElectronKineticEnergy < 0.)
936 return secondaryElectronKineticEnergy;
943 G4int ionizationLevelIndex,
960 G4double maximumEnergyTransfer1 = 0;
961 G4double maximumEnergyTransfer2 = 0;
968 std::vector<double>::iterator k2 = std::upper_bound(
eTdummyVec.begin(),
971 std::vector<double>::iterator k1 = k2 - 1;
988 std::vector<double>::iterator prob12 =
989 std::upper_bound(
eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
993 std::vector<double>::iterator prob11 = prob12 - 1;
995 std::vector<double>::iterator prob22 =
996 std::upper_bound(
eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1000 std::vector<double>::iterator prob21 = prob22 - 1;
1004 valuePROB21 = *prob21;
1005 valuePROB22 = *prob22;
1006 valuePROB12 = *prob12;
1007 valuePROB11 = *prob11;
1016 else nrjTransf11 =
eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1017 if(valuePROB12 == 1)
1019 if ((valueK1+bindingEnergy)/2. > valueK1) maximumEnergyTransfer1=valueK1;
1022 nrjTransf12 = maximumEnergyTransfer1;
1024 else nrjTransf12 =
eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1027 else nrjTransf21 =
eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1028 if(valuePROB22 == 1)
1030 if ((valueK2+bindingEnergy)/2. > valueK2) maximumEnergyTransfer2=valueK2;
1033 nrjTransf22 = maximumEnergyTransfer2;
1035 else nrjTransf22 =
eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1055 std::vector<double>::iterator prob22 =
1056 std::upper_bound(
eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1060 std::vector<double>::iterator prob21 = prob22 - 1;
1064 valuePROB21 = *prob21;
1065 valuePROB22 = *prob22;
1069 nrjTransf21 =
eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1070 nrjTransf22 =
eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1100 std::vector<double>::iterator k2 = std::upper_bound(
pTdummyVec.begin(),
1104 std::vector<double>::iterator k1 = k2 - 1;
1122 std::vector<double>::iterator prob12 =
1123 std::upper_bound(
pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1127 std::vector<double>::iterator prob11 = prob12 - 1;
1129 std::vector<double>::iterator prob22 =
1130 std::upper_bound(
pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1134 std::vector<double>::iterator prob21 = prob22 - 1;
1138 valuePROB21 = *prob21;
1139 valuePROB22 = *prob22;
1140 valuePROB12 = *prob12;
1141 valuePROB11 = *prob11;
1150 else nrjTransf11 =
pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1151 if(valuePROB12 == 1) nrjTransf12 = maximumEnergyTransferP;
1152 else nrjTransf12 =
pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1154 else nrjTransf21 =
pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1155 if(valuePROB22 == 1) nrjTransf22 = maximumEnergyTransferP;
1156 else nrjTransf22 =
pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1177 std::vector<double>::iterator prob22 =
1178 std::upper_bound(
pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1182 std::vector<double>::iterator prob21 = prob22 - 1;
1186 valuePROB21 = *prob21;
1187 valuePROB22 = *prob22;
1191 nrjTransf21 =
pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1192 nrjTransf22 =
pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1219 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1223 if (nrjTransfProduct != 0.)
void set(double x, double y, double z)
TriDimensionMap pNrjTransfData[7]
G4int RandomSelect(G4double energy, const G4String &particle)
std::vector< ExP01TrackerHit * > a
static const G4double pos
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4VAtomDeexcitation * AtomDeexcitation()
static constexpr double MeV
void SetHighEnergyLimit(G4double)
G4double RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
static constexpr double keV
const G4ThreeVector & GetMomentumDirection() const
const G4String & GetParticleName() const
G4double GetPDGCharge() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual G4bool LoadData(const G4String &argFileName)
G4MicroElecSiStructure SiStructure
G4double GetPDGMass() const
G4double Interpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4double LowEnergyLimit() const
static constexpr double proton_mass_c2
G4VAtomDeexcitation * fAtomDeexcitation
void SetDeexcitationFlag(G4bool val)
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4ParticleDefinition * GetDefinition() const
static constexpr double cm2
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
const XML_Char int const XML_Char * value
static G4Proton * ProtonDefinition()
virtual size_t NumberOfComponents(void) const
static constexpr double electron_mass_c2
virtual ~G4MicroElecInelasticModel()
G4ParticleChangeForGamma * fParticleChangeForGamma
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
G4double Energy(G4int level)
static constexpr double eV
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
void SetAngularDistribution(G4VEmAngularDistribution *)
static G4Electron * Electron()
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
static G4Electron * ElectronDefinition()
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
void SetLowEnergyLimit(G4double)
std::vector< double > pTdummyVec
G4MicroElecInelasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MicroElecInelasticModel")
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
G4double GetKineticEnergy() const
static constexpr double cm
G4GLOB_DLL std::ostream G4cout
TriDimensionMap eDiffCrossSectionData[7]
static G4LossTableManager * Instance()
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
G4VEmAngularDistribution * GetAngularDistribution()
const G4Material * GetMaterial() const
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
TriDimensionMap pDiffCrossSectionData[7]
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double GetTotNbOfAtomsPerVolume() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
TriDimensionMap eNrjTransfData[7]
G4double bindingEnergy(G4int A, G4int Z)
G4double HighEnergyLimit() const
static constexpr double GeV
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4Material * GetBaseMaterial() const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
std::vector< double > eTdummyVec
G4double TransferedEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
static G4NistManager * Instance()
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit