Geant4  10.02.p03
G4DNARuddIonisationModel Class Reference

#include <G4DNARuddIonisationModel.hh>

Inheritance diagram for G4DNARuddIonisationModel:
Collaboration diagram for G4DNARuddIonisationModel:

Public Member Functions

 G4DNARuddIonisationModel (const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationModel")
 
virtual ~G4DNARuddIonisationModel ()
 
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 SelectStationary (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 Types

typedef std::map< G4String, G4String, std::less< G4String > > MapFile
 
typedef std::map< G4String, G4DNACrossSectionDataSet *, std::less< G4String > > MapData
 

Private Member Functions

G4double RandomizeEjectedElectronEnergy (G4ParticleDefinition *particleDefinition, G4double incomingParticleEnergy, G4int shell)
 
G4double DifferentialCrossSection (G4ParticleDefinition *particleDefinition, G4double k, G4double energyTransfer, G4int shell)
 
G4double CorrectionFactor (G4ParticleDefinition *particleDefinition, G4double k)
 
G4double S_1s (G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
 
G4double S_2s (G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
 
G4double S_2p (G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
 
G4double R (G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
 
G4double PartialCrossSection (const G4Track &track)
 
G4double Sum (G4double energy, const G4String &particle)
 
G4int RandomSelect (G4double energy, const G4String &particle)
 
G4DNARuddIonisationModeloperator= (const G4DNARuddIonisationModel &right)
 
 G4DNARuddIonisationModel (const G4DNARuddIonisationModel &)
 

Private Attributes

G4bool statCode
 
const std::vector< G4double > * fpWaterDensity
 
G4VAtomDeexcitationfAtomDeexcitation
 
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
 
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
 
G4double lowEnergyLimitForZ1
 
G4double lowEnergyLimitForZ2
 
G4double lowEnergyLimitOfModelForZ1
 
G4double lowEnergyLimitOfModelForZ2
 
G4double killBelowEnergyForZ1
 
G4double killBelowEnergyForZ2
 
G4bool isInitialised
 
G4int verboseLevel
 
MapFile tableFile
 
MapData tableData
 
G4DNAWaterIonisationStructure waterStructure
 
G4double slaterEffectiveCharge [3]
 
G4double sCoefficient [3]
 

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 46 of file G4DNARuddIonisationModel.hh.

Member Typedef Documentation

◆ MapData

Definition at line 104 of file G4DNARuddIonisationModel.hh.

◆ MapFile

typedef std::map<G4String,G4String,std::less<G4String> > G4DNARuddIonisationModel::MapFile
private

Definition at line 101 of file G4DNARuddIonisationModel.hh.

Constructor & Destructor Documentation

◆ G4DNARuddIonisationModel() [1/2]

G4DNARuddIonisationModel::G4DNARuddIonisationModel ( const G4ParticleDefinition p = 0,
const G4String nam = "DNARuddIonisationModel" 
)

Definition at line 46 of file G4DNARuddIonisationModel.cc.

47  :
48  G4VEmModel(nam), isInitialised(false)
49 {
50  fpWaterDensity = 0;
51 
52  slaterEffectiveCharge[0] = 0.;
53  slaterEffectiveCharge[1] = 0.;
54  slaterEffectiveCharge[2] = 0.;
55  sCoefficient[0] = 0.;
56  sCoefficient[1] = 0.;
57  sCoefficient[2] = 0.;
58 
59  lowEnergyLimitForZ1 = 0 * eV;
60  lowEnergyLimitForZ2 = 0 * eV;
65 
66  verboseLevel = 0;
67  // Verbosity scale:
68  // 0 = nothing
69  // 1 = warning for energy non-conservation
70  // 2 = details of energy budget
71  // 3 = calculation of cross sections, file openings, sampling of atoms
72  // 4 = entering in methods
73 
74  if (verboseLevel > 0)
75  {
76  G4cout << "Rudd ionisation model is constructed " << G4endl;
77  }
78 
79  // define default angular generator
81 
82  //Mark this model as "applicable" for atomic deexcitation
83  SetDeexcitationFlag(true);
86 
87  // Selection of stationary mode
88 
89  statCode = false;
90 }
G4ParticleChangeForGamma * fParticleChangeForGamma
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:69
G4GLOB_DLL std::ostream G4cout
const std::vector< G4double > * fpWaterDensity
static const double eV
Definition: G4SIunits.hh:212
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:624
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:213
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:781
G4VAtomDeexcitation * fAtomDeexcitation
Here is the call graph for this function:

◆ ~G4DNARuddIonisationModel()

G4DNARuddIonisationModel::~G4DNARuddIonisationModel ( )
virtual

Definition at line 94 of file G4DNARuddIonisationModel.cc.

95 {
96  // Cross section
97 
98  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
99  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
100  {
101  G4DNACrossSectionDataSet* table = pos->second;
102  delete table;
103  }
104 
105  // The following removal is forbidden since G4VEnergyLossmodel takes care of deletion
106  // Coverity however will signal this as an error
107  //if (fAtomDeexcitation) {delete fAtomDeexcitation;}
108 
109 }
static const G4double pos

◆ G4DNARuddIonisationModel() [2/2]

G4DNARuddIonisationModel::G4DNARuddIonisationModel ( const G4DNARuddIonisationModel )
private

Member Function Documentation

◆ CorrectionFactor()

G4double G4DNARuddIonisationModel::CorrectionFactor ( G4ParticleDefinition particleDefinition,
G4double  k 
)
private

Definition at line 992 of file G4DNARuddIonisationModel.cc.

994 {
997 
998  if (particleDefinition == G4Proton::Proton())
999  {
1000  return (1.);
1001  } else if (particleDefinition == instance->GetIon("hydrogen"))
1002  {
1003  G4double value = (std::log10(k / eV) - 4.2) / 0.5;
1004  // The following values are provided by M. Dingfelder (priv. comm)
1005  return ((0.6 / (1 + std::exp(value))) + 0.9);
1006  } else
1007  {
1008  return (1.);
1009  }
1010 }
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4DNAGenericIonsManager * Instance(void)
static const double eV
Definition: G4SIunits.hh:212
static MCTruthManager * instance
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition * GetIon(const G4String &name)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CrossSectionPerVolume()

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

Reimplemented from G4VEmModel.

Definition at line 288 of file G4DNARuddIonisationModel.cc.

293 {
294  if (verboseLevel > 3)
295  {
296  G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationModel"
297  << G4endl;
298  }
299 
300  // Calculate total cross section for model
301 
304 
305  if (
306  particleDefinition != G4Proton::ProtonDefinition()
307  &&
308  particleDefinition != instance->GetIon("hydrogen")
309  &&
310  particleDefinition != instance->GetIon("alpha++")
311  &&
312  particleDefinition != instance->GetIon("alpha+")
313  &&
314  particleDefinition != instance->GetIon("helium")
315  )
316 
317  return 0;
318 
319  G4double lowLim = 0;
320 
321  if ( particleDefinition == G4Proton::ProtonDefinition()
322  || particleDefinition == instance->GetIon("hydrogen")
323  )
324 
326 
327  if ( particleDefinition == instance->GetIon("alpha++")
328  || particleDefinition == instance->GetIon("alpha+")
329  || particleDefinition == instance->GetIon("helium")
330  )
331 
333 
334  G4double highLim = 0;
335  G4double sigma=0;
336 
337  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
338 
339  if(waterDensity!= 0.0)
340  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
341  {
342  const G4String& particleName = particleDefinition->GetParticleName();
343 
344  // SI - the following is useless since lowLim is already defined
345  /*
346  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
347  pos1 = lowEnergyLimit.find(particleName);
348 
349  if (pos1 != lowEnergyLimit.end())
350  {
351  lowLim = pos1->second;
352  }
353  */
354 
355  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
356  pos2 = highEnergyLimit.find(particleName);
357 
358  if (pos2 != highEnergyLimit.end())
359  {
360  highLim = pos2->second;
361  }
362 
363  if (k <= highLim)
364  {
365  //SI : XS must not be zero otherwise sampling of secondaries method ignored
366 
367  if (k < lowLim) k = lowLim;
368 
369  //
370 
371  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
372  pos = tableData.find(particleName);
373 
374  if (pos != tableData.end())
375  {
376  G4DNACrossSectionDataSet* table = pos->second;
377  if (table != 0)
378  {
379  sigma = table->FindValue(k);
380  }
381  }
382  else
383  {
384  G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume","em0002",
385  FatalException,"Model not applicable to particle type.");
386  }
387 
388  } // if (k >= lowLim && k < highLim)
389 
390  if (verboseLevel > 2)
391  {
392  G4cout << "__________________________________" << G4endl;
393  G4cout << "G4DNARuddIonisationModel - XS INFO START" << G4endl;
394  G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
395  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
396  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
397  // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
398  G4cout << "G4DNARuddIonisationModel - XS INFO END" << G4endl;
399  }
400 
401  } // if (waterMaterial)
402  else
403  {
404  if (verboseLevel > 2)
405  {
406  G4cout << "Warning : RuddIonisationModel: WATER DENSITY IS NULL" << G4endl;
407  }
408  }
409 
410  return sigma*waterDensity;
411  // return sigma*material->GetAtomicNumDensityVector()[1];
412 
413 }
static const double cm
Definition: G4SIunits.hh:118
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
size_t GetIndex() const
Definition: G4Material.hh:262
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4GLOB_DLL std::ostream G4cout
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
static G4DNAGenericIonsManager * Instance(void)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const double eV
Definition: G4SIunits.hh:212
#define G4endl
Definition: G4ios.hh:61
static MCTruthManager * instance
double G4double
Definition: G4Types.hh:76
static const G4double pos
G4ParticleDefinition * GetIon(const G4String &name)
Here is the call graph for this function:

◆ DifferentialCrossSection()

G4double G4DNARuddIonisationModel::DifferentialCrossSection ( G4ParticleDefinition particleDefinition,
G4double  k,
G4double  energyTransfer,
G4int  shell 
)
private

Definition at line 690 of file G4DNARuddIonisationModel.cc.

694 {
695  // Shells ids are 0 1 2 3 4 (4 is k shell)
696  // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
697  // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
698  //
699  // ds S F1(nu) + w * F2(nu)
700  // ---- = G(k) * ---- -------------------------------------------
701  // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
702  //
703  // w is the secondary electron kinetic Energy in eV
704  //
705  // All the other parameters can be found in Rudd's Papers
706  //
707  // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
708  // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
709  //
710 
711  const G4int j = ionizationLevelIndex;
712 
713  G4double A1;
714  G4double B1;
715  G4double C1;
716  G4double D1;
717  G4double E1;
718  G4double A2;
719  G4double B2;
720  G4double C2;
721  G4double D2;
722  G4double alphaConst;
723 
724  // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
725  // The following values are provided by M. dingfelder (priv. comm)
726  const G4double Bj[5] = { 12.60 * eV, 14.70 * eV, 18.40 * eV, 32.20 * eV, 540
727  * eV };
728 
729  if (j == 4)
730  {
731  //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
732  A1 = 1.25;
733  B1 = 0.5;
734  C1 = 1.00;
735  D1 = 1.00;
736  E1 = 3.00;
737  A2 = 1.10;
738  B2 = 1.30;
739  C2 = 1.00;
740  D2 = 0.00;
741  alphaConst = 0.66;
742  } else
743  {
744  //Data For Liquid Water from Dingfelder (Protons in Water)
745  A1 = 1.02;
746  B1 = 82.0;
747  C1 = 0.45;
748  D1 = -0.80;
749  E1 = 0.38;
750  A2 = 1.07;
751  // Value provided by M. Dingfelder (priv. comm)
752  B2 = 11.6;
753  //
754  C2 = 0.60;
755  D2 = 0.04;
756  alphaConst = 0.64;
757  }
758 
759  const G4double n = 2.;
760  const G4double Gj[5] = { 0.99, 1.11, 1.11, 0.52, 1. };
761 
764 
765  G4double wBig = (energyTransfer
766  - waterStructure.IonisationEnergy(ionizationLevelIndex));
767  if (wBig < 0)
768  return 0.;
769 
770  G4double w = wBig / Bj[ionizationLevelIndex];
771  // Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
772  if (j == 4)
773  w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
774 
775  G4double Ry = 13.6 * eV;
776 
777  G4double tau = 0.;
778 
779  G4bool isProtonOrHydrogen = false;
780  G4bool isHelium = false;
781 
782  if (particleDefinition == G4Proton::ProtonDefinition()
783  || particleDefinition == instance->GetIon("hydrogen"))
784  {
785  isProtonOrHydrogen = true;
786  tau = (electron_mass_c2 / proton_mass_c2) * k;
787  }
788 
789  else if (particleDefinition == instance->GetIon("helium")
790  || particleDefinition == instance->GetIon("alpha+")
791  || particleDefinition == instance->GetIon("alpha++"))
792  {
793  isHelium = true;
794  tau = (0.511 / 3728.) * k;
795  }
796 
797  G4double S = 4. * pi * Bohr_radius * Bohr_radius * n
798  * std::pow((Ry / Bj[ionizationLevelIndex]), 2);
799  if (j == 4)
800  S = 4. * pi * Bohr_radius * Bohr_radius * n
801  * std::pow((Ry / waterStructure.IonisationEnergy(ionizationLevelIndex)),
802  2);
803 
804  G4double v2 = tau / Bj[ionizationLevelIndex];
805  if (j == 4)
806  v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
807 
808  G4double v = std::sqrt(v2);
809  G4double wc = 4. * v2 - 2. * v - (Ry / (4. * Bj[ionizationLevelIndex]));
810  if (j == 4)
811  wc = 4. * v2 - 2. * v
812  - (Ry / (4. * waterStructure.IonisationEnergy(ionizationLevelIndex)));
813 
814  G4double L1 = (C1 * std::pow(v, (D1))) / (1. + E1 * std::pow(v, (D1 + 4.)));
815  G4double L2 = C2 * std::pow(v, (D2));
816  G4double H1 = (A1 * std::log(1. + v2)) / (v2 + (B1 / v2));
817  G4double H2 = (A2 / v2) + (B2 / (v2 * v2));
818 
819  G4double F1 = L1 + H1;
820  G4double F2 = (L2 * H2) / (L2 + H2);
821 
822  G4double sigma =
823  CorrectionFactor(particleDefinition, k) * Gj[j]
824  * (S / Bj[ionizationLevelIndex])
825  * ((F1 + w * F2)
826  / (std::pow((1. + w), 3)
827  * (1. + std::exp(alphaConst * (w - wc) / v))));
828 
829  if (j == 4)
830  sigma = CorrectionFactor(particleDefinition, k) * Gj[j]
831  * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
832  * ((F1 + w * F2)
833  / (std::pow((1. + w), 3)
834  * (1. + std::exp(alphaConst * (w - wc) / v))));
835 
836  if ((particleDefinition == instance->GetIon("hydrogen"))
837  && (ionizationLevelIndex == 4))
838  {
839  // sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
840  sigma = Gj[j] * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
841  * ((F1 + w * F2)
842  / (std::pow((1. + w), 3)
843  * (1. + std::exp(alphaConst * (w - wc) / v))));
844  }
845 // if ( particleDefinition == G4Proton::ProtonDefinition()
846 // || particleDefinition == instance->GetIon("hydrogen")
847 // )
848  if (isProtonOrHydrogen)
849  {
850  return (sigma);
851  }
852 
853  if (particleDefinition == instance->GetIon("alpha++"))
854  {
855  slaterEffectiveCharge[0] = 0.;
856  slaterEffectiveCharge[1] = 0.;
857  slaterEffectiveCharge[2] = 0.;
858  sCoefficient[0] = 0.;
859  sCoefficient[1] = 0.;
860  sCoefficient[2] = 0.;
861  }
862 
863  else if (particleDefinition == instance->GetIon("alpha+"))
864  {
865  slaterEffectiveCharge[0] = 2.0;
866  // The following values are provided by M. Dingfelder (priv. comm)
867  slaterEffectiveCharge[1] = 2.0;
868  slaterEffectiveCharge[2] = 2.0;
869  //
870  sCoefficient[0] = 0.7;
871  sCoefficient[1] = 0.15;
872  sCoefficient[2] = 0.15;
873  }
874 
875  else if (particleDefinition == instance->GetIon("helium"))
876  {
877  slaterEffectiveCharge[0] = 1.7;
878  slaterEffectiveCharge[1] = 1.15;
879  slaterEffectiveCharge[2] = 1.15;
880  sCoefficient[0] = 0.5;
881  sCoefficient[1] = 0.25;
882  sCoefficient[2] = 0.25;
883  }
884 
885 // if ( particleDefinition == instance->GetIon("helium")
886 // || particleDefinition == instance->GetIon("alpha+")
887 // || particleDefinition == instance->GetIon("alpha++")
888 // )
889  if (isHelium)
890  {
891  sigma = Gj[j] * (S / Bj[ionizationLevelIndex])
892  * ((F1 + w * F2)
893  / (std::pow((1. + w), 3)
894  * (1. + std::exp(alphaConst * (w - wc) / v))));
895 
896  if (j == 4)
897  sigma = Gj[j]
898  * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
899  * ((F1 + w * F2)
900  / (std::pow((1. + w), 3)
901  * (1. + std::exp(alphaConst * (w - wc) / v))));
902 
903  G4double zEff = particleDefinition->GetPDGCharge() / eplus
904  + particleDefinition->GetLeptonNumber();
905 
906  zEff -= (sCoefficient[0]
907  * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.)
908  + sCoefficient[1]
909  * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.)
910  + sCoefficient[2]
911  * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.));
912 
913  return zEff * zEff * sigma;
914  }
915 
916  return 0;
917 }
const double C2
const double C1
double S(double temp)
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4double S_2s(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
int G4int
Definition: G4Types.hh:78
Char_t n[5]
float Bohr_radius
Definition: hepunit.py:290
bool G4bool
Definition: G4Types.hh:79
float proton_mass_c2
Definition: hepunit.py:275
float electron_mass_c2
Definition: hepunit.py:274
static G4DNAGenericIonsManager * Instance(void)
static const double eV
Definition: G4SIunits.hh:212
static const double pi
Definition: G4SIunits.hh:74
G4double S_1s(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
G4DNAWaterIonisationStructure waterStructure
G4double CorrectionFactor(G4ParticleDefinition *particleDefinition, G4double k)
static MCTruthManager * instance
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:196
G4double S_2p(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
G4double GetPDGCharge() const
G4ParticleDefinition * GetIon(const G4String &name)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Initialise()

void G4DNARuddIonisationModel::Initialise ( const G4ParticleDefinition particle,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 113 of file G4DNARuddIonisationModel.cc.

115 {
116 
117  if (verboseLevel > 3)
118  {
119  G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
120  }
121 
122  // Energy limits
123 
124  G4String fileProton("dna/sigma_ionisation_p_rudd");
125  G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
126  G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
127  G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
128  G4String fileHelium("dna/sigma_ionisation_he_rudd");
129 
133  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
134  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
135  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
136  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
137 
139  G4String hydrogen;
140  G4String alphaPlusPlus;
141  G4String alphaPlus;
142  G4String helium;
143 
144  G4double scaleFactor = 1 * m*m;
145 
146  // LIMITS AND DATA
147 
148  // ********************************************************
149 
150  proton = protonDef->GetParticleName();
151  tableFile[proton] = fileProton;
152 
154  highEnergyLimit[proton] = 500. * keV;
155 
156  // Cross section
157 
159  eV,
160  scaleFactor );
161  tableProton->LoadData(fileProton);
162  tableData[proton] = tableProton;
163 
164  // ********************************************************
165 
166  hydrogen = hydrogenDef->GetParticleName();
167  tableFile[hydrogen] = fileHydrogen;
168 
170  highEnergyLimit[hydrogen] = 100. * MeV;
171 
172  // Cross section
173 
175  eV,
176  scaleFactor );
177  tableHydrogen->LoadData(fileHydrogen);
178 
179  tableData[hydrogen] = tableHydrogen;
180 
181  // ********************************************************
182 
183  alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
184  tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
185 
186  lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
187  highEnergyLimit[alphaPlusPlus] = 400. * MeV;
188 
189  // Cross section
190 
192  eV,
193  scaleFactor );
194  tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
195 
196  tableData[alphaPlusPlus] = tableAlphaPlusPlus;
197 
198  // ********************************************************
199 
200  alphaPlus = alphaPlusDef->GetParticleName();
201  tableFile[alphaPlus] = fileAlphaPlus;
202 
203  lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
204  highEnergyLimit[alphaPlus] = 400. * MeV;
205 
206  // Cross section
207 
209  eV,
210  scaleFactor );
211  tableAlphaPlus->LoadData(fileAlphaPlus);
212  tableData[alphaPlus] = tableAlphaPlus;
213 
214  // ********************************************************
215 
216  helium = heliumDef->GetParticleName();
217  tableFile[helium] = fileHelium;
218 
220  highEnergyLimit[helium] = 400. * MeV;
221 
222  // Cross section
223 
225  eV,
226  scaleFactor );
227  tableHelium->LoadData(fileHelium);
228  tableData[helium] = tableHelium;
229 
230  //
231 
232  if (particle==protonDef)
233  {
236  }
237 
238  if (particle==hydrogenDef)
239  {
242  }
243 
244  if (particle==heliumDef)
245  {
248  }
249 
250  if (particle==alphaPlusDef)
251  {
252  SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
254  }
255 
256  if (particle==alphaPlusPlusDef)
257  {
258  SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
259  SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
260  }
261 
262  if( verboseLevel>0 )
263  {
264  G4cout << "Rudd ionisation model is initialized " << G4endl
265  << "Energy range: "
266  << LowEnergyLimit() / eV << " eV - "
267  << HighEnergyLimit() / keV << " keV for "
268  << particle->GetParticleName()
269  << G4endl;
270  }
271 
272  // Initialize water density pointer
274 
275  //
276 
278 
279  if (isInitialised)
280  { return;}
282  isInitialised = true;
283 
284 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
static const double MeV
Definition: G4SIunits.hh:211
static G4LossTableManager * Instance()
G4ParticleChangeForGamma * fParticleChangeForGamma
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:604
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
virtual G4bool LoadData(const G4String &argFileName)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
const std::vector< G4double > * fpWaterDensity
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
static G4DNAGenericIonsManager * Instance(void)
static const double eV
Definition: G4SIunits.hh:212
static G4DNAMolecularMaterial * Instance()
#define G4endl
Definition: G4ios.hh:61
static const double m
Definition: G4SIunits.hh:128
static MCTruthManager * instance
G4VAtomDeexcitation * AtomDeexcitation()
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
G4VAtomDeexcitation * fAtomDeexcitation
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
G4ParticleDefinition * GetIon(const G4String &name)
Here is the call graph for this function:

◆ operator=()

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

◆ PartialCrossSection()

G4double G4DNARuddIonisationModel::PartialCrossSection ( const G4Track &  track)
private

Definition at line 1081 of file G4DNARuddIonisationModel.cc.

1082 {
1083  G4double sigma = 0.;
1084 
1085  const G4DynamicParticle* particle = track.GetDynamicParticle();
1086  G4double k = particle->GetKineticEnergy();
1087 
1088  G4double lowLim = 0;
1089  G4double highLim = 0;
1090 
1091  const G4String& particleName = particle->GetDefinition()->GetParticleName();
1092 
1093  std::map<G4String, G4double, std::less<G4String> >::iterator pos1;
1094  pos1 = lowEnergyLimit.find(particleName);
1095 
1096  if (pos1 != lowEnergyLimit.end())
1097  {
1098  lowLim = pos1->second;
1099  }
1100 
1101  std::map<G4String, G4double, std::less<G4String> >::iterator pos2;
1102  pos2 = highEnergyLimit.find(particleName);
1103 
1104  if (pos2 != highEnergyLimit.end())
1105  {
1106  highLim = pos2->second;
1107  }
1108 
1109  if (k >= lowLim && k <= highLim)
1110  {
1111  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1112  pos = tableData.find(particleName);
1113 
1114  if (pos != tableData.end())
1115  {
1116  G4DNACrossSectionDataSet* table = pos->second;
1117  if (table != 0)
1118  {
1119  sigma = table->FindValue(k);
1120  }
1121  } else
1122  {
1123  G4Exception("G4DNARuddIonisationModel::PartialCrossSection",
1124  "em0002",
1126  "Model not applicable to particle type.");
1127  }
1128  }
1129 
1130  return sigma;
1131 }
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
G4double GetKineticEnergy() const
const G4String & GetParticleName() const
virtual G4double FindValue(G4double e, G4int componentId=0) const
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ParticleDefinition * GetDefinition() const
double G4double
Definition: G4Types.hh:76
static const G4double pos
Here is the call graph for this function:

◆ R()

G4double G4DNARuddIonisationModel::R ( G4double  t,
G4double  energyTransferred,
G4double  slaterEffectiveChg,
G4double  shellNumber 
)
private

Definition at line 973 of file G4DNARuddIonisationModel.cc.

977 {
978  // tElectron = m_electron / m_alpha * t
979  // Dingfelder, in Chattanooga 2005 proceedings, p 4
980 
981  G4double tElectron = 0.511 / 3728. * t;
982  // The following values are provided by M. Dingfelder (priv. comm)
983  G4double H = 2. * 13.60569172 * eV;
984  G4double value = std::sqrt(2. * tElectron / H) / (energyTransferred / H)
985  * (slaterEffectiveChg / shellNumber);
986 
987  return value;
988 }
static const double eV
Definition: G4SIunits.hh:212
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ RandomizeEjectedElectronEnergy()

G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy ( G4ParticleDefinition particleDefinition,
G4double  incomingParticleEnergy,
G4int  shell 
)
private

Definition at line 601 of file G4DNARuddIonisationModel.cc.

604 {
605  G4double maximumKineticEnergyTransfer = 0.;
606 
609 
610  if (particleDefinition == G4Proton::ProtonDefinition()
611  || particleDefinition == instance->GetIon("hydrogen"))
612  {
613  maximumKineticEnergyTransfer = 4. * (electron_mass_c2 / proton_mass_c2) * k;
614  }
615 
616  else if (particleDefinition == instance->GetIon("helium")
617  || particleDefinition == instance->GetIon("alpha+")
618  || particleDefinition == instance->GetIon("alpha++"))
619  {
620  maximumKineticEnergyTransfer = 4. * (0.511 / 3728) * k;
621  }
622 
623  G4double crossSectionMaximum = 0.;
624 
625  for (G4double value = waterStructure.IonisationEnergy(shell);
626  value <= 5. * waterStructure.IonisationEnergy(shell) && k >= value;
627  value += 0.1 * eV)
628  {
629  G4double differentialCrossSection =
630  DifferentialCrossSection(particleDefinition, k, value, shell);
631  if (differentialCrossSection >= crossSectionMaximum)
632  crossSectionMaximum = differentialCrossSection;
633  }
634 
635  G4double secElecKinetic = 0.;
636 
637  do
638  {
639  secElecKinetic = G4UniformRand()* maximumKineticEnergyTransfer;
640  }while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
641  k,
642  secElecKinetic+waterStructure.IonisationEnergy(shell),
643  shell));
644 
645  return (secElecKinetic);
646 }
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4double DifferentialCrossSection(G4ParticleDefinition *particleDefinition, G4double k, G4double energyTransfer, G4int shell)
#define G4UniformRand()
Definition: Randomize.hh:97
float proton_mass_c2
Definition: hepunit.py:275
float electron_mass_c2
Definition: hepunit.py:274
static G4DNAGenericIonsManager * Instance(void)
static const double eV
Definition: G4SIunits.hh:212
G4DNAWaterIonisationStructure waterStructure
static MCTruthManager * instance
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition * GetIon(const G4String &name)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ RandomSelect()

G4int G4DNARuddIonisationModel::RandomSelect ( G4double  energy,
const G4String particle 
)
private

Definition at line 1014 of file G4DNARuddIonisationModel.cc.

1016 {
1017 
1018  // BEGIN PART 1/2 OF ELECTRON CORRECTION
1019 
1020  // add ONE or TWO electron-water ionisation for alpha+ and helium
1021 
1022  G4int level = 0;
1023 
1024  // Retrieve data table corresponding to the current particle type
1025 
1026  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1027  pos = tableData.find(particle);
1028 
1029  if (pos != tableData.end())
1030  {
1031  G4DNACrossSectionDataSet* table = pos->second;
1032 
1033  if (table != 0)
1034  {
1035  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1036 
1037  const size_t n(table->NumberOfComponents());
1038  size_t i(n);
1039  G4double value = 0.;
1040 
1041  while (i > 0)
1042  {
1043  i--;
1044  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1045  value += valuesBuffer[i];
1046  }
1047 
1048  value *= G4UniformRand();
1049 
1050  i = n;
1051 
1052  while (i > 0)
1053  {
1054  i--;
1055 
1056  if (valuesBuffer[i] > value)
1057  {
1058  delete[] valuesBuffer;
1059  return i;
1060  }
1061  value -= valuesBuffer[i];
1062  }
1063 
1064  if (valuesBuffer)
1065  delete[] valuesBuffer;
1066 
1067  }
1068  } else
1069  {
1070  G4Exception("G4DNARuddIonisationModel::RandomSelect",
1071  "em0002",
1073  "Model not applicable to particle type.");
1074  }
1075 
1076  return level;
1077 }
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
int G4int
Definition: G4Types.hh:78
Char_t n[5]
#define G4UniformRand()
Definition: Randomize.hh:97
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double G4double
Definition: G4Types.hh:76
static const G4double pos
virtual size_t NumberOfComponents(void) const
Here is the call graph for this function:
Here is the caller graph for this function:

◆ S_1s()

G4double G4DNARuddIonisationModel::S_1s ( G4double  t,
G4double  energyTransferred,
G4double  slaterEffectiveChg,
G4double  shellNumber 
)
private

Definition at line 921 of file G4DNARuddIonisationModel.cc.

925 {
926  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
927  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
928 
929  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
930  G4double value = 1. - std::exp(-2 * r) * ((2. * r + 2.) * r + 1.);
931 
932  return value;
933 }
G4double R(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ S_2p()

G4double G4DNARuddIonisationModel::S_2p ( G4double  t,
G4double  energyTransferred,
G4double  slaterEffectiveChg,
G4double  shellNumber 
)
private

Definition at line 955 of file G4DNARuddIonisationModel.cc.

959 {
960  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
961  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
962 
963  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
964  G4double value = 1.
965  - std::exp(-2 * r)
966  * ((((2. / 3. * r + 4. / 3.) * r + 2.) * r + 2.) * r + 1.);
967 
968  return value;
969 }
G4double R(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ S_2s()

G4double G4DNARuddIonisationModel::S_2s ( G4double  t,
G4double  energyTransferred,
G4double  slaterEffectiveChg,
G4double  shellNumber 
)
private

Definition at line 937 of file G4DNARuddIonisationModel.cc.

941 {
942  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
943  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
944 
945  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
946  G4double value = 1.
947  - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
948 
949  return value;
950 
951 }
G4double R(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleSecondaries()

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

Implements G4VEmModel.

Definition at line 417 of file G4DNARuddIonisationModel.cc.

422 {
423  if (verboseLevel > 3)
424  {
425  G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel"
426  << G4endl;
427  }
428 
429  G4double lowLim = 0;
430  G4double highLim = 0;
431 
434 
435  if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
436  || particle->GetDefinition() == instance->GetIon("hydrogen")
437  )
438 
439  lowLim = killBelowEnergyForZ1;
440 
441  if ( particle->GetDefinition() == instance->GetIon("alpha++")
442  || particle->GetDefinition() == instance->GetIon("alpha+")
443  || particle->GetDefinition() == instance->GetIon("helium")
444  )
445 
446  lowLim = killBelowEnergyForZ2;
447 
448  G4double k = particle->GetKineticEnergy();
449 
450  const G4String& particleName = particle->GetDefinition()->GetParticleName();
451 
452  // SI - the following is useless since lowLim is already defined
453  /*
454  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
455  pos1 = lowEnergyLimit.find(particleName);
456 
457  if (pos1 != lowEnergyLimit.end())
458  {
459  lowLim = pos1->second;
460  }
461  */
462 
463  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
464  pos2 = highEnergyLimit.find(particleName);
465 
466  if (pos2 != highEnergyLimit.end())
467  {
468  highLim = pos2->second;
469  }
470 
471  if (k >= lowLim && k <= highLim)
472  {
473  G4ParticleDefinition* definition = particle->GetDefinition();
474  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
475  /*
476  G4double particleMass = definition->GetPDGMass();
477  G4double totalEnergy = k + particleMass;
478  G4double pSquare = k*(totalEnergy+particleMass);
479  G4double totalMomentum = std::sqrt(pSquare);
480  */
481 
482  G4int ionizationShell = RandomSelect(k,particleName);
483 
484  // sample deexcitation
485  // here we assume that H_{2}O electronic levels are the same of Oxigen.
486  // this can be considered true with a rough 10% error in energy on K-shell,
487 
488  G4int secNumberInit = 0;// need to know at a certain point the enrgy of secondaries
489  G4int secNumberFinal = 0;// So I'll make the diference and then sum the energies
490 
491  G4double bindingEnergy = 0;
492  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
493 
494  //SI: additional protection if tcs interpolation method is modified
495  if (k<bindingEnergy) return;
496  //
497 
498  G4int Z = 8;
500  {
502 
503  if (ionizationShell <5 && ionizationShell >1)
504  {
505  as = G4AtomicShellEnumerator(4-ionizationShell);
506  }
507  else if (ionizationShell <2)
508  {
509  as = G4AtomicShellEnumerator(3);
510  }
511 
512  // DEBUG
513  // if (ionizationShell == 4) {
514  //
515  // G4cout << "Z: " << Z << " as: " << as
516  // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
517  // G4cout << "Press <Enter> key to continue..." << G4endl;
518  // G4cin.ignore();
519  // }
520 
521  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
522  secNumberInit = fvect->size();
523  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
524  secNumberFinal = fvect->size();
525  }
526 
527  G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
528 
529  G4ThreeVector deltaDirection =
530  GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
531  Z, ionizationShell,
532  couple->GetMaterial());
533 
534  // Ignored for ions on electrons
535  /*
536  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
537 
538  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
539  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
540  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
541  G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
542  finalPx /= finalMomentum;
543  finalPy /= finalMomentum;
544  finalPz /= finalMomentum;
545 
546  G4ThreeVector direction;
547  direction.set(finalPx,finalPy,finalPz);
548 
549  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
550  */
551  fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
552 
553  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
554  G4double deexSecEnergy = 0;
555  for (G4int j=secNumberInit; j < secNumberFinal; j++)
556  {
557 
558  deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
559 
560  }
561 
562  if (!statCode)
563  {
564  fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
565  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
566  }
567  else
568  {
569  fParticleChangeForGamma->SetProposedKineticEnergy(k);
570  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
571  }
572 
573  // debug
574  // k-scatteredEnergy-secondaryKinetic-deexSecEnergy = k-(k-bindingEnergy-secondaryKinetic)-secondaryKinetic-deexSecEnergy =
575  // = k-k+bindingEnergy+secondaryKinetic-secondaryKinetic-deexSecEnergy=
576  // = bindingEnergy-deexSecEnergy
577  // SO deexSecEnergy=0 => LocalEnergyDeposit = bindingEnergy
578 
579  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
580  fvect->push_back(dp);
581 
582  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
584  ionizationShell,
585  theIncomingTrack);
586  }
587 
588  // SI - not useful since low energy of model is 0 eV
589 
590  if (k < lowLim)
591  {
592  fParticleChangeForGamma->SetProposedKineticEnergy(0.);
593  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
594  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
595  }
596 
597 }
const G4Material * GetMaterial() const
G4ParticleChangeForGamma * fParticleChangeForGamma
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:617
int G4int
Definition: G4Types.hh:78
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4double GetKineticEnergy() const
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
Float_t Z
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
G4int RandomSelect(G4double energy, const G4String &particle)
static G4DNAGenericIonsManager * Instance(void)
static G4DNAChemistryManager * Instance()
const G4ThreeVector & GetMomentumDirection() const
void GenerateParticles(std::vector< G4DynamicParticle *> *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
G4DNAWaterIonisationStructure waterStructure
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *particleDefinition, G4double incomingParticleEnergy, G4int shell)
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleDefinition * GetDefinition() const
#define G4endl
Definition: G4ios.hh:61
static MCTruthManager * instance
double G4double
Definition: G4Types.hh:76
G4VAtomDeexcitation * fAtomDeexcitation
G4AtomicShellEnumerator
G4ParticleDefinition * GetIon(const G4String &name)
Here is the call graph for this function:

◆ SelectStationary()

void G4DNARuddIonisationModel::SelectStationary ( G4bool  input)
inline

Definition at line 163 of file G4DNARuddIonisationModel.hh.

164 {
165  statCode = input;
166 }

◆ Sum()

G4double G4DNARuddIonisationModel::Sum ( G4double  energy,
const G4String particle 
)
private

Definition at line 1135 of file G4DNARuddIonisationModel.cc.

1137 {
1138  return 0;
1139 }

Member Data Documentation

◆ fAtomDeexcitation

G4VAtomDeexcitation* G4DNARuddIonisationModel::fAtomDeexcitation
private

Definition at line 84 of file G4DNARuddIonisationModel.hh.

◆ fParticleChangeForGamma

G4ParticleChangeForGamma* G4DNARuddIonisationModel::fParticleChangeForGamma
protected

Definition at line 74 of file G4DNARuddIonisationModel.hh.

◆ fpWaterDensity

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

Definition at line 81 of file G4DNARuddIonisationModel.hh.

◆ highEnergyLimit

std::map<G4String,G4double,std::less<G4String> > G4DNARuddIonisationModel::highEnergyLimit
private

Definition at line 87 of file G4DNARuddIonisationModel.hh.

◆ isInitialised

G4bool G4DNARuddIonisationModel::isInitialised
private

