Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Protected Attributes | List of all members
G4MicroElecInelasticModel Class Reference

#include <G4MicroElecInelasticModel.hh>

+ Inheritance diagram for G4MicroElecInelasticModel:

Public Member Functions

 G4MicroElecInelasticModel (const G4ParticleDefinition *p=0, const G4String &nam="MicroElecInelasticModel")
 
virtual ~G4MicroElecInelasticModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual G4double CrossSectionPerVolume (const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
double DifferentialCrossSection (G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
 
G4double TransferedEnergy (G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
 
void SelectFasterComputation (G4bool input)
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)=0
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
virtual G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectTargetAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
G4int SelectIsotopeNumber (const G4Element *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
G4VEmModel * GetTripletModel ()
 
void SetTripletModel (G4VEmModel *)
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool LPMFlag () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy) const
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
void SetUseBaseMaterials (G4bool val)
 
G4bool UseBaseMaterials () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4Element * GetCurrentElement () const
 
const G4Isotope * GetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 
G4VEmModel & operator= (const G4VEmModel &right)=delete
 
 G4VEmModel (const G4VEmModel &)=delete
 

Protected Attributes

G4ParticleChangeForGamma * fParticleChangeForGamma
 
- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const G4Material * pBaseMaterial
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 
G4double inveplus
 
G4double pFactor
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCouple * CurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 

Detailed Description

Definition at line 64 of file G4MicroElecInelasticModel.hh.

Constructor & Destructor Documentation

◆ G4MicroElecInelasticModel()

G4MicroElecInelasticModel::G4MicroElecInelasticModel ( const G4ParticleDefinition *  p = 0,
const G4String &  nam = "MicroElecInelasticModel" 
)

Definition at line 62 of file G4MicroElecInelasticModel.cc.

64:G4VEmModel(nam),isInitialised(false)
65{
66 nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
67
68 verboseLevel= 0;
69 // Verbosity scale:
70 // 0 = nothing
71 // 1 = warning for energy non-conservation
72 // 2 = details of energy budget
73 // 3 = calculation of cross sections, file openings, sampling of atoms
74 // 4 = entering in methods
75
76 if( verboseLevel>0 )
77 {
78 G4cout << "MicroElec inelastic model is constructed " << G4endl;
79 }
80
81 //Mark this model as "applicable" for atomic deexcitation
82 SetDeexcitationFlag(true);
83 fAtomDeexcitation = 0;
84 fParticleChangeForGamma = 0;
85
86 // default generator
87 SetAngularDistribution(new G4DeltaAngle());
88
89 // Selection of computation method
90 fasterCode = true; //false;
91}
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4MicroElecInelasticModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4MicroElecInelasticModel.hh:97
G4NistManager::FindOrBuildMaterial
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
Definition: G4NistManager.hh:491
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:70
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:813
G4VEmModel::SetAngularDistribution
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:618

◆ ~G4MicroElecInelasticModel()

G4MicroElecInelasticModel::~G4MicroElecInelasticModel ( )
virtual

Definition at line 95 of file G4MicroElecInelasticModel.cc.

96{
97 // Cross section
98
99 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
100 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
101 {
102 G4MicroElecCrossSectionDataSet* table = pos->second;
103 delete table;
104 }
105
106 // Final state
107
108 eVecm.clear();
109 pVecm.clear();
110
111}

Member Function Documentation

◆ CrossSectionPerVolume()

G4double G4MicroElecInelasticModel::CrossSectionPerVolume ( const G4Material *  material,
const G4ParticleDefinition *  p,
G4double  ekin,
G4double  emin,
G4double  emax 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 324 of file G4MicroElecInelasticModel.cc.

329{
330 if (verboseLevel > 3)
331 G4cout << "Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" << G4endl;
332
333 G4double density = material->GetTotNbOfAtomsPerVolume();
334
335 /* if (
336 particleDefinition != G4Proton::ProtonDefinition()
337 &&
338 particleDefinition != G4Electron::ElectronDefinition()
339 &&
340 particleDefinition != G4GenericIon::GenericIonDefinition()
341 )
342
343 return 0;*/
344
345 // Calculate total cross section for model
346
347 G4double lowLim = 0;
348 G4double highLim = 0;
349 G4double sigma=0;
350
351 const G4String& particleName = particleDefinition->GetParticleName();
352 G4String nameLocal = particleName ;
353
354 G4double Zeff2 = 1.0;
355 G4double Mion_c2 = particleDefinition->GetPDGMass();
356
357 if (Mion_c2 > proton_mass_c2)
358 {
359 G4ionEffectiveCharge EffCharge ;
360 G4double Zeff = EffCharge.EffectiveCharge(particleDefinition, material,ekin);
361 Zeff2 = Zeff*Zeff;
362
363 if (verboseLevel > 3)
364 G4cout << "Before scaling : " << G4endl
365 << "Particle : " << nameLocal << ", mass : " << Mion_c2/proton_mass_c2 << "*mp, charge " << Zeff
366 << ", Ekin (eV) = " << ekin/eV << G4endl ;
367
368 ekin *= proton_mass_c2/Mion_c2 ;
369 nameLocal = "proton" ;
370
371 if (verboseLevel > 3)
372 G4cout << "After scaling : " << G4endl
373 << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl ;
374 }
375
376 if (material == nistSi || material->GetBaseMaterial() == nistSi)
377 {
378
379 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
380 pos1 = lowEnergyLimit.find(nameLocal);
381 if (pos1 != lowEnergyLimit.end())
382 {
383 lowLim = pos1->second;
384 }
385
386 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
387 pos2 = highEnergyLimit.find(nameLocal);
388 if (pos2 != highEnergyLimit.end())
389 {
390 highLim = pos2->second;
391 }
392
393 if (ekin >= lowLim && ekin < highLim)
394 {
395 std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
396 pos = tableData.find(nameLocal);
397
398 if (pos != tableData.end())
399 {
400 G4MicroElecCrossSectionDataSet* table = pos->second;
401 if (table != 0)
402 {
403 sigma = table->FindValue(ekin);
404 }
405 }
406 else
407 {
408 G4Exception("G4MicroElecInelasticModel::CrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
409 }
410 }
411 else
412 {
413 if (nameLocal!="e-")
414 {
415 // G4cout << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl;
416 // G4cout << "### Warning: particle energy out of bounds! ###" << G4endl;
417 }
418 }
419
420 if (verboseLevel > 3)
421 {
422 G4cout << "---> Kinetic energy (eV)=" << ekin/eV << G4endl;
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;
425 }
426
427 } // if (SiMaterial)
428 return sigma*density*Zeff2;
429
430
431}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4double
double G4double
Definition: G4Types.hh:83
G4Material::GetBaseMaterial
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:231
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
G4MicroElecCrossSectionDataSet::FindValue
virtual G4double FindValue(G4double e, G4int componentId=0) const
Definition: G4MicroElecCrossSectionDataSet.cc:458
G4ionEffectiveCharge::EffectiveCharge
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
Definition: G4ionEffectiveCharge.cc:87

◆ DifferentialCrossSection()

double G4MicroElecInelasticModel::DifferentialCrossSection ( G4ParticleDefinition *  aParticleDefinition,
G4double  k,
G4double  energyTransfer,
G4int  shell 
)

Definition at line 708 of file G4MicroElecInelasticModel.cc.

712{
713 G4double sigma = 0.;
714
715 if (energyTransfer >= SiStructure.Energy(LevelIndex))
716 {
717 G4double valueT1 = 0;
718 G4double valueT2 = 0;
719 G4double valueE21 = 0;
720 G4double valueE22 = 0;
721 G4double valueE12 = 0;
722 G4double valueE11 = 0;
723
724 G4double xs11 = 0;
725 G4double xs12 = 0;
726 G4double xs21 = 0;
727 G4double xs22 = 0;
728
729 if (particleDefinition == G4Electron::ElectronDefinition())
730 {
731 // k should be in eV and energy transfer eV also
732
733 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
734 std::vector<double>::iterator t1 = t2-1;
735 // SI : the following condition avoids situations where energyTransfer >last vector element
736 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
737 {
738 std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
739 std::vector<double>::iterator e11 = e12-1;
740
741 std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
742 std::vector<double>::iterator e21 = e22-1;
743
744 valueT1 =*t1;
745 valueT2 =*t2;
746 valueE21 =*e21;
747 valueE22 =*e22;
748 valueE12 =*e12;
749 valueE11 =*e11;
750
751 xs11 = eDiffCrossSectionData[LevelIndex][valueT1][valueE11];
752 xs12 = eDiffCrossSectionData[LevelIndex][valueT1][valueE12];
753 xs21 = eDiffCrossSectionData[LevelIndex][valueT2][valueE21];
754 xs22 = eDiffCrossSectionData[LevelIndex][valueT2][valueE22];
755 }
756
757 }
758
759 if (particleDefinition == G4Proton::ProtonDefinition())
760 {
761 // k should be in eV and energy transfer eV also
762 std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
763 std::vector<double>::iterator t1 = t2-1;
764 if (energyTransfer <= pVecm[(*t1)].back() && energyTransfer <= pVecm[(*t2)].back() )
765 {
766 std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
767 std::vector<double>::iterator e11 = e12-1;
768
769 std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
770 std::vector<double>::iterator e21 = e22-1;
771
772 valueT1 =*t1;
773 valueT2 =*t2;
774 valueE21 =*e21;
775 valueE22 =*e22;
776 valueE12 =*e12;
777 valueE11 =*e11;
778
779 xs11 = pDiffCrossSectionData[LevelIndex][valueT1][valueE11];
780 xs12 = pDiffCrossSectionData[LevelIndex][valueT1][valueE12];
781 xs21 = pDiffCrossSectionData[LevelIndex][valueT2][valueE21];
782 xs22 = pDiffCrossSectionData[LevelIndex][valueT2][valueE22];
783 }
784 }
785
786 // G4double xsProduct = xs11 * xs12 * xs21 * xs22;
787 // if (xsProduct != 0.)
788 // {
789 sigma = QuadInterpolator( valueE11, valueE12,
790 valueE21, valueE22,
791 xs11, xs12,
792 xs21, xs22,
793 valueT1, valueT2,
794 k, energyTransfer);
795 // }
796
797 }
798
799 return sigma;
800}
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:88
G4MicroElecSiStructure::Energy
G4double Energy(G4int level)
Definition: G4MicroElecSiStructure.cc:63
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:87

◆ Initialise()

void G4MicroElecInelasticModel::Initialise ( const G4ParticleDefinition *  particle,
const G4DataVector &   
)
virtual

Implements G4VEmModel.

Definition at line 115 of file G4MicroElecInelasticModel.cc.

117{
118
119 if (verboseLevel > 3)
120 G4cout << "Calling G4MicroElecInelasticModel::Initialise()" << G4endl;
121
122 // Energy limits
123
124 G4String fileElectron("microelec/sigma_inelastic_e_Si");
125 G4String fileProton("microelec/sigma_inelastic_p_Si");
126
127 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
128 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
129
130 G4String electron;
131 G4String proton;
132
133 G4double scaleFactor = 1e-18 * cm *cm;
134
135 char *path = std::getenv("G4LEDATA");
136
137 // *** ELECTRON
138 electron = electronDef->GetParticleName();
139
140 tableFile[electron] = fileElectron;
141
142 lowEnergyLimit[electron] = 16.7 * eV;
143 highEnergyLimit[electron] = 100.0 * MeV;
144
145 // Cross section
146
147 G4MicroElecCrossSectionDataSet* tableE = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
148 tableE->LoadData(fileElectron);
149
150 tableData[electron] = tableE;
151
152 // Final state
153
154 std::ostringstream eFullFileName;
155
156 if (fasterCode) eFullFileName << path << "/microelec/sigmadiff_cumulated_inelastic_e_Si.dat";
157 else eFullFileName << path << "/microelec/sigmadiff_inelastic_e_Si.dat";
158
159 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
160
161 if (!eDiffCrossSection)
162 {
163 if (fasterCode) G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
164 FatalException,"Missing data file:/microelec/sigmadiff_cumulated_inelastic_e_Si.dat");
165
166 else G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
167 FatalException,"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
168
169 }
170
171 //
172
173 // Clear the arrays for re-initialization case (MT mode)
174 // Octobre 22nd, 2014 - Melanie Raine
175
176 eTdummyVec.clear();
177 pTdummyVec.clear();
178
179 eVecm.clear();
180 pVecm.clear();
181
182 for (int j=0; j<6; j++)
183 {
184 eProbaShellMap[j].clear();
185 pProbaShellMap[j].clear();
186
187 eDiffCrossSectionData[j].clear();
188 pDiffCrossSectionData[j].clear();
189
190 eNrjTransfData[j].clear();
191 pNrjTransfData[j].clear();
192 }
193
194 //
195
196
197 eTdummyVec.push_back(0.);
198 while(!eDiffCrossSection.eof())
199 {
200 double tDummy;
201 double eDummy;
202 eDiffCrossSection>>tDummy>>eDummy;
203 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
204
205 double tmp;
206 for (int j=0; j<6; j++)
207 {
208 eDiffCrossSection>> tmp;
209
210 eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
211
212 if (fasterCode)
213 {
214 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
215 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
216 }
217 else
218 {
219 // SI - only if eof is not reached !
220 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
221 eVecm[tDummy].push_back(eDummy);
222 }
223
224 }
225 }
226 //
227
228 // *** PROTON
229
230 proton = protonDef->GetParticleName();
231
232 tableFile[proton] = fileProton;
233
234 lowEnergyLimit[proton] = 50. * keV;
235 highEnergyLimit[proton] = 10. * GeV;
236
237 // Cross section
238
239 G4MicroElecCrossSectionDataSet* tableP = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
240 tableP->LoadData(fileProton);
241
242 tableData[proton] = tableP;
243
244 // Final state
245
246 std::ostringstream pFullFileName;
247
248 if (fasterCode) pFullFileName << path << "/microelec/sigmadiff_cumulated_inelastic_p_Si.dat";
249 else pFullFileName << path << "/microelec/sigmadiff_inelastic_p_Si.dat";
250
251 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
252
253 if (!pDiffCrossSection)
254 {
255 if (fasterCode) G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
256 FatalException,"Missing data file:/microelec/sigmadiff_cumulated_inelastic_p_Si.dat");
257
258 else G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
259 FatalException,"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
260 }
261
262 pTdummyVec.push_back(0.);
263 while(!pDiffCrossSection.eof())
264 {
265 double tDummy;
266 double eDummy;
267 pDiffCrossSection>>tDummy>>eDummy;
268 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
269 for (int j=0; j<6; j++)
270 {
271 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
272
273 if (fasterCode)
274 {
275 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
276 pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
277 }
278 else
279 {
280 // SI - only if eof is not reached !
281 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
282 pVecm[tDummy].push_back(eDummy);
283 }
284 }
285 }
286
287
288 if (particle==electronDef)
289 {
290 SetLowEnergyLimit(lowEnergyLimit[electron]);
291 SetHighEnergyLimit(highEnergyLimit[electron]);
292 }
293
294 if (particle==protonDef)
295 {
296 SetLowEnergyLimit(lowEnergyLimit[proton]);
297 SetHighEnergyLimit(highEnergyLimit[proton]);
298 }
299
300 if( verboseLevel>0 )
301 {
302 G4cout << "MicroElec Inelastic model is initialized " << G4endl
303 << "Energy range: "
304 << LowEnergyLimit() / keV << " keV - "
305 << HighEnergyLimit() / MeV << " MeV for "
306 << particle->GetParticleName()
307 << " with mass (amu) " << particle->GetPDGMass()/proton_mass_c2
308 << " and charge " << particle->GetPDGCharge()
309 << G4endl << G4endl ;
310 }
311
312 //
313
314 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
315
316 if (isInitialised) { return; }
317 fParticleChangeForGamma = GetParticleChangeForGamma();
318 isInitialised = true;
319
320}
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:87
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:404
G4MicroElecCrossSectionDataSet::LoadData
virtual G4bool LoadData(const G4String &argFileName)
Definition: G4MicroElecCrossSectionDataSet.cc:58
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:112
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:108
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:757
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:133
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:652
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:645
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:764

◆ SampleSecondaries()

void G4MicroElecInelasticModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  fvect,
const G4MaterialCutsCouple *  couple,
const G4DynamicParticle *  particle,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 435 of file G4MicroElecInelasticModel.cc.

440{
441
442 if (verboseLevel > 3)
443 G4cout << "Calling SampleSecondaries() of G4MicroElecInelasticModel" << G4endl;
444
445 G4double lowLim = 0;
446 G4double highLim = 0;
447
448 G4double ekin = particle->GetKineticEnergy();
449 G4double k = ekin ;
450
451 G4ParticleDefinition* PartDef = particle->GetDefinition();
452 const G4String& particleName = PartDef->GetParticleName();
453 G4String nameLocal2 = particleName ;
454 G4double particleMass = particle->GetDefinition()->GetPDGMass();
455
456 if (particleMass > proton_mass_c2)
457 {
458 k *= proton_mass_c2/particleMass ;
459 PartDef = G4Proton::ProtonDefinition();
460 nameLocal2 = "proton" ;
461 }
462
463 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
464 pos1 = lowEnergyLimit.find(nameLocal2);
465
466 if (pos1 != lowEnergyLimit.end())
467 {
468 lowLim = pos1->second;
469 }
470
471 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
472 pos2 = highEnergyLimit.find(nameLocal2);
473
474 if (pos2 != highEnergyLimit.end())
475 {
476 highLim = pos2->second;
477 }
478
479 if (k >= lowLim && k < highLim)
480 {
481 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
482 G4double totalEnergy = ekin + particleMass;
483 G4double pSquare = ekin * (totalEnergy + particleMass);
484 G4double totalMomentum = std::sqrt(pSquare);
485
486 G4int Shell = 0;
487
488 /* if (!fasterCode)*/ Shell = RandomSelect(k,nameLocal2);
489
490 // SI: The following protection is necessary to avoid infinite loops :
491 // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
492 // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
493 // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
494 // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
495
496 /*if (fasterCode)
497 do
498 {
499 Shell = RandomSelect(k,nameLocal2);
500 }while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());*/
501
502 G4double bindingEnergy = SiStructure.Energy(Shell);
503
504 if (verboseLevel > 3)
505 {
506 G4cout << "---> Kinetic energy (eV)=" << k/eV << G4endl ;
507 G4cout << "Shell: " << Shell << ", energy: " << bindingEnergy/eV << G4endl;
508 }
509
510 // sample deexcitation
511
512 G4int secNumberInit = 0; // need to know at a certain point the energy of secondaries
513 G4int secNumberFinal = 0; // So I'll make the difference and then sum the energies
514
515 //SI: additional protection if tcs interpolation method is modified
516 if (k<bindingEnergy) return;
517
518 G4int Z = 14;
519
520 if(fAtomDeexcitation && Shell > 2) {
521
522 G4AtomicShellEnumerator as = fKShell;
523
524 if (Shell == 4)
525 {
526 as = G4AtomicShellEnumerator(1);
527 }
528 else if (Shell == 3)
529 {
530 as = G4AtomicShellEnumerator(3);
531 }
532
533 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
534 secNumberInit = fvect->size();
535 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
536 secNumberFinal = fvect->size();
537 }
538
539 G4double secondaryKinetic=-1000*eV;
540
541 if (!fasterCode)
542 {
543 secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef,k,Shell);
544 }
545 else
546 {
547 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(PartDef,k,Shell);
548 }
549
550
551 if (verboseLevel > 3)
552 {
553 G4cout << "Ionisation process" << G4endl;
554 G4cout << "Shell: " << Shell << " Kin. energy (eV)=" << k/eV
555 << " Sec. energy (eV)=" << secondaryKinetic/eV << G4endl;
556 }
557
558 G4ThreeVector deltaDirection =
559 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
560 Z, Shell,
561 couple->GetMaterial());
562
563 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
564 {
565 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
566
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;
574
575 G4ThreeVector direction;
576 direction.set(finalPx,finalPy,finalPz);
577
578 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
579 }
580 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
581
582 // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
583 G4double deexSecEnergy = 0;
584 for (G4int j=secNumberInit; j < secNumberFinal; j++) {
585 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
586
587 fParticleChangeForGamma->SetProposedKineticEnergy(ekin-bindingEnergy-secondaryKinetic);
588 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy-deexSecEnergy);
589
590 if (secondaryKinetic>0)
591 {
592 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
593 fvect->push_back(dp);
594 }
595
596 }
597}
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
G4int
int G4int
Definition: G4Types.hh:85
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:119
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
Definition: G4ParticleChangeForGamma.hh:139
G4VAtomDeexcitation::GetAtomicShell
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.cc:254
G4VEmAngularDistribution::SampleDirectionForShell
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
Definition: G4VEmAngularDistribution.cc:69
G4VEmModel::GetAngularDistribution
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:611
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:36

