Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Protected Attributes | List of all members
G4MicroElecInelasticModel Class Reference

#include <G4MicroElecInelasticModel.hh>

Inheritance diagram for G4MicroElecInelasticModel:
Inheritance graph
[legend]
Collaboration diagram for G4MicroElecInelasticModel:
Collaboration graph
[legend]

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 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 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)
 
G4int SelectIsotopeNumber (const G4Element *)
 
const G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, 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 *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
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)
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () 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)
 

Protected Attributes

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

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 *)
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 64 of file G4MicroElecInelasticModel.hh.

Constructor & Destructor Documentation

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 }
G4NistManager::FindOrBuildMaterial
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
Definition: G4NistManager.hh:474
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:68
G4cout
G4GLOB_DLL std::ostream G4cout
G4MicroElecInelasticModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4MicroElecInelasticModel.hh:97
G4VEmModel::SetAngularDistribution
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:623
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:780

Here is the call graph for this function:

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 }
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66

Member Function Documentation

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 }
CLHEP::proton_mass_c2
static constexpr double proton_mass_c2
Definition: PhysicalConstants.h:80
cm2
static constexpr double cm2
Definition: G4SIunits.hh:120
G4cout
G4GLOB_DLL std::ostream G4cout
cm
static constexpr double cm
Definition: G4SIunits.hh:119
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4MicroElecCrossSectionDataSet::FindValue
virtual G4double FindValue(G4double e, G4int componentId=0) const
Definition: G4MicroElecCrossSectionDataSet.cc:474
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ionEffectiveCharge::EffectiveCharge
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
Definition: G4ionEffectiveCharge.cc:89
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
FatalException
Definition: G4ExceptionSeverity.hh:60
G4Material::GetBaseMaterial
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66

Here is the call graph for this function:

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:89
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4MicroElecSiStructure::Energy
G4double Energy(G4int level)
Definition: G4MicroElecSiStructure.cc:63
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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 = 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 }
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:112
CLHEP::proton_mass_c2
static constexpr double proton_mass_c2
Definition: PhysicalConstants.h:80
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:633
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:724
G4cout
G4GLOB_DLL std::ostream G4cout
cm
static constexpr double cm
Definition: G4SIunits.hh:119
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4MicroElecInelasticModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4MicroElecInelasticModel.hh:97
G4MicroElecCrossSectionDataSet::LoadData
virtual G4bool LoadData(const G4String &argFileName)
Definition: G4MicroElecCrossSectionDataSet.cc:58
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
FatalException
Definition: G4ExceptionSeverity.hh:60
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:444
G4double
double G4double
Definition: G4Types.hh:76
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:731
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
keV
static constexpr double keV
Definition: G4SIunits.hh:216
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132

Here is the call graph for this function:

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 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
CLHEP::proton_mass_c2
static constexpr double proton_mass_c2
Definition: PhysicalConstants.h:80
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
CLHEP::Hep3Vector::x
double x() const
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4VEmModel::GetAngularDistribution
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:616
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
Definition: G4ParticleChangeForGamma.hh:147
fKShell
Definition: G4AtomicShellEnumerator.hh:53
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
CLHEP::Hep3Vector::z
double z() const
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4MicroElecSiStructure::Energy
G4double Energy(G4int level)
Definition: G4MicroElecSiStructure.cc:63
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4VEmAngularDistribution::SampleDirectionForShell
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
Definition: G4VEmAngularDistribution.cc:60
G4VAtomDeexcitation::GetAtomicShell
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4cout
G4GLOB_DLL std::ostream G4cout
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4MicroElecInelasticModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4MicroElecInelasticModel.hh:97
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::y
double y() const
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:129
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.hh:274
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:37
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150

Here is the call graph for this function:

void G4MicroElecInelasticModel::SelectFasterComputation ( G4bool  input)
inline

Definition at line 177 of file G4MicroElecInelasticModel.hh.

178 {
179  fasterCode = input;
180 }
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 explicitely 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 explicitely 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 }
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
CLHEP::proton_mass_c2
static constexpr double proton_mass_c2
Definition: PhysicalConstants.h:80
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4MicroElecSiStructure::Energy
G4double Energy(G4int level)
Definition: G4MicroElecSiStructure.cc:63
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
value
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4double
double G4double
Definition: G4Types.hh:76
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:37

Here is the call graph for this function:

Member Data Documentation

G4ParticleChangeForGamma* G4MicroElecInelasticModel::fParticleChangeForGamma
protected

Definition at line 97 of file G4MicroElecInelasticModel.hh.

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