Definition at line 96 of file G4DNARuddIonisationModel.hh.

◆ killBelowEnergyForZ1

G4double G4DNARuddIonisationModel::killBelowEnergyForZ1
private

Definition at line 93 of file G4DNARuddIonisationModel.hh.

◆ killBelowEnergyForZ2

G4double G4DNARuddIonisationModel::killBelowEnergyForZ2
private

Definition at line 94 of file G4DNARuddIonisationModel.hh.

◆ lowEnergyLimit

std::map<G4String,G4double,std::less<G4String> > G4DNARuddIonisationModel::lowEnergyLimit
private

Definition at line 86 of file G4DNARuddIonisationModel.hh.

◆ lowEnergyLimitForZ1

G4double G4DNARuddIonisationModel::lowEnergyLimitForZ1
private

Definition at line 89 of file G4DNARuddIonisationModel.hh.

◆ lowEnergyLimitForZ2

G4double G4DNARuddIonisationModel::lowEnergyLimitForZ2
private

Definition at line 90 of file G4DNARuddIonisationModel.hh.

◆ lowEnergyLimitOfModelForZ1

G4double G4DNARuddIonisationModel::lowEnergyLimitOfModelForZ1
private

Definition at line 91 of file G4DNARuddIonisationModel.hh.

◆ lowEnergyLimitOfModelForZ2