◆ SelectFasterComputation()

void G4MicroElecInelasticModel::SelectFasterComputation ( G4bool  input)
inline

Definition at line 177 of file G4MicroElecInelasticModel.hh.

178{
179 fasterCode = input;
180}

◆ TransferedEnergy()

G4double G4MicroElecInelasticModel::TransferedEnergy ( G4ParticleDefinition *  aParticleDefinition,
G4double  incomingParticleEnergy,
G4int  shell,
G4double  random 
)

Definition at line 941 of file G4MicroElecInelasticModel.cc.

945{
946 G4double nrj = 0.;
947
948 G4double valueK1 = 0;
949 G4double valueK2 = 0;
950 G4double valuePROB21 = 0;
951 G4double valuePROB22 = 0;
952 G4double valuePROB12 = 0;
953 G4double valuePROB11 = 0;
954
955 G4double nrjTransf11 = 0;
956 G4double nrjTransf12 = 0;
957 G4double nrjTransf21 = 0;
958 G4double nrjTransf22 = 0;
959
960 G4double maximumEnergyTransfer1 = 0;
961 G4double maximumEnergyTransfer2 = 0;
962 G4double maximumEnergyTransferP = 4.* (electron_mass_c2 / proton_mass_c2) * k;
963 G4double bindingEnergy = SiStructure.Energy(ionizationLevelIndex)*1e6;
964
965 if (particleDefinition == G4Electron::ElectronDefinition())
966 {
967 // k should be in eV
968 std::vector<double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),
969 eTdummyVec.end(),
970 k);
971 std::vector<double>::iterator k1 = k2 - 1;
972
973 /*
974 G4cout << "----> k=" << k
975 << " " << *k1
976 << " " << *k2
977 << " " << random
978 << " " << ionizationLevelIndex
979 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
980 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
981 << G4endl;
982 */
983
984 // SI : the following condition avoids situations where random >last vector element
985 if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
986 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
987 {
988 std::vector<double>::iterator prob12 =
989 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
990 eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
991 random);
992
993 std::vector<double>::iterator prob11 = prob12 - 1;
994
995 std::vector<double>::iterator prob22 =
996 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
997 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
998 random);
999
1000 std::vector<double>::iterator prob21 = prob22 - 1;
1001
1002 valueK1 = *k1;
1003 valueK2 = *k2;
1004 valuePROB21 = *prob21;
1005 valuePROB22 = *prob22;
1006 valuePROB12 = *prob12;
1007 valuePROB11 = *prob11;
1008
1009 /*
1010 G4cout << " " << random << " " << valuePROB11 << " "
1011 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1012 */
1013
1014 // The following condition avoid getting transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
1015 if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
1016 else nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1017 if(valuePROB12 == 1)
1018 {
1019 if ((valueK1+bindingEnergy)/2. > valueK1) maximumEnergyTransfer1=valueK1;
1020 else maximumEnergyTransfer1 = (valueK1+bindingEnergy)/2.;
1021
1022 nrjTransf12 = maximumEnergyTransfer1;
1023 }
1024 else nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1025
1026 if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
1027 else nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1028 if(valuePROB22 == 1)
1029 {
1030 if ((valueK2+bindingEnergy)/2. > valueK2) maximumEnergyTransfer2=valueK2;
1031 else maximumEnergyTransfer2 = (valueK2+bindingEnergy)/2.;
1032
1033 nrjTransf22 = maximumEnergyTransfer2;
1034 }
1035 else nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1036
1037
1038 /*nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1039 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1040 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1041 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];*/
1042
1043 /*
1044 G4cout << " " << ionizationLevelIndex << " "
1045 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1046
1047 G4cout << " " << random << " " << nrjTransf11 << " "
1048 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1049 */
1050
1051 }
1052 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1053 if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
1054 {
1055 std::vector<double>::iterator prob22 =
1056 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1057 eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1058 random);
1059
1060 std::vector<double>::iterator prob21 = prob22 - 1;
1061
1062 valueK1 = *k1;
1063 valueK2 = *k2;
1064 valuePROB21 = *prob21;
1065 valuePROB22 = *prob22;
1066
1067 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1068
1069 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1070 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1071
1072 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1073 valuePROB22,
1074 random,
1075 nrjTransf21,
1076 nrjTransf22);
1077
1078 // zeros are explicitly set
1079
1080 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1081
1082 /*
1083 G4cout << " " << ionizationLevelIndex << " "
1084 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1085
1086 G4cout << " " << random << " " << nrjTransf11 << " "
1087 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1088
1089 G4cout << "ici" << " " << value << G4endl;
1090 */
1091
1092 return value;
1093 }
1094 }
1095 //
1096 else if (particleDefinition == G4Proton::ProtonDefinition())
1097 {
1098 // k should be in eV
1099
1100 std::vector<double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),
1101 pTdummyVec.end(),
1102 k);
1103
1104 std::vector<double>::iterator k1 = k2 - 1;
1105
1106 /*
1107 G4cout << "----> k=" << k
1108 << " " << *k1
1109 << " " << *k2
1110 << " " << random
1111 << " " << ionizationLevelIndex
1112 << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1113 << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1114 << G4endl;
1115 */
1116
1117 // SI : the following condition avoids situations where random > last vector element,
1118 // for eg. when the last element is zero
1119 if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1120 && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1121 {
1122 std::vector<double>::iterator prob12 =
1123 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1124 pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1125 random);
1126
1127 std::vector<double>::iterator prob11 = prob12 - 1;
1128
1129 std::vector<double>::iterator prob22 =
1130 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1131 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1132 random);
1133
1134 std::vector<double>::iterator prob21 = prob22 - 1;
1135
1136 valueK1 = *k1;
1137 valueK2 = *k2;
1138 valuePROB21 = *prob21;
1139 valuePROB22 = *prob22;
1140 valuePROB12 = *prob12;
1141 valuePROB11 = *prob11;
1142
1143 /*
1144 G4cout << " " << random << " " << valuePROB11 << " "
1145 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1146 */
1147
1148 // The following condition avoid getting transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
1149 if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
1150 else nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1151 if(valuePROB12 == 1) nrjTransf12 = maximumEnergyTransferP;
1152 else nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1153 if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
1154 else nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1155 if(valuePROB22 == 1) nrjTransf22 = maximumEnergyTransferP;
1156 else nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1157
1158
1159 /* nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1160 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1161 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1162 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];*/
1163
1164 /*
1165 G4cout << " " << ionizationLevelIndex << " "
1166 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1167
1168 G4cout << " " << random << " " << nrjTransf11 << " "
1169 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1170 */
1171 }
1172
1173 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1174
1175 if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1176 {
1177 std::vector<double>::iterator prob22 =
1178 std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1179 pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1180 random);
1181
1182 std::vector<double>::iterator prob21 = prob22 - 1;
1183
1184 valueK1 = *k1;
1185 valueK2 = *k2;
1186 valuePROB21 = *prob21;
1187 valuePROB22 = *prob22;
1188
1189 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1190
1191 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1192 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1193
1194 G4double interpolatedvalue2 = Interpolate(valuePROB21,
1195 valuePROB22,
1196 random,
1197 nrjTransf21,
1198 nrjTransf22);
1199
1200 // zeros are explicitly set
1201
1202 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1203
1204 /*
1205 G4cout << " " << ionizationLevelIndex << " "
1206 << random << " " <<valueK1 << " " << valueK2 << G4endl;
1207
1208 G4cout << " " << random << " " << nrjTransf11 << " "
1209 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1210
1211 G4cout << "ici" << " " << value << G4endl;
1212 */
1213
1214 return value;
1215 }
1216 }
1217 // End electron and proton cases
1218
1219 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1220 * nrjTransf22;
1221 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1222
1223 if (nrjTransfProduct != 0.)
1224 {
1225 nrj = QuadInterpolator(valuePROB11,
1226 valuePROB12,
1227 valuePROB21,
1228 valuePROB22,
1229 nrjTransf11,
1230 nrjTransf12,
1231 nrjTransf21,
1232 nrjTransf22,
1233 valueK1,
1234 valueK2,
1235 k,
1236 random);
1237 }
1238 //G4cout << nrj << endl;
1239
1240 return nrj;
1241}

Member Data Documentation

◆ fParticleChangeForGamma

G4ParticleChangeForGamma* G4MicroElecInelasticModel::fParticleChangeForGamma
protected

Definition at line 97 of file G4MicroElecInelasticModel.hh.

Referenced by G4MicroElecInelasticModel(), Initialise(), and SampleSecondaries().

The documentation for this class was generated from the following files: