Geant4  10.02.p03
G4DNAScreenedRutherfordElasticModel Class Reference

#include <G4DNAScreenedRutherfordElasticModel.hh>

Inheritance diagram for G4DNAScreenedRutherfordElasticModel:
Collaboration diagram for G4DNAScreenedRutherfordElasticModel:

Public Member Functions

 G4DNAScreenedRutherfordElasticModel (const G4ParticleDefinition *p=0, const G4String &nam="DNAScreenedRutherfordElasticModel")
 
virtual ~G4DNAScreenedRutherfordElasticModel ()
 
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)
 
void SetKillBelowThreshold (G4double threshold)
 
G4double GetKillBelowThreshold ()
 
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, const G4ParticleDefinition *, G4double)
 
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 *> *)
 
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 G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
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 G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
const G4IsotopeGetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Protected Attributes

G4ParticleChangeForGamma * fParticleChangeForGamma
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 

Private Member Functions

G4double RutherfordCrossSection (G4double energy, G4double z)
 
G4double ScreeningFactor (G4double energy, G4double z)
 
G4double BrennerZaiderRandomizeCosTheta (G4double k)
 
G4double CalculatePolynomial (G4double k, std::vector< G4double > &vec)
 
G4double ScreenedRutherfordRandomizeCosTheta (G4double k, G4double z)
 
G4DNAScreenedRutherfordElasticModeloperator= (const G4DNAScreenedRutherfordElasticModel &right)
 
 G4DNAScreenedRutherfordElasticModel (const G4DNAScreenedRutherfordElasticModel &)
 

Private Attributes

G4bool fasterCode
 
const std::vector< G4double > * fpWaterDensity
 
G4double killBelowEnergy
 
G4double lowEnergyLimit
 
G4double intermediateEnergyLimit
 
G4double highEnergyLimit
 
G4bool isInitialised
 
G4int verboseLevel
 
std::vector< G4doublebetaCoeff
 
std::vector< G4doubledeltaCoeff
 
std::vector< G4doublegamma035_10Coeff
 
std::vector< G4doublegamma10_100Coeff
 
std::vector< G4doublegamma100_200Coeff
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 39 of file G4DNAScreenedRutherfordElasticModel.hh.

Constructor & Destructor Documentation

◆ G4DNAScreenedRutherfordElasticModel() [1/2]

G4DNAScreenedRutherfordElasticModel::G4DNAScreenedRutherfordElasticModel ( const G4ParticleDefinition p = 0,
const G4String nam = "DNAScreenedRutherfordElasticModel" 
)

Definition at line 40 of file G4DNAScreenedRutherfordElasticModel.cc.

41  :
42  G4VEmModel(nam), isInitialised(false)
43 {
44  // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
45  fpWaterDensity = 0;
46 
47  killBelowEnergy = 9 * eV;
48  lowEnergyLimit = 0 * eV;
49  intermediateEnergyLimit = 200 * eV; // Switch between two final state models
50  highEnergyLimit = 1. * MeV;
53 
54  verboseLevel = 0;
55  // Verbosity scale:
56  // 0 = nothing
57  // 1 = warning for energy non-conservation
58  // 2 = details of energy budget
59  // 3 = calculation of cross sections, file openings, sampling of atoms
60  // 4 = entering in methods
61 
62  if (verboseLevel > 0)
63  {
64  G4cout << "Screened Rutherford Elastic model is constructed " << G4endl<< "Energy range: "
65  << lowEnergyLimit / eV << " eV - "
66  << highEnergyLimit / MeV << " MeV"
67  << G4endl;
68  }
70 
71  // Selection of computation method
72  // We do not recommend "true" usage with the current cumul. proba. settings
73  fasterCode = false;
74 }
static const double MeV
Definition: G4SIunits.hh:211
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:69
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
G4GLOB_DLL std::ostream G4cout
static const double eV
Definition: G4SIunits.hh:212
#define G4endl
Definition: G4ios.hh:61
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
Here is the call graph for this function:

◆ ~G4DNAScreenedRutherfordElasticModel()

G4DNAScreenedRutherfordElasticModel::~G4DNAScreenedRutherfordElasticModel ( )
virtual

Definition at line 78 of file G4DNAScreenedRutherfordElasticModel.cc.

79 {
80 }

◆ G4DNAScreenedRutherfordElasticModel() [2/2]

G4DNAScreenedRutherfordElasticModel::G4DNAScreenedRutherfordElasticModel ( const G4DNAScreenedRutherfordElasticModel )
private

Member Function Documentation

◆ BrennerZaiderRandomizeCosTheta()

G4double G4DNAScreenedRutherfordElasticModel::BrennerZaiderRandomizeCosTheta ( G4double  k)
private

Definition at line 336 of file G4DNAScreenedRutherfordElasticModel.cc.

337 {
338  // d sigma_el 1 beta(K)
339  // ------------ (K) ~ --------------------------------- + ---------------------------------
340  // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2
341  //
342  // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
343  //
344  // Phys. Med. Biol. 29 N.4 (1983) 443-447
345 
346  // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
347 
348  k /= eV;
349 
350  G4double beta = std::exp(CalculatePolynomial(k, betaCoeff));
351  G4double delta = std::exp(CalculatePolynomial(k, deltaCoeff));
352  G4double gamma;
353 
354  if (k > 100.)
355  {
357  // Only in this case it is not the exponent of the polynomial
358  }
359  else
360  {
361  if (k > 10)
362  {
363  gamma = std::exp(CalculatePolynomial(k, gamma10_100Coeff));
364  }
365  else
366  {
367  gamma = std::exp(CalculatePolynomial(k, gamma035_10Coeff));
368  }
369  }
370 
371  // ***** Original method
372 
373  if (!fasterCode)
374  {
375 
376  G4double oneOverMax = 1.
377  / (1. / (4. * gamma * gamma) + beta
378  / ((2. + 2. * delta) * (2. + 2. * delta)));
379 
380  G4double cosTheta = 0.;
381  G4double leftDenominator = 0.;
382  G4double rightDenominator = 0.;
383  G4double fCosTheta = 0.;
384 
385  do
386  {
387  cosTheta = 2. * G4UniformRand()- 1.;
388 
389  leftDenominator = (1. + 2.*gamma - cosTheta);
390  rightDenominator = (1. + 2.*delta + cosTheta);
391  if ( (leftDenominator * rightDenominator) != 0. )
392  {
393  fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
394  }
395  }
396  while (fCosTheta < G4UniformRand());
397 
398  return cosTheta;
399  }
400 
401  // ***** Alternative method using cumulative probability
402 
403  if (fasterCode)
404  {
405 
406  G4double cosTheta = -1;
407  G4double cumul = 0;
408  G4double value = 0;
409  G4double leftDenominator = 0.;
410  G4double rightDenominator = 0.;
411 
412  // Number of integration steps in the -1,1 range
413  G4int iMax=200;
414 
415  G4double random = G4UniformRand();
416 
417  // Cumulate differential cross section
418  for (G4int i=0; i<iMax; i++)
419  {
420  cosTheta = -1 + i*2./(iMax-1);
421  leftDenominator = (1. + 2.*gamma - cosTheta);
422  rightDenominator = (1. + 2.*delta + cosTheta);
423  if ( (leftDenominator * rightDenominator) != 0. )
424  {
425  cumul = cumul + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
426  }
427  }
428 
429  // Select cosTheta
430  for (G4int i=0; i<iMax; i++)
431  {
432  cosTheta = -1 + i*2./(iMax-1);
433  leftDenominator = (1. + 2.*gamma - cosTheta);
434  rightDenominator = (1. + 2.*delta + cosTheta);
435  if (cumul !=0 && (leftDenominator * rightDenominator) != 0.)
436  value = value + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)) / cumul;
437  if (random < value) break;
438  }
439 
440  return cosTheta;
441  }
442 
443  return 0.;
444 }
G4double CalculatePolynomial(G4double k, std::vector< G4double > &vec)
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
static const double eV
Definition: G4SIunits.hh:212
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalculatePolynomial()