G4double G4DNARuddIonisationModel::lowEnergyLimitOfModelForZ2
private

Definition at line 92 of file G4DNARuddIonisationModel.hh.

◆ sCoefficient

G4double G4DNARuddIonisationModel::sCoefficient[3]
private

Definition at line 144 of file G4DNARuddIonisationModel.hh.

◆ slaterEffectiveCharge

G4double G4DNARuddIonisationModel::slaterEffectiveCharge[3]
private

Definition at line 143 of file G4DNARuddIonisationModel.hh.

◆ statCode

G4bool G4DNARuddIonisationModel::statCode
private

Definition at line 78 of file G4DNARuddIonisationModel.hh.

◆ tableData

MapData G4DNARuddIonisationModel::tableData
private

Definition at line 105 of file G4DNARuddIonisationModel.hh.

◆ tableFile

MapFile G4DNARuddIonisationModel::tableFile
private

Definition at line 102 of file G4DNARuddIonisationModel.hh.

◆ verboseLevel

G4int G4DNARuddIonisationModel::verboseLevel
private

Definition at line 97 of file G4DNARuddIonisationModel.hh.

◆ waterStructure

G4DNAWaterIonisationStructure G4DNARuddIonisationModel::waterStructure
private

Definition at line 109 of file G4DNARuddIonisationModel.hh.


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