Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4IonParametrisedLossModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4IonParametrisedLossModel.cc 95950 2016-03-03 10:42:48Z gcosmo $
27 //
28 // ===========================================================================
29 // GEANT4 class source file
30 //
31 // Class: G4IonParametrisedLossModel
32 //
33 // Base class: G4VEmModel (utils)
34 //
35 // Author: Anton Lechner (Anton.Lechner@cern.ch)
36 //
37 // First implementation: 10. 11. 2008
38 //
39 // Modifications: 03. 02. 2009 - Bug fix iterators (AL)
40 // 11. 03. 2009 - Introduced new table handler(G4IonDEDXHandler)
41 // and modified method to add/remove tables
42 // (tables are now built in init. phase),
43 // Minor bug fix in ComputeDEDXPerVolume (AL)
44 // 11. 05. 2009 - Introduced scaling algorithm for heavier ions:
45 // G4IonDEDXScalingICRU73 (AL)
46 // 12. 11. 2009 - Moved from original ICRU 73 classes to new
47 // class (G4IonStoppingData), which is capable
48 // of reading stopping power data files stored
49 // in G4LEDATA (requires G4EMLOW6.8 or higher).
50 // Simultanesouly, the upper energy limit of
51 // ICRU 73 is increased to 1 GeV/nucleon.
52 // - Removed nuclear stopping from Corrections-
53 // AlongStep since dedicated process was created.
54 // - Added function for switching off scaling
55 // of heavy ions from ICRU 73 data
56 // - Minor fix in ComputeLossForStep function
57 // - Minor fix in ComputeDEDXPerVolume (AL)
58 // 23. 11. 2009 - Changed energy loss limit from 0.15 to 0.01
59 // to improve accuracy for large steps (AL)
60 // 24. 11. 2009 - Bug fix: Range calculation corrected if same
61 // materials appears with different cuts in diff.
62 // regions (added UpdateRangeCache function and
63 // modified BuildRangeVector, ComputeLossForStep
64 // functions accordingly, added new cache param.)
65 // - Removed GetRange function (AL)
66 // 04. 11. 2010 - Moved virtual methods to the source (VI)
67 //
68 //
69 // Class description:
70 // Model for computing the energy loss of ions by employing a
71 // parameterisation of dE/dx tables (by default ICRU 73 tables). For
72 // ion-material combinations and/or projectile energies not covered
73 // by this model, the G4BraggIonModel and G4BetheBloch models are
74 // employed.
75 //
76 // Comments:
77 //
78 // ===========================================================================
79 
80 
82 #include "G4PhysicalConstants.hh"
83 #include "G4SystemOfUnits.hh"
84 #include "G4LPhysicsFreeVector.hh"
85 #include "G4IonStoppingData.hh"
86 #include "G4VIonDEDXTable.hh"
89 #include "G4BraggIonModel.hh"
90 #include "G4BetheBlochModel.hh"
91 #include "G4ProductionCutsTable.hh"
93 #include "G4LossTableManager.hh"
94 #include "G4GenericIon.hh"
95 #include "G4Electron.hh"
96 #include "G4DeltaAngle.hh"
97 #include "Randomize.hh"
98 #include "G4Exp.hh"
99 
100 //#define PRINT_TABLE_BUILT
101 
102 
103 // #########################################################################
104 
106  const G4ParticleDefinition*,
107  const G4String& nam)
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 }
162 
163 // #########################################################################
164 
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 }
187 
188 // #########################################################################
189 
191  const G4ParticleDefinition*,
192  const G4MaterialCutsCouple* couple) {
193 
194  return couple -> GetMaterial() -> GetIonisation() ->
195  GetMeanExcitationEnergy();
196 }
197 
198 // #########################################################################
199 
201  const G4ParticleDefinition* particle,
202  G4double kineticEnergy) {
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 }
225 
226 // #########################################################################
227 
229  const G4ParticleDefinition* particle,
230  const G4Material* material,
231  G4double kineticEnergy) { // 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 }
243 
244 // #########################################################################
245 
247  const G4ParticleDefinition* particle,
248  const G4Material* material,
249  G4double kineticEnergy) { // Kinetic energy
250 
251  return corrections -> GetParticleCharge(particle, material, kineticEnergy);
252 }
253 
254 // #########################################################################
255 
257  const G4ParticleDefinition* particle,
258  const G4DataVector& cuts) {
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 }
371 
372 // #########################################################################
373 
375  const G4ParticleDefinition* particle,
376  G4double kineticEnergy,
377  G4double atomicNumber,
378  G4double,
379  G4double cutEnergy,
380  G4double maxKinEnergy) {
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 }
436 
437 // #########################################################################
438 
440  const G4Material* material,
441  const G4ParticleDefinition* particle,
442  G4double kineticEnergy,
443  G4double cutEnergy,
444  G4double maxEnergy) {
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 }
455 
456 // #########################################################################
457 
459  const G4Material* material,
460  const G4ParticleDefinition* particle,
461  G4double kineticEnergy,
462  G4double cutEnergy) {
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 }
604 
605 // #########################################################################
606 
608  const G4ParticleDefinition* particle, // Projectile (ion)
609  const G4Material* material, // Absorber material
610  G4double lowerBoundary, // Minimum energy per nucleon
611  G4double upperBoundary, // Maximum energy per nucleon
612  G4int numBins, // Number of bins
613  G4bool logScaleEnergy) { // 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 }
668 
669 // #########################################################################
670 
672  const G4ParticleDefinition* particle, // Projectile (ion)
673  const G4Material* material, // Absorber material
674  G4double lowerBoundary, // Minimum energy per nucleon
675  G4double upperBoundary, // Maximum energy per nucleon
676  G4int numBins, // Number of bins
677  G4bool logScaleEnergy) { // 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 }
692 
693 // #########################################################################
694 
696  std::vector<G4DynamicParticle*>* secondaries,
697  const G4MaterialCutsCouple* couple,
698  const G4DynamicParticle* particle,
699  G4double cutKinEnergySec,
700  G4double userMaxKinEnergySec) {
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 }
788 
789 // #########################################################################
790 
791 void G4IonParametrisedLossModel::UpdateRangeCache(
792  const G4ParticleDefinition* particle,
793  const G4MaterialCutsCouple* matCutsCouple) {
794 
795  // ############## Caching ##################################################
796  // If the ion-material-cut combination is covered by any native ion data
797  // parameterisation (for low energies), range vectors are computed
798 
799  if(particle == rangeCacheParticle &&
800  matCutsCouple == rangeCacheMatCutsCouple) {
801  }
802  else{
803  rangeCacheParticle = particle;
804  rangeCacheMatCutsCouple = matCutsCouple;
805 
806  const G4Material* material = matCutsCouple -> GetMaterial();
807  LossTableList::iterator iter = IsApplicable(particle, material);
808 
809  // If any table is applicable, the transition factor is computed:
810  if(iter != lossTableList.end()) {
811 
812  // Build range-energy and energy-range vectors if they don't exist
813  IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
814  RangeEnergyTable::iterator iterRange = r.find(ionMatCouple);
815 
816  if(iterRange == r.end()) BuildRangeVector(particle, matCutsCouple);
817 
818  rangeCacheEnergyRange = E[ionMatCouple];
819  rangeCacheRangeEnergy = r[ionMatCouple];
820  }
821  else {
822  rangeCacheEnergyRange = 0;
823  rangeCacheRangeEnergy = 0;
824  }
825  }
826 }
827 
828 // #########################################################################
829 
830 void G4IonParametrisedLossModel::UpdateDEDXCache(
831  const G4ParticleDefinition* particle,
832  const G4Material* material,
833  G4double cutEnergy) {
834 
835  // ############## Caching ##################################################
836  // If the ion-material combination is covered by any native ion data
837  // parameterisation (for low energies), a transition factor is computed
838  // which is applied to Bethe-Bloch results at higher energies to
839  // guarantee a smooth transition.
840  // This factor only needs to be calculated for the first step an ion
841  // performs inside a certain material.
842 
843  if(particle == dedxCacheParticle &&
844  material == dedxCacheMaterial &&
845  cutEnergy == dedxCacheEnergyCut) {
846  }
847  else {
848 
849  dedxCacheParticle = particle;
850  dedxCacheMaterial = material;
851  dedxCacheEnergyCut = cutEnergy;
852 
853  G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
854  dedxCacheGenIonMassRatio = massRatio;
855 
856  LossTableList::iterator iter = IsApplicable(particle, material);
857  dedxCacheIter = iter;
858 
859  // If any table is applicable, the transition factor is computed:
860  if(iter != lossTableList.end()) {
861 
862  // Retrieving the transition energy from the parameterisation table
863  G4double transitionEnergy =
864  (*iter) -> GetUpperEnergyEdge(particle, material);
865  dedxCacheTransitionEnergy = transitionEnergy;
866 
867  // Computing dE/dx from low-energy parameterisation at
868  // transition energy
869  G4double dEdxParam = (*iter) -> GetDEDX(particle, material,
870  transitionEnergy);
871 
872  G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material,
873  particle,
874  transitionEnergy,
875  cutEnergy);
876  dEdxParam -= dEdxDeltaRays;
877 
878  // Computing dE/dx from Bethe-Bloch formula at transition
879  // energy
880  G4double transitionChargeSquare =
881  GetChargeSquareRatio(particle, material, transitionEnergy);
882 
883  G4double scaledTransitionEnergy = transitionEnergy * massRatio;
884 
885  G4double dEdxBetheBloch =
886  betheBlochModel -> ComputeDEDXPerVolume(
887  material, genericIon,
888  scaledTransitionEnergy, cutEnergy);
889  dEdxBetheBloch *= transitionChargeSquare;
890 
891  // Additionally, high order corrections are added
892  dEdxBetheBloch +=
893  corrections -> ComputeIonCorrections(particle,
894  material, transitionEnergy);
895 
896  // Computing transition factor from both dE/dx values
897  dedxCacheTransitionFactor =
898  (dEdxParam - dEdxBetheBloch)/dEdxBetheBloch
899  * transitionEnergy;
900  }
901  else {
902 
903  dedxCacheParticle = particle;
904  dedxCacheMaterial = material;
905  dedxCacheEnergyCut = cutEnergy;
906 
907  dedxCacheGenIonMassRatio =
908  genericIonPDGMass / particle -> GetPDGMass();
909 
910  dedxCacheTransitionEnergy = 0.0;
911  dedxCacheTransitionFactor = 0.0;
912  }
913  }
914 }
915 
916 // #########################################################################
917 
919  const G4MaterialCutsCouple* couple,
920  const G4DynamicParticle* dynamicParticle,
921  G4double& eloss,
922  G4double&,
923  G4double length) {
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 }
1062 
1063 // #########################################################################
1064 
1065 void G4IonParametrisedLossModel::BuildRangeVector(
1066  const G4ParticleDefinition* particle,
1067  const G4MaterialCutsCouple* matCutsCouple) {
1068 
1069  G4double cutEnergy = DBL_MAX;
1070  size_t cutIndex = matCutsCouple -> GetIndex();
1071  cutEnergy = cutEnergies[cutIndex];
1072 
1073  const G4Material* material = matCutsCouple -> GetMaterial();
1074 
1075  G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
1076 
1077  G4double lowerEnergy = lowerEnergyEdgeIntegr / massRatio;
1078  G4double upperEnergy = upperEnergyEdgeIntegr / massRatio;
1079 
1080  G4double logLowerEnergyEdge = std::log(lowerEnergy);
1081  G4double logUpperEnergyEdge = std::log(upperEnergy);
1082 
1083  G4double logDeltaEnergy = (logUpperEnergyEdge - logLowerEnergyEdge) /
1084  G4double(nmbBins);
1085 
1086  G4double logDeltaIntegr = logDeltaEnergy / G4double(nmbSubBins);
1087 
1088  G4LPhysicsFreeVector* energyRangeVector =
1089  new G4LPhysicsFreeVector(nmbBins+1,
1090  lowerEnergy,
1091  upperEnergy);
1092 
1093  G4double dedxLow = ComputeDEDXPerVolume(material,
1094  particle,
1095  lowerEnergy,
1096  cutEnergy);
1097 
1098  G4double range = 2.0 * lowerEnergy / dedxLow;
1099 
1100  energyRangeVector -> PutValues(0, lowerEnergy, range);
1101 
1102  G4double logEnergy = std::log(lowerEnergy);
1103  for(size_t i = 1; i < nmbBins+1; i++) {
1104 
1105  G4double logEnergyIntegr = logEnergy;
1106 
1107  for(size_t j = 0; j < nmbSubBins; j++) {
1108 
1109  G4double binLowerBoundary = G4Exp(logEnergyIntegr);
1110  logEnergyIntegr += logDeltaIntegr;
1111 
1112  G4double binUpperBoundary = G4Exp(logEnergyIntegr);
1113  G4double deltaIntegr = binUpperBoundary - binLowerBoundary;
1114 
1115  G4double energyIntegr = binLowerBoundary + 0.5 * deltaIntegr;
1116 
1117  G4double dedxValue = ComputeDEDXPerVolume(material,
1118  particle,
1119  energyIntegr,
1120  cutEnergy);
1121 
1122  if(dedxValue > 0.0) range += deltaIntegr / dedxValue;
1123 
1124 #ifdef PRINT_DEBUG_DETAILS
1125  G4cout << " E = "<< energyIntegr/MeV
1126  << " MeV -> dE = " << deltaIntegr/MeV
1127  << " MeV -> dE/dx = " << dedxValue/MeV*mm
1128  << " MeV/mm -> dE/(dE/dx) = " << deltaIntegr /
1129  dedxValue / mm
1130  << " mm -> range = " << range / mm
1131  << " mm " << G4endl;
1132 #endif
1133  }
1134 
1135  logEnergy += logDeltaEnergy;
1136 
1137  G4double energy = G4Exp(logEnergy);
1138 
1139  energyRangeVector -> PutValues(i, energy, range);
1140 
1141 #ifdef PRINT_DEBUG_DETAILS
1142  G4cout << "G4IonParametrisedLossModel::BuildRangeVector() bin = "
1143  << i <<", E = "
1144  << energy / MeV << " MeV, R = "
1145  << range / mm << " mm"
1146  << G4endl;
1147 #endif
1148 
1149  }
1150  energyRangeVector -> SetSpline(true);
1151 
1152  G4double lowerRangeEdge =
1153  energyRangeVector -> Value(lowerEnergy);
1154  G4double upperRangeEdge =
1155  energyRangeVector -> Value(upperEnergy);
1156 
1157  G4LPhysicsFreeVector* rangeEnergyVector
1158  = new G4LPhysicsFreeVector(nmbBins+1,
1159  lowerRangeEdge,
1160  upperRangeEdge);
1161 
1162  for(size_t i = 0; i < nmbBins+1; i++) {
1163  G4double energy = energyRangeVector -> Energy(i);
1164  rangeEnergyVector ->
1165  PutValues(i, energyRangeVector -> Value(energy), energy);
1166  }
1167 
1168  rangeEnergyVector -> SetSpline(true);
1169 
1170 #ifdef PRINT_DEBUG_TABLES
1171  G4cout << *energyLossVector
1172  << *energyRangeVector
1173  << *rangeEnergyVector << G4endl;
1174 #endif
1175 
1176  IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
1177 
1178  E[ionMatCouple] = energyRangeVector;
1179  r[ionMatCouple] = rangeEnergyVector;
1180 }
1181 
1182 // #########################################################################
1183 
1185  const G4MaterialCutsCouple* matCutsCouple,
1186  const G4ParticleDefinition* particle,
1187  G4double kineticEnergy,
1188  G4double stepLength) {
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 }
1240 
1241 // #########################################################################
1242 
1244  const G4String& nam,
1245  G4VIonDEDXTable* table,
1246  G4VIonDEDXScalingAlgorithm* algorithm) {
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 }
1283 
1284 // #########################################################################
1285 
1287  const G4String& nam) {
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 }
1322 
1323 // #########################################################################
1324 
1326 
1327  RemoveDEDXTable("ICRU73");
1328  AddDEDXTable("ICRU73", new G4IonStoppingData("ion_stopping_data/icru73"));
1329 }
1330 
1331 // #########################################################################
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:484
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
static constexpr double mm
Definition: G4SIunits.hh:115
static constexpr double cm2
Definition: G4SIunits.hh:120
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
Definition: G4VEmModel.cc:418
static G4GenericIon * Definition()
Definition: G4GenericIon.cc:49
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:619
G4bool AddDEDXTable(const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
static constexpr double second
Definition: G4SIunits.hh:157
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:612
int G4int
Definition: G4Types.hh:78
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
static constexpr double TeV
Definition: G4SIunits.hh:218
string material
Definition: eplot.py:19
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
static constexpr double um
Definition: G4SIunits.hh:113
const G4ParticleDefinition const G4Material *G4double range
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
static constexpr double cm
Definition: G4SIunits.hh:119
void SetProposedKineticEnergy(G4double proposedKinEnergy)
float electron_mass_c2
Definition: hepunit.py:274
void SetProposedMomentumDirection(const G4ThreeVector &dir)
static constexpr double cm3
Definition: G4SIunits.hh:121
G4double DeltaRayMeanEnergyTransferRate(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
void PrintDEDXTableHandlers(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
static G4ProductionCutsTable * GetProductionCutsTable()
G4bool RemoveDEDXTable(const G4String &name)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double)
G4double energy(const ThreeVector &p, const G4double m)
std::pair< const G4ParticleDefinition *, const G4MaterialCutsCouple * > IonMatCouple
Hep3Vector unit() const
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:626
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double, G4double)
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
const G4String & GetName() const
Definition: G4VEmModel.hh:802
G4IonParametrisedLossModel(const G4ParticleDefinition *particle=0, const G4String &name="ParamICRU73")
LossTableList::iterator IsApplicable(const G4ParticleDefinition *, const G4Material *)
double G4double
Definition: G4Types.hh:76
void PrintDEDXTable(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
G4double ComputeLossForStep(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double, G4double)
#define DBL_MAX
Definition: templates.hh:83
G4ThreeVector GetMomentum() const
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:377
const G4Material * GetMaterial() const
G4GLOB_DLL std::ostream G4cerr
G4int SelectRandomAtomNumber(const G4Material *)
Definition: G4VEmModel.hh:564
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &, G4double &, G4double)