G4double G4DNAScreenedRutherfordElasticModel::CalculatePolynomial ( G4double  k,
std::vector< G4double > &  vec 
)
private

Definition at line 448 of file G4DNAScreenedRutherfordElasticModel.cc.

451 {
452  // Sum_{i=0}^{size-1} vector_i k^i
453  //
454  // Phys. Med. Biol. 29 N.4 (1983) 443-447
455 
456  G4double result = 0.;
457  size_t size = vec.size();
458 
459  while (size > 0)
460  {
461  size--;
462 
463  result *= k;
464  result += vec[size];
465  }
466 
467  return result;
468 }
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ CrossSectionPerVolume()

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

Reimplemented from G4VEmModel.

Definition at line 169 of file G4DNAScreenedRutherfordElasticModel.cc.

174 {
175  if (verboseLevel > 3)
176  G4cout << "Calling CrossSectionPerVolume() of G4DNAScreenedRutherfordElasticModel"
177  << G4endl;
178 
179  // Calculate total cross section for model
180 
181  G4double sigma=0;
182 
183  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
184 
185  if(waterDensity!= 0.0)
186  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
187  {
188 
189  if (ekin < highEnergyLimit)
190  {
191 
192  if (ekin < killBelowEnergy) return DBL_MAX;
193 
194  G4double z = 10.;
195  G4double n = ScreeningFactor(ekin,z);
196  G4double crossSection = RutherfordCrossSection(ekin, z);
197  sigma = pi * crossSection / (n * (n + 1.));
198  }
199 
200  if (verboseLevel > 2)
201  {
202  G4cout << "__________________________________" << G4endl;
203  G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO START" << G4endl;
204  G4cout << "=== Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
205  G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
206  G4cout << "=== Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
207  // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
208  G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO END" << G4endl;
209  }
210 
211  }
212 
213  return sigma*waterDensity;
214 }
static const double cm
Definition: G4SIunits.hh:118
G4double RutherfordCrossSection(G4double energy, G4double z)
size_t GetIndex() const
Definition: G4Material.hh:262
Char_t n[5]
G4GLOB_DLL std::ostream G4cout
static const double eV
Definition: G4SIunits.hh:212
G4double ScreeningFactor(G4double energy, G4double z)
static const double pi
Definition: G4SIunits.hh:74
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
Here is the call graph for this function:

◆ GetKillBelowThreshold()

G4double G4DNAScreenedRutherfordElasticModel::GetKillBelowThreshold ( )
inline

Definition at line 64 of file G4DNAScreenedRutherfordElasticModel.hh.

Here is the call graph for this function:

◆ Initialise()

