Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4IonParametrisedLossModel Class Reference

#include <G4IonParametrisedLossModel.hh>

Inheritance diagram for G4IonParametrisedLossModel:
Collaboration diagram for G4IonParametrisedLossModel:

Public Member Functions

 G4IonParametrisedLossModel (const G4ParticleDefinition *particle=0, const G4String &name="ParamICRU73")
 
virtual ~G4IonParametrisedLossModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double, G4double, G4double)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double, G4double)
 
G4double ComputeLossForStep (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double, G4double)
 
G4double DeltaRayMeanEnergyTransferRate (const G4Material *, const G4ParticleDefinition *, G4double, G4double)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &, G4double &, G4double)
 
G4bool AddDEDXTable (const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
 
G4bool RemoveDEDXTable (const G4String &name)
 
void DeactivateICRU73Scaling ()
 
LossTableList::iterator IsApplicable (const G4ParticleDefinition *, const G4Material *)
 
void PrintDEDXTable (const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
 
void PrintDEDXTableHandlers (const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
 
void SetEnergyLossLimit (G4double ionEnergyLossLimit)
 
- 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 GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
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 void StartTracking (G4Track *)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
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 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=nullptr)
 
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 SetFluctuationFlag (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 Member Functions

virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double)
 
- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 

Additional Inherited Members

- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChangepParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 93 of file G4IonParametrisedLossModel.hh.

Constructor & Destructor Documentation

G4IonParametrisedLossModel::G4IonParametrisedLossModel ( const G4ParticleDefinition particle = 0,
const G4String name = "ParamICRU73" 
)

Definition at line 105 of file G4IonParametrisedLossModel.cc.

108  : G4VEmModel(nam),
109  braggIonModel(0),
110  betheBlochModel(0),
111  nmbBins(90),
112  nmbSubBins(100),
113  particleChangeLoss(0),
114  corrFactor(1.0),
115  energyLossLimit(0.01),
116  cutEnergies(0)
117 {
118  genericIon = G4GenericIon::Definition();
119  genericIonPDGMass = genericIon -> GetPDGMass();
120  corrections = G4LossTableManager::Instance() -> EmCorrections();
121 
122  // The upper limit of the current model is set to 100 TeV
123  SetHighEnergyLimit(100.0 * TeV);
124 
125  // The Bragg ion and Bethe Bloch models are instantiated
126  braggIonModel = new G4BraggIonModel();
127  betheBlochModel = new G4BetheBlochModel();
128 
129  // By default ICRU 73 stopping power tables are loaded:
130  AddDEDXTable("ICRU73",
131  new G4IonStoppingData("ion_stopping_data/icru73"),
132  new G4IonDEDXScalingICRU73());
133 
134  // The boundaries for the range tables are set
135  lowerEnergyEdgeIntegr = 0.025 * MeV;
136  upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit();
137 
138  // Cache parameters are set
139  cacheParticle = 0;
140  cacheMass = 0;
141  cacheElecMassRatio = 0;
142  cacheChargeSquare = 0;
143 
144  // Cache parameters are set
145  rangeCacheParticle = 0;
146  rangeCacheMatCutsCouple = 0;
147  rangeCacheEnergyRange = 0;
148  rangeCacheRangeEnergy = 0;
149 
150  // Cache parameters are set
151  dedxCacheParticle = 0;
152  dedxCacheMaterial = 0;
153  dedxCacheEnergyCut = 0;
154  dedxCacheIter = lossTableList.end();
155  dedxCacheTransitionEnergy = 0.0;
156  dedxCacheTransitionFactor = 0.0;
157  dedxCacheGenIonMassRatio = 0.0;
158 
159  // default generator
161 }
static G4LossTableManager * Instance()
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
static G4GenericIon * Definition()
Definition: G4GenericIon.cc:49
G4bool AddDEDXTable(const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:68
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
static constexpr double TeV
Definition: G4SIunits.hh:218
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:626
static constexpr double MeV
Definition: G4SIunits.hh:214

Here is the call graph for this function:

G4IonParametrisedLossModel::~G4IonParametrisedLossModel ( )
virtual

Definition at line 165 of file G4IonParametrisedLossModel.cc.

165  {
166 
167  // dE/dx table objects are deleted and the container is cleared
168  LossTableList::iterator iterTables = lossTableList.begin();
169  LossTableList::iterator iterTables_end = lossTableList.end();
170 
171  for(;iterTables != iterTables_end; ++iterTables) { delete *iterTables; }
172  lossTableList.clear();
173 
174  // range table
175  RangeEnergyTable::iterator itr = r.begin();
176  RangeEnergyTable::iterator itr_end = r.end();
177  for(;itr != itr_end; ++itr) { delete itr->second; }
178  r.clear();
179 
180  // inverse range
181  EnergyRangeTable::iterator ite = E.begin();
182  EnergyRangeTable::iterator ite_end = E.end();
183  for(;ite != ite_end; ++ite) { delete ite->second; }
184  E.clear();
185 
186 }

Member Function Documentation

G4bool G4IonParametrisedLossModel::AddDEDXTable ( const G4String name,
G4VIonDEDXTable table,
G4VIonDEDXScalingAlgorithm algorithm = 0 
)

Definition at line 1243 of file G4IonParametrisedLossModel.cc.

1246  {
1247 
1248  if(table == 0) {
1249  G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1250  << " add table: Invalid pointer."
1251  << G4endl;
1252 
1253  return false;
1254  }
1255 
1256  // Checking uniqueness of name
1257  LossTableList::iterator iter = lossTableList.begin();
1258  LossTableList::iterator iter_end = lossTableList.end();
1259 
1260  for(;iter != iter_end; iter++) {
1261  G4String tableName = (*iter) -> GetName();
1262 
1263  if(tableName == nam) {
1264  G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1265  << " add table: Name already exists."
1266  << G4endl;
1267 
1268  return false;
1269  }
1270  }
1271 
1272  G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm;
1273  if(scalingAlgorithm == 0)
1274  scalingAlgorithm = new G4VIonDEDXScalingAlgorithm;
1275 
1276  G4IonDEDXHandler* handler =
1277  new G4IonDEDXHandler(table, scalingAlgorithm, nam);
1278 
1279  lossTableList.push_front(handler);
1280 
1281  return true;
1282 }
#define G4endl
Definition: G4ios.hh:61
const G4String & GetName() const
Definition: G4VEmModel.hh:802
G4GLOB_DLL std::ostream G4cerr

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4IonParametrisedLossModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition particle,
G4double  kineticEnergy,
G4double  atomicNumber,
G4double  ,
G4double  cutEnergy,
G4double  maxKinEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 374 of file G4IonParametrisedLossModel.cc.

380  {
381 
382  // ############## Cross section per atom ################################
383  // Function computes ionization cross section per atom
384  //
385  // See Geant4 physics reference manual (version 9.1), section 9.1.3
386  //
387  // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
388  // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
389  //
390  // (Implementation adapted from G4BraggIonModel)
391 
392  G4double crosssection = 0.0;
393  G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
394  G4double maxEnergy = std::min(tmax, maxKinEnergy);
395 
396  if(cutEnergy < tmax) {
397 
398  G4double energy = kineticEnergy + cacheMass;
399  G4double betaSquared = kineticEnergy *
400  (energy + cacheMass) / (energy * energy);
401 
402  crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
403  betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
404 
405  crosssection *= twopi_mc2_rcl2 * cacheChargeSquare / betaSquared;
406  }
407 
408 #ifdef PRINT_DEBUG_CS
409  G4cout << "########################################################"
410  << G4endl
411  << "# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
412  << G4endl
413  << "# particle =" << particle -> GetParticleName()
414  << G4endl
415  << "# cut(MeV) = " << cutEnergy/MeV
416  << G4endl;
417 
418  G4cout << "#"
419  << std::setw(13) << std::right << "E(MeV)"
420  << std::setw(14) << "CS(um)"
421  << std::setw(14) << "E_max_sec(MeV)"
422  << G4endl
423  << "# ------------------------------------------------------"
424  << G4endl;
425 
426  G4cout << std::setw(14) << std::right << kineticEnergy / MeV
427  << std::setw(14) << crosssection / (um * um)
428  << std::setw(14) << tmax / MeV
429  << G4endl;
430 #endif
431 
432  crosssection *= atomicNumber;
433 
434  return crosssection;
435 }
G4GLOB_DLL std::ostream G4cout
static constexpr double um
Definition: G4SIunits.hh:113
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double)
G4double energy(const ThreeVector &p, const G4double m)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4IonParametrisedLossModel::ComputeDEDXPerVolume ( const G4Material material,
const G4ParticleDefinition particle,
G4double  kineticEnergy,
G4double  cutEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 458 of file G4IonParametrisedLossModel.cc.

462  {
463 
464  // ############## dE/dx ##################################################
465  // Function computes dE/dx values, where following rules are adopted:
466  // A. If the ion-material pair is covered by any native ion data
467  // parameterisation, then:
468  // * This parameterization is used for energies below a given energy
469  // limit,
470  // * whereas above the limit the Bethe-Bloch model is applied, in
471  // combination with an effective charge estimate and high order
472  // correction terms.
473  // A smoothing procedure is applied to dE/dx values computed with
474  // the second approach. The smoothing factor is based on the dE/dx
475  // values of both approaches at the transition energy (high order
476  // correction terms are included in the calculation of the transition
477  // factor).
478  // B. If the particle is a generic ion, the BraggIon and Bethe-Bloch
479  // models are used and a smoothing procedure is applied to values
480  // obtained with the second approach.
481  // C. If the ion-material is not covered by any ion data parameterization
482  // then:
483  // * The BraggIon model is used for energies below a given energy
484  // limit,
485  // * whereas above the limit the Bethe-Bloch model is applied, in
486  // combination with an effective charge estimate and high order
487  // correction terms.
488  // Also in this case, a smoothing procedure is applied to dE/dx values
489  // computed with the second model.
490 
491  G4double dEdx = 0.0;
492  UpdateDEDXCache(particle, material, cutEnergy);
493 
494  LossTableList::iterator iter = dedxCacheIter;
495 
496  if(iter != lossTableList.end()) {
497 
498  G4double transitionEnergy = dedxCacheTransitionEnergy;
499 
500  if(transitionEnergy > kineticEnergy) {
501 
502  dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
503 
504  G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material,
505  particle,
506  kineticEnergy,
507  cutEnergy);
508  dEdx -= dEdxDeltaRays;
509  }
510  else {
511  G4double massRatio = dedxCacheGenIonMassRatio;
512 
513  G4double chargeSquare =
514  GetChargeSquareRatio(particle, material, kineticEnergy);
515 
516  G4double scaledKineticEnergy = kineticEnergy * massRatio;
517  G4double scaledTransitionEnergy = transitionEnergy * massRatio;
518 
519  G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
520 
521  if(scaledTransitionEnergy >= lowEnergyLimit) {
522 
523  dEdx = betheBlochModel -> ComputeDEDXPerVolume(
524  material, genericIon,
525  scaledKineticEnergy, cutEnergy);
526 
527  dEdx *= chargeSquare;
528 
529  dEdx += corrections -> ComputeIonCorrections(particle,
530  material, kineticEnergy);
531 
532  G4double factor = 1.0 + dedxCacheTransitionFactor /
533  kineticEnergy;
534 
535  dEdx *= factor;
536  }
537  }
538  }
539  else {
540  G4double massRatio = 1.0;
541  G4double chargeSquare = 1.0;
542 
543  if(particle != genericIon) {
544 
545  chargeSquare = GetChargeSquareRatio(particle, material, kineticEnergy);
546  massRatio = genericIonPDGMass / particle -> GetPDGMass();
547  }
548 
549  G4double scaledKineticEnergy = kineticEnergy * massRatio;
550 
551  G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
552  if(scaledKineticEnergy < lowEnergyLimit) {
553  dEdx = braggIonModel -> ComputeDEDXPerVolume(
554  material, genericIon,
555  scaledKineticEnergy, cutEnergy);
556 
557  dEdx *= chargeSquare;
558  }
559  else {
560  G4double dEdxLimitParam = braggIonModel -> ComputeDEDXPerVolume(
561  material, genericIon,
562  lowEnergyLimit, cutEnergy);
563 
564  G4double dEdxLimitBetheBloch = betheBlochModel -> ComputeDEDXPerVolume(
565  material, genericIon,
566  lowEnergyLimit, cutEnergy);
567 
568  if(particle != genericIon) {
569  G4double chargeSquareLowEnergyLimit =
570  GetChargeSquareRatio(particle, material,
571  lowEnergyLimit / massRatio);
572 
573  dEdxLimitParam *= chargeSquareLowEnergyLimit;
574  dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
575 
576  dEdxLimitBetheBloch +=
577  corrections -> ComputeIonCorrections(particle,
578  material, lowEnergyLimit / massRatio);
579  }
580 
581  G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
582  * lowEnergyLimit / scaledKineticEnergy);
583 
584  dEdx = betheBlochModel -> ComputeDEDXPerVolume(
585  material, genericIon,
586  scaledKineticEnergy, cutEnergy);
587 
588  dEdx *= chargeSquare;
589 
590  if(particle != genericIon) {
591  dEdx += corrections -> ComputeIonCorrections(particle,
592  material, kineticEnergy);
593  }
594 
595  dEdx *= factor;
596  }
597 
598  }
599 
600  if (dEdx < 0.0) dEdx = 0.0;
601 
602  return dEdx;
603 }
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
G4double DeltaRayMeanEnergyTransferRate(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4IonParametrisedLossModel::ComputeLossForStep ( const G4MaterialCutsCouple matCutsCouple,
const G4ParticleDefinition particle,
G4double  kineticEnergy,
G4double  stepLength 
)

Definition at line 1184 of file G4IonParametrisedLossModel.cc.

1188  {
1189 
1190  G4double loss = 0.0;
1191 
1192  UpdateRangeCache(particle, matCutsCouple);
1193 
1194  G4PhysicsVector* energyRange = rangeCacheEnergyRange;
1195  G4PhysicsVector* rangeEnergy = rangeCacheRangeEnergy;
1196 
1197  if(energyRange != 0 && rangeEnergy != 0) {
1198 
1199  G4double lowerEnEdge = energyRange -> Energy( 0 );
1200  G4double lowerRangeEdge = rangeEnergy -> Energy( 0 );
1201 
1202  // Computing range for pre-step kinetic energy:
1203  G4double range = energyRange -> Value(kineticEnergy);
1204 
1205  // Energy below vector boundary:
1206  if(kineticEnergy < lowerEnEdge) {
1207 
1208  range = energyRange -> Value(lowerEnEdge);
1209  range *= std::sqrt(kineticEnergy / lowerEnEdge);
1210  }
1211 
1212 #ifdef PRINT_DEBUG
1213  G4cout << "G4IonParametrisedLossModel::ComputeLossForStep() range = "
1214  << range / mm << " mm, step = " << stepLength / mm << " mm"
1215  << G4endl;
1216 #endif
1217 
1218  // Remaining range:
1219  G4double remRange = range - stepLength;
1220 
1221  // If range is smaller than step length, the loss is set to kinetic
1222  // energy
1223  if(remRange < 0.0) loss = kineticEnergy;
1224  else if(remRange < lowerRangeEdge) {
1225 
1226  G4double ratio = remRange / lowerRangeEdge;
1227  loss = kineticEnergy - ratio * ratio * lowerEnEdge;
1228  }
1229  else {
1230 
1231  G4double energy = rangeEnergy -> Value(range - stepLength);
1232  loss = kineticEnergy - energy;
1233  }
1234  }
1235 
1236  if(loss < 0.0) loss = 0.0;
1237 
1238  return loss;
1239 }
static constexpr double mm
Definition: G4SIunits.hh:115
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition const G4Material *G4double range
G4double energy(const ThreeVector &p, const G4double m)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:377

Here is the call graph for this function:

Here is the caller graph for this function:

void G4IonParametrisedLossModel::CorrectionsAlongStep ( const G4MaterialCutsCouple couple,
const G4DynamicParticle dynamicParticle,
G4double eloss,
G4double ,
G4double  length 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 918 of file G4IonParametrisedLossModel.cc.

923  {
924 
925  // ############## Corrections for along step energy loss calculation ######
926  // The computed energy loss (due to electronic stopping) is overwritten
927  // by this function if an ion data parameterization is available for the
928  // current ion-material pair.
929  // No action on the energy loss (due to electronic stopping) is performed
930  // if no parameterization is available. In this case the original
931  // generic ion tables (in combination with the effective charge) are used
932  // in the along step DoIt function.
933  //
934  //
935  // (Implementation partly adapted from G4BraggIonModel/G4BetheBlochModel)
936 
937  const G4ParticleDefinition* particle = dynamicParticle -> GetDefinition();
938  const G4Material* material = couple -> GetMaterial();
939 
940  G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
941 
942  if(kineticEnergy == eloss) { return; }
943 
944  G4double cutEnergy = DBL_MAX;
945  size_t cutIndex = couple -> GetIndex();
946  cutEnergy = cutEnergies[cutIndex];
947 
948  UpdateDEDXCache(particle, material, cutEnergy);
949 
950  LossTableList::iterator iter = dedxCacheIter;
951 
952  // If parameterization for ions is available the electronic energy loss
953  // is overwritten
954  if(iter != lossTableList.end()) {
955 
956  // The energy loss is calculated using the ComputeDEDXPerVolume function
957  // and the step length (it is assumed that dE/dx does not change
958  // considerably along the step)
959  eloss =
960  length * ComputeDEDXPerVolume(material, particle,
961  kineticEnergy, cutEnergy);
962 
963 #ifdef PRINT_DEBUG
964  G4cout.precision(6);
965  G4cout << "########################################################"
966  << G4endl
967  << "# G4IonParametrisedLossModel::CorrectionsAlongStep"
968  << G4endl
969  << "# cut(MeV) = " << cutEnergy/MeV
970  << G4endl;
971 
972  G4cout << "#"
973  << std::setw(13) << std::right << "E(MeV)"
974  << std::setw(14) << "l(um)"
975  << std::setw(14) << "l*dE/dx(MeV)"
976  << std::setw(14) << "(l*dE/dx)/E"
977  << G4endl
978  << "# ------------------------------------------------------"
979  << G4endl;
980 
981  G4cout << std::setw(14) << std::right << kineticEnergy / MeV
982  << std::setw(14) << length / um
983  << std::setw(14) << eloss / MeV
984  << std::setw(14) << eloss / kineticEnergy * 100.0
985  << G4endl;
986 #endif
987 
988  // If the energy loss exceeds a certain fraction of the kinetic energy
989  // (the fraction is indicated by the parameter "energyLossLimit") then
990  // the range tables are used to derive a more accurate value of the
991  // energy loss
992  if(eloss > energyLossLimit * kineticEnergy) {
993 
994  eloss = ComputeLossForStep(couple, particle,
995  kineticEnergy,length);
996 
997 #ifdef PRINT_DEBUG
998  G4cout << "# Correction applied:"
999  << G4endl;
1000 
1001  G4cout << std::setw(14) << std::right << kineticEnergy / MeV
1002  << std::setw(14) << length / um
1003  << std::setw(14) << eloss / MeV
1004  << std::setw(14) << eloss / kineticEnergy * 100.0
1005  << G4endl;
1006 #endif
1007 
1008  }
1009  }
1010 
1011  // For all corrections below a kinetic energy between the Pre- and
1012  // Post-step energy values is used
1013  G4double energy = kineticEnergy - eloss * 0.5;
1014  if(energy < 0.0) energy = kineticEnergy * 0.5;
1015 
1016  G4double chargeSquareRatio = corrections ->
1017  EffectiveChargeSquareRatio(particle,
1018  material,
1019  energy);
1020  GetModelOfFluctuations() -> SetParticleAndCharge(particle,
1021  chargeSquareRatio);
1022 
1023  // A correction is applied considering the change of the effective charge
1024  // along the step (the parameter "corrFactor" refers to the effective
1025  // charge at the beginning of the step). Note: the correction is not
1026  // applied for energy loss values deriving directly from parameterized
1027  // ion stopping power tables
1028  G4double transitionEnergy = dedxCacheTransitionEnergy;
1029 
1030  if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
1031  chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1032  material,
1033  energy);
1034 
1035  G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1036  eloss *= chargeSquareRatioCorr;
1037  }
1038  else if (iter == lossTableList.end()) {
1039 
1040  chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1041  material,
1042  energy);
1043 
1044  G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1045  eloss *= chargeSquareRatioCorr;
1046  }
1047 
1048  // Ion high order corrections are applied if the current model does not
1049  // overwrite the energy loss (i.e. when the effective charge approach is
1050  // used)
1051  if(iter == lossTableList.end()) {
1052 
1053  G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
1054  G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
1055 
1056  // Corrections are only applied in the Bethe-Bloch energy region
1057  if(scaledKineticEnergy > lowEnergyLimit)
1058  eloss += length *
1059  corrections -> IonHighOrderCorrections(particle, couple, energy);
1060  }
1061 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:612
string material
Definition: eplot.py:19
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
G4GLOB_DLL std::ostream G4cout
static constexpr double um
Definition: G4SIunits.hh:113
G4double energy(const ThreeVector &p, const G4double m)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
G4double ComputeLossForStep(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double, G4double)
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

G4double G4IonParametrisedLossModel::CrossSectionPerVolume ( const G4Material material,
const G4ParticleDefinition particle,
G4double  kineticEnergy,
G4double  cutEnergy,
G4double  maxEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 439 of file G4IonParametrisedLossModel.cc.

444  {
445 
446  G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
447  G4double cross = ComputeCrossSectionPerAtom(particle,
448  kineticEnergy,
449  nbElecPerVolume, 0,
450  cutEnergy,
451  maxEnergy);
452 
453  return cross;
454 }
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4IonParametrisedLossModel::DeactivateICRU73Scaling ( )

Definition at line 1325 of file G4IonParametrisedLossModel.cc.

1325  {
1326 
1327  RemoveDEDXTable("ICRU73");
1328  AddDEDXTable("ICRU73", new G4IonStoppingData("ion_stopping_data/icru73"));
1329 }
G4bool AddDEDXTable(const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
G4bool RemoveDEDXTable(const G4String &name)

Here is the call graph for this function:

G4double G4IonParametrisedLossModel::DeltaRayMeanEnergyTransferRate ( const G4Material ,
const G4ParticleDefinition ,
G4double  ,
G4double   
)
inline

Here is the caller graph for this function:

G4double G4IonParametrisedLossModel::GetChargeSquareRatio ( const G4ParticleDefinition particle,
const G4Material material,
G4double  kineticEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 228 of file G4IonParametrisedLossModel.cc.

231  { // Kinetic energy
232 
233  G4double chargeSquareRatio = corrections ->
234  EffectiveChargeSquareRatio(particle,
235  material,
236  kineticEnergy);
237  corrFactor = chargeSquareRatio *
238  corrections -> EffectiveChargeCorrection(particle,
239  material,
240  kineticEnergy);
241  return corrFactor;
242 }
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4IonParametrisedLossModel::GetParticleCharge ( const G4ParticleDefinition particle,
const G4Material material,
G4double  kineticEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 246 of file G4IonParametrisedLossModel.cc.

249  { // Kinetic energy
250 
251  return corrections -> GetParticleCharge(particle, material, kineticEnergy);
252 }
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double)
void G4IonParametrisedLossModel::Initialise ( const G4ParticleDefinition particle,
const G4DataVector cuts 
)
virtual

Implements G4VEmModel.

Definition at line 256 of file G4IonParametrisedLossModel.cc.

258  {
259 
260  // Cached parameters are reset
261  cacheParticle = 0;
262  cacheMass = 0;
263  cacheElecMassRatio = 0;
264  cacheChargeSquare = 0;
265 
266  // Cached parameters are reset
267  rangeCacheParticle = 0;
268  rangeCacheMatCutsCouple = 0;
269  rangeCacheEnergyRange = 0;
270  rangeCacheRangeEnergy = 0;
271 
272  // Cached parameters are reset
273  dedxCacheParticle = 0;
274  dedxCacheMaterial = 0;
275  dedxCacheEnergyCut = 0;
276  dedxCacheIter = lossTableList.end();
277  dedxCacheTransitionEnergy = 0.0;
278  dedxCacheTransitionFactor = 0.0;
279  dedxCacheGenIonMassRatio = 0.0;
280 
281  // The cache of loss tables is cleared
282  LossTableList::iterator iterTables = lossTableList.begin();
283  LossTableList::iterator iterTables_end = lossTableList.end();
284 
285  for(;iterTables != iterTables_end; iterTables++)
286  (*iterTables) -> ClearCache();
287 
288  // Range vs energy and energy vs range vectors from previous runs are
289  // cleared
290  RangeEnergyTable::iterator iterRange = r.begin();
291  RangeEnergyTable::iterator iterRange_end = r.end();
292 
293  for(;iterRange != iterRange_end; iterRange++) {
294  delete iterRange->second;
295  }
296  r.clear();
297 
298  EnergyRangeTable::iterator iterEnergy = E.begin();
299  EnergyRangeTable::iterator iterEnergy_end = E.end();
300 
301  for(;iterEnergy != iterEnergy_end; iterEnergy++) {
302  delete iterEnergy->second;
303  }
304  E.clear();
305 
306  // The cut energies are (re)loaded
307  size_t size = cuts.size();
308  cutEnergies.clear();
309  for(size_t i = 0; i < size; i++) cutEnergies.push_back(cuts[i]);
310 
311  // All dE/dx vectors are built
312  const G4ProductionCutsTable* coupleTable=
314  size_t nmbCouples = coupleTable -> GetTableSize();
315 
316 #ifdef PRINT_TABLE_BUILT
317  G4cout << "G4IonParametrisedLossModel::Initialise():"
318  << " Building dE/dx vectors:"
319  << G4endl;
320 #endif
321 
322  for (size_t i = 0; i < nmbCouples; i++) {
323 
324  const G4MaterialCutsCouple* couple =
325  coupleTable -> GetMaterialCutsCouple(i);
326 
327  const G4Material* material = couple -> GetMaterial();
328  // G4ProductionCuts* productionCuts = couple -> GetProductionCuts();
329 
330  for(G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) {
331 
332  LossTableList::iterator iter = lossTableList.begin();
333  LossTableList::iterator iter_end = lossTableList.end();
334 
335  for(;iter != iter_end; iter++) {
336 
337  if(*iter == 0) {
338  G4cout << "G4IonParametrisedLossModel::Initialise():"
339  << " Skipping illegal table."
340  << G4endl;
341  }
342 
343  G4bool isApplicable =
344  (*iter) -> BuildDEDXTable(atomicNumberIon, material);
345  if(isApplicable) {
346 
347 #ifdef PRINT_TABLE_BUILT
348  G4cout << " Atomic Number Ion = " << atomicNumberIon
349  << ", Material = " << material -> GetName()
350  << ", Table = " << (*iter) -> GetName()
351  << G4endl;
352 #endif
353  break;
354  }
355  }
356  }
357  }
358 
359  // The particle change object
360  if(! particleChangeLoss) {
361  particleChangeLoss = GetParticleChangeForLoss();
362  braggIonModel -> SetParticleChange(particleChangeLoss, 0);
363  betheBlochModel -> SetParticleChange(particleChangeLoss, 0);
364  }
365 
366  // The G4BraggIonModel and G4BetheBlochModel instances are initialised with
367  // the same settings as the current model:
368  braggIonModel -> Initialise(particle, cuts);
369  betheBlochModel -> Initialise(particle, cuts);
370 }
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
Definition: G4VEmModel.cc:418
int G4int
Definition: G4Types.hh:78
string material
Definition: eplot.py:19
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
static G4ProductionCutsTable * GetProductionCutsTable()
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
#define G4endl
Definition: G4ios.hh:61
const G4String & GetName() const
Definition: G4VEmModel.hh:802

Here is the call graph for this function:

LossTableList::iterator G4IonParametrisedLossModel::IsApplicable ( const G4ParticleDefinition ,
const G4Material  
)
inline

Here is the caller graph for this function:

G4double G4IonParametrisedLossModel::MaxSecondaryEnergy ( const G4ParticleDefinition particle,
G4double  kineticEnergy 
)
protectedvirtual

Reimplemented from G4VEmModel.

Definition at line 200 of file G4IonParametrisedLossModel.cc.

202  {
203 
204  // ############## Maximum energy of secondaries ##########################
205  // Function computes maximum energy of secondary electrons which are
206  // released by an ion
207  //
208  // See Geant4 physics reference manual (version 9.1), section 9.1.1
209  //
210  // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
211  // C.Caso et al. (Part. Data Group), Europ. Phys. Jour. C 3 1 (1998).
212  // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
213  //
214  // (Implementation adapted from G4BraggIonModel)
215 
216  if(particle != cacheParticle) UpdateCache(particle);
217 
218  G4double tau = kineticEnergy/cacheMass;
219  G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
220  (1. + 2.0 * (tau + 1.) * cacheElecMassRatio +
221  cacheElecMassRatio * cacheElecMassRatio);
222 
223  return tmax;
224 }
float electron_mass_c2
Definition: hepunit.py:274
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4IonParametrisedLossModel::MinEnergyCut ( const G4ParticleDefinition ,
const G4MaterialCutsCouple couple 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 190 of file G4IonParametrisedLossModel.cc.

192  {
193 
194  return couple -> GetMaterial() -> GetIonisation() ->
195  GetMeanExcitationEnergy();
196 }
void G4IonParametrisedLossModel::PrintDEDXTable ( const G4ParticleDefinition particle,
const G4Material material,
G4double  lowerBoundary,
G4double  upperBoundary,
G4int  numBins,
G4bool  logScaleEnergy 
)

Definition at line 607 of file G4IonParametrisedLossModel.cc.

613  { // Logarithmic scaling of energy
614 
615  G4double atomicMassNumber = particle -> GetAtomicMass();
616  G4double materialDensity = material -> GetDensity();
617 
618  G4cout << "# dE/dx table for " << particle -> GetParticleName()
619  << " in material " << material -> GetName()
620  << " of density " << materialDensity / g * cm3
621  << " g/cm3"
622  << G4endl
623  << "# Projectile mass number A1 = " << atomicMassNumber
624  << G4endl
625  << "# ------------------------------------------------------"
626  << G4endl;
627  G4cout << "#"
628  << std::setw(13) << std::right << "E"
629  << std::setw(14) << "E/A1"
630  << std::setw(14) << "dE/dx"
631  << std::setw(14) << "1/rho*dE/dx"
632  << G4endl;
633  G4cout << "#"
634  << std::setw(13) << std::right << "(MeV)"
635  << std::setw(14) << "(MeV)"
636  << std::setw(14) << "(MeV/cm)"
637  << std::setw(14) << "(MeV*cm2/mg)"
638  << G4endl
639  << "# ------------------------------------------------------"
640  << G4endl;
641 
642  G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
643  G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
644 
645  if(logScaleEnergy) {
646 
647  energyLowerBoundary = std::log(energyLowerBoundary);
648  energyUpperBoundary = std::log(energyUpperBoundary);
649  }
650 
651  G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
652  G4double(nmbBins);
653 
654  for(int i = 0; i < numBins + 1; i++) {
655 
656  G4double energy = energyLowerBoundary + i * deltaEnergy;
657  if(logScaleEnergy) energy = G4Exp(energy);
658 
659  G4double dedx = ComputeDEDXPerVolume(material, particle, energy, DBL_MAX);
660  G4cout.precision(6);
661  G4cout << std::setw(14) << std::right << energy / MeV
662  << std::setw(14) << energy / atomicMassNumber / MeV
663  << std::setw(14) << dedx / MeV * cm
664  << std::setw(14) << dedx / materialDensity / (MeV*cm2/(0.001*g))
665  << G4endl;
666  }
667 }
static constexpr double cm2
Definition: G4SIunits.hh:120
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
G4GLOB_DLL std::ostream G4cout
static constexpr double cm
Definition: G4SIunits.hh:119
static constexpr double cm3
Definition: G4SIunits.hh:121
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double energy(const ThreeVector &p, const G4double m)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
const G4String & GetName() const
Definition: G4VEmModel.hh:802
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

Here is the caller graph for this function:

void G4IonParametrisedLossModel::PrintDEDXTableHandlers ( const G4ParticleDefinition particle,
const G4Material material,
G4double  lowerBoundary,
G4double  upperBoundary,
G4int  numBins,
G4bool  logScaleEnergy 
)

Definition at line 671 of file G4IonParametrisedLossModel.cc.

677  { // Logarithmic scaling of energy
678 
679  LossTableList::iterator iter = lossTableList.begin();
680  LossTableList::iterator iter_end = lossTableList.end();
681 
682  for(;iter != iter_end; iter++) {
683  G4bool isApplicable = (*iter) -> IsApplicable(particle, material);
684  if(isApplicable) {
685  (*iter) -> PrintDEDXTable(particle, material,
686  lowerBoundary, upperBoundary,
687  numBins,logScaleEnergy);
688  break;
689  }
690  }
691 }
bool G4bool
Definition: G4Types.hh:79
LossTableList::iterator IsApplicable(const G4ParticleDefinition *, const G4Material *)
void PrintDEDXTable(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)

Here is the call graph for this function:

G4bool G4IonParametrisedLossModel::RemoveDEDXTable ( const G4String name)

Definition at line 1286 of file G4IonParametrisedLossModel.cc.

1287  {
1288 
1289  LossTableList::iterator iter = lossTableList.begin();
1290  LossTableList::iterator iter_end = lossTableList.end();
1291 
1292  for(;iter != iter_end; iter++) {
1293  G4String tableName = (*iter) -> GetName();
1294 
1295  if(tableName == nam) {
1296  delete (*iter);
1297 
1298  // Remove from table list
1299  lossTableList.erase(iter);
1300 
1301  // Range vs energy and energy vs range vectors are cleared
1302  RangeEnergyTable::iterator iterRange = r.begin();
1303  RangeEnergyTable::iterator iterRange_end = r.end();
1304 
1305  for(;iterRange != iterRange_end; iterRange++)
1306  delete iterRange -> second;
1307  r.clear();
1308 
1309  EnergyRangeTable::iterator iterEnergy = E.begin();
1310  EnergyRangeTable::iterator iterEnergy_end = E.end();
1311 
1312  for(;iterEnergy != iterEnergy_end; iterEnergy++)
1313  delete iterEnergy -> second;
1314  E.clear();
1315 
1316  return true;
1317  }
1318  }
1319 
1320  return false;
1321 }
static constexpr double second
Definition: G4SIunits.hh:157
const G4String & GetName() const
Definition: G4VEmModel.hh:802

Here is the call graph for this function:

Here is the caller graph for this function:

void G4IonParametrisedLossModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  secondaries,
const G4MaterialCutsCouple couple,
const G4DynamicParticle particle,
G4double  cutKinEnergySec,
G4double  userMaxKinEnergySec 
)
virtual

Implements G4VEmModel.

Definition at line 695 of file G4IonParametrisedLossModel.cc.

700  {
701 
702 
703  // ############## Sampling of secondaries #################################
704  // The probability density function (pdf) of the kinetic energy T of a
705  // secondary electron may be written as:
706  // pdf(T) = f(T) * g(T)
707  // where
708  // f(T) = (Tmax - Tcut) / (Tmax * Tcut) * (1 / T^2)
709  // g(T) = 1 - beta^2 * T / Tmax
710  // where Tmax is the maximum kinetic energy of the secondary, Tcut
711  // is the lower energy cut and beta is the kinetic energy of the
712  // projectile.
713  //
714  // Sampling of the kinetic energy of a secondary electron:
715  // 1) T0 is sampled from f(T) using the cumulated distribution function
716  // F(T) = int_Tcut^T f(T')dT'
717  // 2) T is accepted or rejected by evaluating the rejection function g(T)
718  // at the sampled energy T0 against a randomly sampled value
719  //
720  //
721  // See Geant4 physics reference manual (version 9.1), section 9.1.4
722  //
723  //
724  // Reference pdf: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
725  //
726  // (Implementation adapted from G4BraggIonModel)
727 
728  G4double rossiMaxKinEnergySec = MaxSecondaryKinEnergy(particle);
729  G4double maxKinEnergySec =
730  std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
731 
732  if(cutKinEnergySec >= maxKinEnergySec) return;
733 
734  G4double kineticEnergy = particle -> GetKineticEnergy();
735 
736  G4double energy = kineticEnergy + cacheMass;
737  G4double betaSquared = kineticEnergy * (energy + cacheMass)
738  / (energy * energy);
739 
740  G4double kinEnergySec;
741  G4double grej;
742 
743  do {
744 
745  // Sampling kinetic energy from f(T) (using F(T)):
746  G4double xi = G4UniformRand();
747  kinEnergySec = cutKinEnergySec * maxKinEnergySec /
748  (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
749 
750  // Deriving the value of the rejection function at the obtained kinetic
751  // energy:
752  grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
753 
754  if(grej > 1.0) {
755  G4cout << "G4IonParametrisedLossModel::SampleSecondary Warning: "
756  << "Majorant 1.0 < "
757  << grej << " for e= " << kinEnergySec
758  << G4endl;
759  }
760 
761  } while( G4UniformRand() >= grej );
762 
763  const G4Material* mat = couple->GetMaterial();
765 
767 
768  G4DynamicParticle* delta = new G4DynamicParticle(electron,
769  GetAngularDistribution()->SampleDirection(particle, kinEnergySec,
770  Z, mat),
771  kinEnergySec);
772 
773 
774  secondaries->push_back(delta);
775 
776  // Change kinematics of primary particle
777  G4ThreeVector direction = particle ->GetMomentumDirection();
778  G4double totalMomentum = std::sqrt(kineticEnergy*(energy + cacheMass));
779 
780  G4ThreeVector finalP = totalMomentum*direction - delta->GetMomentum();
781  finalP = finalP.unit();
782 
783  kineticEnergy -= kinEnergySec;
784 
785  particleChangeLoss->SetProposedKineticEnergy(kineticEnergy);
786  particleChangeLoss->SetProposedMomentumDirection(finalP);
787 }
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:484
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:619
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
G4double energy(const ThreeVector &p, const G4double m)
Hep3Vector unit() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4ThreeVector GetMomentum() const
const G4Material * GetMaterial() const
G4int SelectRandomAtomNumber(const G4Material *)
Definition: G4VEmModel.hh:564

Here is the call graph for this function:

void G4IonParametrisedLossModel::SetEnergyLossLimit ( G4double  ionEnergyLossLimit)
inline

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