void G4DNAScreenedRutherfordElasticModel::Initialise ( const G4ParticleDefinition ,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 84 of file G4DNAScreenedRutherfordElasticModel.cc.

86 {
87 
88  if (verboseLevel > 3)
89  G4cout << "Calling G4DNAScreenedRutherfordElasticModel::Initialise()"
90  << G4endl;
91 
92  // Energy limits
93 
95  {
96  G4cout << "G4DNAScreenedRutherfordElasticModel: low energy limit increased from " <<
97  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
99  }
100 
102  {
103  G4cout << "G4DNAScreenedRutherfordElasticModel: high energy limit decreased from " <<
104  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
106  }
107 
108  // Constants for final state by Brenner & Zaider
109  // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
110  // Added clear for MT
111 
112  betaCoeff.clear();
113  betaCoeff.push_back(7.51525);
114  betaCoeff.push_back(-0.41912);
115  betaCoeff.push_back(7.2017E-3);
116  betaCoeff.push_back(-4.646E-5);
117  betaCoeff.push_back(1.02897E-7);
118 
119  deltaCoeff.clear();
120  deltaCoeff.push_back(2.9612);
121  deltaCoeff.push_back(-0.26376);
122  deltaCoeff.push_back(4.307E-3);
123  deltaCoeff.push_back(-2.6895E-5);
124  deltaCoeff.push_back(5.83505E-8);
125 
126  gamma035_10Coeff.clear();
127  gamma035_10Coeff.push_back(-1.7013);
128  gamma035_10Coeff.push_back(-1.48284);
129  gamma035_10Coeff.push_back(0.6331);
130  gamma035_10Coeff.push_back(-0.10911);
131  gamma035_10Coeff.push_back(8.358E-3);
132  gamma035_10Coeff.push_back(-2.388E-4);
133 
134  gamma10_100Coeff.clear();
135  gamma10_100Coeff.push_back(-3.32517);
136  gamma10_100Coeff.push_back(0.10996);
137  gamma10_100Coeff.push_back(-4.5255E-3);
138  gamma10_100Coeff.push_back(5.8372E-5);
139  gamma10_100Coeff.push_back(-2.4659E-7);
140 
141  gamma100_200Coeff.clear();
142  gamma100_200Coeff.push_back(2.4775E-2);
143  gamma100_200Coeff.push_back(-2.96264E-5);
144  gamma100_200Coeff.push_back(-1.20655E-7);
145 
146  //
147 
148  if( verboseLevel>0 )
149  {
150  G4cout << "Screened Rutherford elastic model is initialized " << G4endl
151  << "Energy range: "
152  << LowEnergyLimit() / eV << " eV - "
153  << HighEnergyLimit() / MeV << " MeV"
154  << G4endl;
155  }
156 
157  // Initialize water density pointer
159 
160  if (isInitialised)
161  { return;}
163  isInitialised = true;
164 
165 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
static const double MeV
Definition: G4SIunits.hh:211
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:604
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
G4GLOB_DLL std::ostream G4cout
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
static const double eV
Definition: G4SIunits.hh:212
static G4DNAMolecularMaterial * Instance()
#define G4endl
Definition: G4ios.hh:61
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
Here is the call graph for this function:

◆ operator=()

G4DNAScreenedRutherfordElasticModel& G4DNAScreenedRutherfordElasticModel::operator= ( const G4DNAScreenedRutherfordElasticModel right)
private

◆ RutherfordCrossSection()

G4double G4DNAScreenedRutherfordElasticModel::RutherfordCrossSection ( G4double  energy,
G4double  z 
)
private

Definition at line 218 of file G4DNAScreenedRutherfordElasticModel.cc.

220 {
221  //
222  // e^4 / K + m_e c^2 \^2
223  // sigma_Ruth(K) = Z (Z+1) -------------------- | --------------------- |
224  // (4 pi epsilon_0)^2 \ K * (K + 2 m_e c^2) /
225  //
226  // Where K is the electron non-relativistic kinetic energy
227  //
228  // NIM 155, pp. 145-156, 1978
229 
230  G4double length = (e_squared * (k + electron_mass_c2))
231  / (4 * pi * epsilon0 * k * (k + 2 * electron_mass_c2));
232  G4double cross = z * (z + 1) * length * length;
233 
234  return cross;
235 }
float electron_mass_c2
Definition: hepunit.py:274
static const double pi
Definition: G4SIunits.hh:74
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ SampleSecondaries()

void G4DNAScreenedRutherfordElasticModel::SampleSecondaries ( std::vector< G4DynamicParticle *> *  ,
const G4MaterialCutsCouple ,
const G4DynamicParticle aDynamicElectron,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 273 of file G4DNAScreenedRutherfordElasticModel.cc.

279 {
280 
281  if (verboseLevel > 3)
282  {
283  G4cout << "Calling SampleSecondaries() of G4DNAScreenedRutherfordElasticModel"
284  << G4endl;
285  }
286 
287  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
288 
289  if (electronEnergy0 < killBelowEnergy)
290  {
291  fParticleChangeForGamma->SetProposedKineticEnergy(0.);
292  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
293  fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
294  return;
295  }
296 
297  G4double cosTheta = 0.;
298 
299  if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
300  {
301  if (electronEnergy0<intermediateEnergyLimit)
302  {
303  if (verboseLevel > 3) G4cout << "---> Using Brenner & Zaider model" << G4endl;
304  cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
305  }
306 
307  if (electronEnergy0>=intermediateEnergyLimit)
308  {
309  if (verboseLevel > 3) G4cout << "---> Using Screened Rutherford model" << G4endl;
310  G4double z = 10.;
311  cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
312  }
313 
314  G4double phi = 2. * pi * G4UniformRand();
315 
316  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
317  G4ThreeVector xVers = zVers.orthogonal();
318  G4ThreeVector yVers = zVers.cross(xVers);
319 
320  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
321  G4double yDir = xDir;
322  xDir *= std::cos(phi);
323  yDir *= std::sin(phi);
324 
325  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
326 
327  fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
328 
329  fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
330  }
331 
332 }
G4double GetKineticEnergy() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
Hep3Vector cross(const Hep3Vector &) const
Hep3Vector orthogonal() const
static const double pi
Definition: G4SIunits.hh:74
const G4ThreeVector & GetMomentumDirection() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double ScreenedRutherfordRandomizeCosTheta(G4double k, G4double z)
Here is the call graph for this function:

◆ ScreenedRutherfordRandomizeCosTheta()

G4double G4DNAScreenedRutherfordElasticModel::ScreenedRutherfordRandomizeCosTheta ( G4double  k,
G4double  z 
)
private

Definition at line 472 of file G4DNAScreenedRutherfordElasticModel.cc.

474 {
475 
476  // d sigma_el sigma_Ruth(K)
477  // ------------ (K) ~ -----------------------------
478  // d Omega (1 + 2 n(K) - cos(theta))^2
479  //
480  // We extract cos(theta) distributed as (1 + 2 n(K) - cos(theta))^-2
481  //
482  // Maximum is for theta=0: 1/(4 n(K)^2) (When n(K) is positive, that is always satisfied within the validity of the process)
483  //
484  // Phys. Med. Biol. 45 (2000) 3171-3194
485 
486  // ***** Original method
487 
488  if (!fasterCode)
489  {
490 
491  G4double n = ScreeningFactor(k, z);
492 
493  G4double oneOverMax = (4. * n * n);
494 
495  G4double cosTheta = 0.;
496  G4double fCosTheta;
497 
498  do
499  {
500  cosTheta = 2. * G4UniformRand()- 1.;
501  fCosTheta = (1 + 2.*n - cosTheta);
502  if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta);
503  }
504  while (fCosTheta < G4UniformRand());
505 
506  return cosTheta;
507  }
508 
509  // ***** Alternative method using cumulative probability
510  if (fasterCode)
511  {
512 
513  G4double cosTheta = -1;
514  G4double cumul = 0;
515  G4double value = 0;
516  G4double n = ScreeningFactor(k, z);
517  G4double fCosTheta;
518 
519  // Number of integration steps in the -1,1 range
520  G4int iMax=200;
521 
522  G4double random = G4UniformRand();
523 
524  // Cumulate differential cross section
525  for (G4int i=0; i<iMax; i++)
526  {
527  cosTheta = -1 + i*2./(iMax-1);
528  fCosTheta = (1 + 2.*n - cosTheta);
529  if (fCosTheta !=0.) cumul = cumul + 1./(fCosTheta*fCosTheta);
530  }
531 
532  // Select cosTheta
533  for (G4int i=0; i<iMax; i++)
534  {
535  cosTheta = -1 + i*2./(iMax-1);
536  fCosTheta = (1 + 2.*n - cosTheta);
537  if (cumul !=0.) value = value + (1./(fCosTheta*fCosTheta)) / cumul;
538  if (random < value) break;
539  }
540  return cosTheta;
541  }
542 
543  return 0.;
544 }
int G4int
Definition: G4Types.hh:78
Char_t n[5]
#define G4UniformRand()
Definition: Randomize.hh:97
G4double ScreeningFactor(G4double energy, G4double z)
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ScreeningFactor()

G4double G4DNAScreenedRutherfordElasticModel::ScreeningFactor ( G4double  energy,
G4double  z 
)
private

Definition at line 239 of file G4DNAScreenedRutherfordElasticModel.cc.

241 {
242  //
243  // alpha_1 + beta_1 ln(K/eV) constK Z^(2/3)
244  // n(T) = -------------------------- -----------------
245  // K/(m_e c^2) 2 + K/(m_e c^2)
246  //
247  // Where K is the electron non-relativistic kinetic energy
248  //
249  // n(T) > 0 for T < ~ 400 MeV
250  //
251  // NIM 155, pp. 145-156, 1978
252  // Formulae (2) and (5)
253 
254  const G4double alpha_1(1.64);
255  const G4double beta_1(-0.0825);
256  const G4double constK(1.7E-5);
257 
258  G4double numerator = (alpha_1 + beta_1 * std::log(k / eV)) * constK
259  * std::pow(z, 2. / 3.);
260 
261  k /= electron_mass_c2;
262 
263  G4double denominator = k * (2 + k);
264 
265  G4double value = 0.;
266  if (denominator > 0.) value = numerator / denominator;
267 
268  return value;
269 }
float electron_mass_c2
Definition: hepunit.py:274
static const double eV
Definition: G4SIunits.hh:212
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ SelectFasterComputation()

void G4DNAScreenedRutherfordElasticModel::SelectFasterComputation ( G4bool  input)
inline

Definition at line 120 of file G4DNAScreenedRutherfordElasticModel.hh.

Here is the caller graph for this function:

◆ SetKillBelowThreshold()

void G4DNAScreenedRutherfordElasticModel::SetKillBelowThreshold ( G4double  threshold)
inline

Definition at line 113 of file G4DNAScreenedRutherfordElasticModel.hh.

114 {
115  killBelowEnergy = threshold;
116  if (threshold < 9*CLHEP::eV)
117  G4Exception ("*** WARNING : the G4DNAScreenedRutherfordElasticModel class is not validated below 9 eV !","",JustWarning,"") ;
118 }
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const double eV
Here is the call graph for this function:

Member Data Documentation

◆ betaCoeff

std::vector<G4double> G4DNAScreenedRutherfordElasticModel::betaCoeff
private

Definition at line 96 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ deltaCoeff

std::vector<G4double> G4DNAScreenedRutherfordElasticModel::deltaCoeff
private

Definition at line 97 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ fasterCode

G4bool G4DNAScreenedRutherfordElasticModel::fasterCode
private

Definition at line 74 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ fParticleChangeForGamma

G4ParticleChangeForGamma* G4DNAScreenedRutherfordElasticModel::fParticleChangeForGamma
protected

Definition at line 70 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ fpWaterDensity

const std::vector<G4double>* G4DNAScreenedRutherfordElasticModel::fpWaterDensity
private

Definition at line 77 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ gamma035_10Coeff

std::vector<G4double> G4DNAScreenedRutherfordElasticModel::gamma035_10Coeff
private

Definition at line 98 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ gamma100_200Coeff

std::vector<G4double> G4DNAScreenedRutherfordElasticModel::gamma100_200Coeff
private

Definition at line 100 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ gamma10_100Coeff

std::vector<G4double> G4DNAScreenedRutherfordElasticModel::gamma10_100Coeff
private

Definition at line 99 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ highEnergyLimit

G4double G4DNAScreenedRutherfordElasticModel::highEnergyLimit
private

Definition at line 82 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ intermediateEnergyLimit

G4double G4DNAScreenedRutherfordElasticModel::intermediateEnergyLimit
private

Definition at line 81 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ isInitialised

G4bool G4DNAScreenedRutherfordElasticModel::isInitialised
private

Definition at line 83 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ killBelowEnergy

G4double G4DNAScreenedRutherfordElasticModel::killBelowEnergy
private

Definition at line 79 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ lowEnergyLimit

G4double G4DNAScreenedRutherfordElasticModel::lowEnergyLimit
private

Definition at line 80 of file G4DNAScreenedRutherfordElasticModel.hh.

◆ verboseLevel

G4int G4DNAScreenedRutherfordElasticModel::verboseLevel
private

Definition at line 84 of file G4DNAScreenedRutherfordElasticModel.hh.


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