Geant4_10
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 76243 2013-11-08 11:11:38Z 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 "Randomize.hh"
97 
98 //#define PRINT_TABLE_BUILT
99 
100 
101 // #########################################################################
102 
104  const G4ParticleDefinition*,
105  const G4String& nam)
106  : G4VEmModel(nam),
107  braggIonModel(0),
108  betheBlochModel(0),
109  nmbBins(90),
110  nmbSubBins(100),
111  particleChangeLoss(0),
112  corrFactor(1.0),
113  energyLossLimit(0.01),
114  cutEnergies(0)
115 {
116  genericIon = G4GenericIon::Definition();
117  genericIonPDGMass = genericIon -> GetPDGMass();
118  corrections = G4LossTableManager::Instance() -> EmCorrections();
119 
120  // The upper limit of the current model is set to 100 TeV
121  SetHighEnergyLimit(100.0 * TeV);
122 
123  // The Bragg ion and Bethe Bloch models are instantiated
124  braggIonModel = new G4BraggIonModel();
125  betheBlochModel = new G4BetheBlochModel();
126 
127  // By default ICRU 73 stopping power tables are loaded:
128  AddDEDXTable("ICRU73",
129  new G4IonStoppingData("ion_stopping_data/icru73"),
130  new G4IonDEDXScalingICRU73());
131 
132  // The boundaries for the range tables are set
133  lowerEnergyEdgeIntegr = 0.025 * MeV;
134  upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit();
135 
136  // Cache parameters are set
137  cacheParticle = 0;
138  cacheMass = 0;
139  cacheElecMassRatio = 0;
140  cacheChargeSquare = 0;
141 
142  // Cache parameters are set
143  rangeCacheParticle = 0;
144  rangeCacheMatCutsCouple = 0;
145  rangeCacheEnergyRange = 0;
146  rangeCacheRangeEnergy = 0;
147 
148  // Cache parameters are set
149  dedxCacheParticle = 0;
150  dedxCacheMaterial = 0;
151  dedxCacheEnergyCut = 0;
152  dedxCacheIter = lossTableList.end();
153  dedxCacheTransitionEnergy = 0.0;
154  dedxCacheTransitionFactor = 0.0;
155  dedxCacheGenIonMassRatio = 0.0;
156 }
157 
158 // #########################################################################
159 
161 
162  // dE/dx table objects are deleted and the container is cleared
163  LossTableList::iterator iterTables = lossTableList.begin();
164  LossTableList::iterator iterTables_end = lossTableList.end();
165 
166  for(;iterTables != iterTables_end; iterTables++) delete *iterTables;
167  lossTableList.clear();
168 }
169 
170 // #########################################################################
171 
173  const G4ParticleDefinition*,
174  const G4MaterialCutsCouple* couple) {
175 
176  return couple -> GetMaterial() -> GetIonisation() ->
177  GetMeanExcitationEnergy();
178 }
179 
180 // #########################################################################
181 
183  const G4ParticleDefinition* particle,
184  G4double kineticEnergy) {
185 
186  // ############## Maximum energy of secondaries ##########################
187  // Function computes maximum energy of secondary electrons which are
188  // released by an ion
189  //
190  // See Geant4 physics reference manual (version 9.1), section 9.1.1
191  //
192  // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
193  // C.Caso et al. (Part. Data Group), Europ. Phys. Jour. C 3 1 (1998).
194  // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
195  //
196  // (Implementation adapted from G4BraggIonModel)
197 
198  if(particle != cacheParticle) UpdateCache(particle);
199 
200  G4double tau = kineticEnergy/cacheMass;
201  G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
202  (1. + 2.0 * (tau + 1.) * cacheElecMassRatio +
203  cacheElecMassRatio * cacheElecMassRatio);
204 
205  return tmax;
206 }
207 
208 // #########################################################################
209 
211  const G4ParticleDefinition* particle,
212  const G4Material* material,
213  G4double kineticEnergy) { // Kinetic energy
214 
215  G4double chargeSquareRatio = corrections ->
216  EffectiveChargeSquareRatio(particle,
217  material,
218  kineticEnergy);
219  corrFactor = chargeSquareRatio *
220  corrections -> EffectiveChargeCorrection(particle,
221  material,
222  kineticEnergy);
223  return corrFactor;
224 }
225 
226 // #########################################################################
227 
229  const G4ParticleDefinition* particle,
230  const G4Material* material,
231  G4double kineticEnergy) { // Kinetic energy
232 
233  return corrections -> GetParticleCharge(particle, material, kineticEnergy);
234 }
235 
236 // #########################################################################
237 
239  const G4ParticleDefinition* particle,
240  const G4DataVector& cuts) {
241 
242  // Cached parameters are reset
243  cacheParticle = 0;
244  cacheMass = 0;
245  cacheElecMassRatio = 0;
246  cacheChargeSquare = 0;
247 
248  // Cached parameters are reset
249  rangeCacheParticle = 0;
250  rangeCacheMatCutsCouple = 0;
251  rangeCacheEnergyRange = 0;
252  rangeCacheRangeEnergy = 0;
253 
254  // Cached parameters are reset
255  dedxCacheParticle = 0;
256  dedxCacheMaterial = 0;
257  dedxCacheEnergyCut = 0;
258  dedxCacheIter = lossTableList.end();
259  dedxCacheTransitionEnergy = 0.0;
260  dedxCacheTransitionFactor = 0.0;
261  dedxCacheGenIonMassRatio = 0.0;
262 
263  // The cache of loss tables is cleared
264  LossTableList::iterator iterTables = lossTableList.begin();
265  LossTableList::iterator iterTables_end = lossTableList.end();
266 
267  for(;iterTables != iterTables_end; iterTables++)
268  (*iterTables) -> ClearCache();
269 
270  // Range vs energy and energy vs range vectors from previous runs are
271  // cleared
272  RangeEnergyTable::iterator iterRange = r.begin();
273  RangeEnergyTable::iterator iterRange_end = r.end();
274 
275  for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second;
276  r.clear();
277 
278  EnergyRangeTable::iterator iterEnergy = E.begin();
279  EnergyRangeTable::iterator iterEnergy_end = E.end();
280 
281  for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second;
282  E.clear();
283 
284  // The cut energies are (re)loaded
285  size_t size = cuts.size();
286  cutEnergies.clear();
287  for(size_t i = 0; i < size; i++) cutEnergies.push_back(cuts[i]);
288 
289  // All dE/dx vectors are built
290  const G4ProductionCutsTable* coupleTable=
292  size_t nmbCouples = coupleTable -> GetTableSize();
293 
294 #ifdef PRINT_TABLE_BUILT
295  G4cout << "G4IonParametrisedLossModel::Initialise():"
296  << " Building dE/dx vectors:"
297  << G4endl;
298 #endif
299 
300  for (size_t i = 0; i < nmbCouples; i++) {
301 
302  const G4MaterialCutsCouple* couple =
303  coupleTable -> GetMaterialCutsCouple(i);
304 
305  const G4Material* material = couple -> GetMaterial();
306  // G4ProductionCuts* productionCuts = couple -> GetProductionCuts();
307 
308  for(G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) {
309 
310  LossTableList::iterator iter = lossTableList.begin();
311  LossTableList::iterator iter_end = lossTableList.end();
312 
313  for(;iter != iter_end; iter++) {
314 
315  if(*iter == 0) {
316  G4cout << "G4IonParametrisedLossModel::Initialise():"
317  << " Skipping illegal table."
318  << G4endl;
319  }
320 
321  G4bool isApplicable =
322  (*iter) -> BuildDEDXTable(atomicNumberIon, material);
323  if(isApplicable) {
324 
325 #ifdef PRINT_TABLE_BUILT
326  G4cout << " Atomic Number Ion = " << atomicNumberIon
327  << ", Material = " << material -> GetName()
328  << ", Table = " << (*iter) -> GetName()
329  << G4endl;
330 #endif
331  break;
332  }
333  }
334  }
335  }
336 
337  // The particle change object
338  if(! particleChangeLoss) {
339  particleChangeLoss = GetParticleChangeForLoss();
340  braggIonModel -> SetParticleChange(particleChangeLoss, 0);
341  betheBlochModel -> SetParticleChange(particleChangeLoss, 0);
342  }
343 
344  // The G4BraggIonModel and G4BetheBlochModel instances are initialised with
345  // the same settings as the current model:
346  braggIonModel -> Initialise(particle, cuts);
347  betheBlochModel -> Initialise(particle, cuts);
348 }
349 
350 // #########################################################################
351 
353  const G4ParticleDefinition* particle,
354  G4double kineticEnergy,
355  G4double atomicNumber,
356  G4double,
357  G4double cutEnergy,
358  G4double maxKinEnergy) {
359 
360  // ############## Cross section per atom ################################
361  // Function computes ionization cross section per atom
362  //
363  // See Geant4 physics reference manual (version 9.1), section 9.1.3
364  //
365  // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
366  // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
367  //
368  // (Implementation adapted from G4BraggIonModel)
369 
370  G4double crosssection = 0.0;
371  G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
372  G4double maxEnergy = std::min(tmax, maxKinEnergy);
373 
374  if(cutEnergy < tmax) {
375 
376  G4double energy = kineticEnergy + cacheMass;
377  G4double betaSquared = kineticEnergy *
378  (energy + cacheMass) / (energy * energy);
379 
380  crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
381  betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
382 
383  crosssection *= twopi_mc2_rcl2 * cacheChargeSquare / betaSquared;
384  }
385 
386 #ifdef PRINT_DEBUG_CS
387  G4cout << "########################################################"
388  << G4endl
389  << "# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
390  << G4endl
391  << "# particle =" << particle -> GetParticleName()
392  << G4endl
393  << "# cut(MeV) = " << cutEnergy/MeV
394  << G4endl;
395 
396  G4cout << "#"
397  << std::setw(13) << std::right << "E(MeV)"
398  << std::setw(14) << "CS(um)"
399  << std::setw(14) << "E_max_sec(MeV)"
400  << G4endl
401  << "# ------------------------------------------------------"
402  << G4endl;
403 
404  G4cout << std::setw(14) << std::right << kineticEnergy / MeV
405  << std::setw(14) << crosssection / (um * um)
406  << std::setw(14) << tmax / MeV
407  << G4endl;
408 #endif
409 
410  crosssection *= atomicNumber;
411 
412  return crosssection;
413 }
414 
415 // #########################################################################
416 
418  const G4Material* material,
419  const G4ParticleDefinition* particle,
420  G4double kineticEnergy,
421  G4double cutEnergy,
422  G4double maxEnergy) {
423 
424  G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
425  G4double cross = ComputeCrossSectionPerAtom(particle,
426  kineticEnergy,
427  nbElecPerVolume, 0,
428  cutEnergy,
429  maxEnergy);
430 
431  return cross;
432 }
433 
434 // #########################################################################
435 
437  const G4Material* material,
438  const G4ParticleDefinition* particle,
439  G4double kineticEnergy,
440  G4double cutEnergy) {
441 
442  // ############## dE/dx ##################################################
443  // Function computes dE/dx values, where following rules are adopted:
444  // A. If the ion-material pair is covered by any native ion data
445  // parameterisation, then:
446  // * This parameterization is used for energies below a given energy
447  // limit,
448  // * whereas above the limit the Bethe-Bloch model is applied, in
449  // combination with an effective charge estimate and high order
450  // correction terms.
451  // A smoothing procedure is applied to dE/dx values computed with
452  // the second approach. The smoothing factor is based on the dE/dx
453  // values of both approaches at the transition energy (high order
454  // correction terms are included in the calculation of the transition
455  // factor).
456  // B. If the particle is a generic ion, the BraggIon and Bethe-Bloch
457  // models are used and a smoothing procedure is applied to values
458  // obtained with the second approach.
459  // C. If the ion-material is not covered by any ion data parameterization
460  // then:
461  // * The BraggIon model is used for energies below a given energy
462  // limit,
463  // * whereas above the limit the Bethe-Bloch model is applied, in
464  // combination with an effective charge estimate and high order
465  // correction terms.
466  // Also in this case, a smoothing procedure is applied to dE/dx values
467  // computed with the second model.
468 
469  G4double dEdx = 0.0;
470  UpdateDEDXCache(particle, material, cutEnergy);
471 
472  LossTableList::iterator iter = dedxCacheIter;
473 
474  if(iter != lossTableList.end()) {
475 
476  G4double transitionEnergy = dedxCacheTransitionEnergy;
477 
478  if(transitionEnergy > kineticEnergy) {
479 
480  dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
481 
482  G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material,
483  particle,
484  kineticEnergy,
485  cutEnergy);
486  dEdx -= dEdxDeltaRays;
487  }
488  else {
489  G4double massRatio = dedxCacheGenIonMassRatio;
490 
491  G4double chargeSquare =
492  GetChargeSquareRatio(particle, material, kineticEnergy);
493 
494  G4double scaledKineticEnergy = kineticEnergy * massRatio;
495  G4double scaledTransitionEnergy = transitionEnergy * massRatio;
496 
497  G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
498 
499  if(scaledTransitionEnergy >= lowEnergyLimit) {
500 
501  dEdx = betheBlochModel -> ComputeDEDXPerVolume(
502  material, genericIon,
503  scaledKineticEnergy, cutEnergy);
504 
505  dEdx *= chargeSquare;
506 
507  dEdx += corrections -> ComputeIonCorrections(particle,
508  material, kineticEnergy);
509 
510  G4double factor = 1.0 + dedxCacheTransitionFactor /
511  kineticEnergy;
512 
513  dEdx *= factor;
514  }
515  }
516  }
517  else {
518  G4double massRatio = 1.0;
519  G4double chargeSquare = 1.0;
520 
521  if(particle != genericIon) {
522 
523  chargeSquare = GetChargeSquareRatio(particle, material, kineticEnergy);
524  massRatio = genericIonPDGMass / particle -> GetPDGMass();
525  }
526 
527  G4double scaledKineticEnergy = kineticEnergy * massRatio;
528 
529  G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
530  if(scaledKineticEnergy < lowEnergyLimit) {
531  dEdx = braggIonModel -> ComputeDEDXPerVolume(
532  material, genericIon,
533  scaledKineticEnergy, cutEnergy);
534 
535  dEdx *= chargeSquare;
536  }
537  else {
538  G4double dEdxLimitParam = braggIonModel -> ComputeDEDXPerVolume(
539  material, genericIon,
540  lowEnergyLimit, cutEnergy);
541 
542  G4double dEdxLimitBetheBloch = betheBlochModel -> ComputeDEDXPerVolume(
543  material, genericIon,
544  lowEnergyLimit, cutEnergy);
545 
546  if(particle != genericIon) {
547  G4double chargeSquareLowEnergyLimit =
548  GetChargeSquareRatio(particle, material,
549  lowEnergyLimit / massRatio);
550 
551  dEdxLimitParam *= chargeSquareLowEnergyLimit;
552  dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
553 
554  dEdxLimitBetheBloch +=
555  corrections -> ComputeIonCorrections(particle,
556  material, lowEnergyLimit / massRatio);
557  }
558 
559  G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
560  * lowEnergyLimit / scaledKineticEnergy);
561 
562  dEdx = betheBlochModel -> ComputeDEDXPerVolume(
563  material, genericIon,
564  scaledKineticEnergy, cutEnergy);
565 
566  dEdx *= chargeSquare;
567 
568  if(particle != genericIon) {
569  dEdx += corrections -> ComputeIonCorrections(particle,
570  material, kineticEnergy);
571  }
572 
573  dEdx *= factor;
574  }
575 
576  }
577 
578  if (dEdx < 0.0) dEdx = 0.0;
579 
580  return dEdx;
581 }
582 
583 // #########################################################################
584 
586  const G4ParticleDefinition* particle, // Projectile (ion)
587  const G4Material* material, // Absorber material
588  G4double lowerBoundary, // Minimum energy per nucleon
589  G4double upperBoundary, // Maximum energy per nucleon
590  G4int numBins, // Number of bins
591  G4bool logScaleEnergy) { // Logarithmic scaling of energy
592 
593  G4double atomicMassNumber = particle -> GetAtomicMass();
594  G4double materialDensity = material -> GetDensity();
595 
596  G4cout << "# dE/dx table for " << particle -> GetParticleName()
597  << " in material " << material -> GetName()
598  << " of density " << materialDensity / g * cm3
599  << " g/cm3"
600  << G4endl
601  << "# Projectile mass number A1 = " << atomicMassNumber
602  << G4endl
603  << "# ------------------------------------------------------"
604  << G4endl;
605  G4cout << "#"
606  << std::setw(13) << std::right << "E"
607  << std::setw(14) << "E/A1"
608  << std::setw(14) << "dE/dx"
609  << std::setw(14) << "1/rho*dE/dx"
610  << G4endl;
611  G4cout << "#"
612  << std::setw(13) << std::right << "(MeV)"
613  << std::setw(14) << "(MeV)"
614  << std::setw(14) << "(MeV/cm)"
615  << std::setw(14) << "(MeV*cm2/mg)"
616  << G4endl
617  << "# ------------------------------------------------------"
618  << G4endl;
619 
620  G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
621  G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
622 
623  if(logScaleEnergy) {
624 
625  energyLowerBoundary = std::log(energyLowerBoundary);
626  energyUpperBoundary = std::log(energyUpperBoundary);
627  }
628 
629  G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
630  G4double(nmbBins);
631 
632  for(int i = 0; i < numBins + 1; i++) {
633 
634  G4double energy = energyLowerBoundary + i * deltaEnergy;
635  if(logScaleEnergy) energy = std::exp(energy);
636 
637  G4double dedx = ComputeDEDXPerVolume(material, particle, energy, DBL_MAX);
638  G4cout.precision(6);
639  G4cout << std::setw(14) << std::right << energy / MeV
640  << std::setw(14) << energy / atomicMassNumber / MeV
641  << std::setw(14) << dedx / MeV * cm
642  << std::setw(14) << dedx / materialDensity / (MeV*cm2/(0.001*g))
643  << G4endl;
644  }
645 }
646 
647 // #########################################################################
648 
650  const G4ParticleDefinition* particle, // Projectile (ion)
651  const G4Material* material, // Absorber material
652  G4double lowerBoundary, // Minimum energy per nucleon
653  G4double upperBoundary, // Maximum energy per nucleon
654  G4int numBins, // Number of bins
655  G4bool logScaleEnergy) { // Logarithmic scaling of energy
656 
657  LossTableList::iterator iter = lossTableList.begin();
658  LossTableList::iterator iter_end = lossTableList.end();
659 
660  for(;iter != iter_end; iter++) {
661  G4bool isApplicable = (*iter) -> IsApplicable(particle, material);
662  if(isApplicable) {
663  (*iter) -> PrintDEDXTable(particle, material,
664  lowerBoundary, upperBoundary,
665  numBins,logScaleEnergy);
666  break;
667  }
668  }
669 }
670 
671 // #########################################################################
672 
674  std::vector<G4DynamicParticle*>* secondaries,
675  const G4MaterialCutsCouple*,
676  const G4DynamicParticle* particle,
677  G4double cutKinEnergySec,
678  G4double userMaxKinEnergySec) {
679 
680 
681  // ############## Sampling of secondaries #################################
682  // The probability density function (pdf) of the kinetic energy T of a
683  // secondary electron may be written as:
684  // pdf(T) = f(T) * g(T)
685  // where
686  // f(T) = (Tmax - Tcut) / (Tmax * Tcut) * (1 / T^2)
687  // g(T) = 1 - beta^2 * T / Tmax
688  // where Tmax is the maximum kinetic energy of the secondary, Tcut
689  // is the lower energy cut and beta is the kinetic energy of the
690  // projectile.
691  //
692  // Sampling of the kinetic energy of a secondary electron:
693  // 1) T0 is sampled from f(T) using the cumulated distribution function
694  // F(T) = int_Tcut^T f(T')dT'
695  // 2) T is accepted or rejected by evaluating the rejection function g(T)
696  // at the sampled energy T0 against a randomly sampled value
697  //
698  //
699  // See Geant4 physics reference manual (version 9.1), section 9.1.4
700  //
701  //
702  // Reference pdf: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
703  //
704  // (Implementation adapted from G4BraggIonModel)
705 
706  G4double rossiMaxKinEnergySec = MaxSecondaryKinEnergy(particle);
707  G4double maxKinEnergySec =
708  std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
709 
710  if(cutKinEnergySec >= maxKinEnergySec) return;
711 
712  G4double kineticEnergy = particle -> GetKineticEnergy();
713  G4ThreeVector direction = particle ->GetMomentumDirection();
714 
715  G4double energy = kineticEnergy + cacheMass;
716  G4double betaSquared = kineticEnergy *
717  (energy + cacheMass) / (energy * energy);
718 
719  G4double kinEnergySec;
720  G4double grej;
721 
722  do {
723 
724  // Sampling kinetic energy from f(T) (using F(T)):
725  G4double xi = G4UniformRand();
726  kinEnergySec = cutKinEnergySec * maxKinEnergySec /
727  (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
728 
729  // Deriving the value of the rejection function at the obtained kinetic
730  // energy:
731  grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
732 
733  if(grej > 1.0) {
734  G4cout << "G4IonParametrisedLossModel::SampleSecondary Warning: "
735  << "Majorant 1.0 < "
736  << grej << " for e= " << kinEnergySec
737  << G4endl;
738  }
739 
740  } while( G4UniformRand() >= grej );
741 
742  G4double momentumSec =
743  std::sqrt(kinEnergySec * (kinEnergySec + 2.0 * electron_mass_c2));
744 
745  G4double totMomentum = energy*std::sqrt(betaSquared);
746  G4double cost = kinEnergySec * (energy + electron_mass_c2) /
747  (momentumSec * totMomentum);
748  if(cost > 1.0) cost = 1.0;
749  G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
750 
751  G4double phi = twopi * G4UniformRand() ;
752 
753  G4ThreeVector directionSec(sint*std::cos(phi),sint*std::sin(phi), cost) ;
754  directionSec.rotateUz(direction);
755 
756  // create G4DynamicParticle object for delta ray
758  directionSec,
759  kinEnergySec);
760 
761  secondaries -> push_back(delta);
762 
763  // Change kinematics of primary particle
764  kineticEnergy -= kinEnergySec;
765  G4ThreeVector finalP = direction*totMomentum - directionSec*momentumSec;
766  finalP = finalP.unit();
767 
768  particleChangeLoss -> SetProposedKineticEnergy(kineticEnergy);
769  particleChangeLoss -> SetProposedMomentumDirection(finalP);
770 }
771 
772 // #########################################################################
773 
774 void G4IonParametrisedLossModel::UpdateRangeCache(
775  const G4ParticleDefinition* particle,
776  const G4MaterialCutsCouple* matCutsCouple) {
777 
778  // ############## Caching ##################################################
779  // If the ion-material-cut combination is covered by any native ion data
780  // parameterisation (for low energies), range vectors are computed
781 
782  if(particle == rangeCacheParticle &&
783  matCutsCouple == rangeCacheMatCutsCouple) {
784  }
785  else{
786  rangeCacheParticle = particle;
787  rangeCacheMatCutsCouple = matCutsCouple;
788 
789  const G4Material* material = matCutsCouple -> GetMaterial();
790  LossTableList::iterator iter = IsApplicable(particle, material);
791 
792  // If any table is applicable, the transition factor is computed:
793  if(iter != lossTableList.end()) {
794 
795  // Build range-energy and energy-range vectors if they don't exist
796  IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
797  RangeEnergyTable::iterator iterRange = r.find(ionMatCouple);
798 
799  if(iterRange == r.end()) BuildRangeVector(particle, matCutsCouple);
800 
801  rangeCacheEnergyRange = E[ionMatCouple];
802  rangeCacheRangeEnergy = r[ionMatCouple];
803  }
804  else {
805  rangeCacheEnergyRange = 0;
806  rangeCacheRangeEnergy = 0;
807  }
808  }
809 }
810 
811 // #########################################################################
812 
813 void G4IonParametrisedLossModel::UpdateDEDXCache(
814  const G4ParticleDefinition* particle,
815  const G4Material* material,
816  G4double cutEnergy) {
817 
818  // ############## Caching ##################################################
819  // If the ion-material combination is covered by any native ion data
820  // parameterisation (for low energies), a transition factor is computed
821  // which is applied to Bethe-Bloch results at higher energies to
822  // guarantee a smooth transition.
823  // This factor only needs to be calculated for the first step an ion
824  // performs inside a certain material.
825 
826  if(particle == dedxCacheParticle &&
827  material == dedxCacheMaterial &&
828  cutEnergy == dedxCacheEnergyCut) {
829  }
830  else {
831 
832  dedxCacheParticle = particle;
833  dedxCacheMaterial = material;
834  dedxCacheEnergyCut = cutEnergy;
835 
836  G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
837  dedxCacheGenIonMassRatio = massRatio;
838 
839  LossTableList::iterator iter = IsApplicable(particle, material);
840  dedxCacheIter = iter;
841 
842  // If any table is applicable, the transition factor is computed:
843  if(iter != lossTableList.end()) {
844 
845  // Retrieving the transition energy from the parameterisation table
846  G4double transitionEnergy =
847  (*iter) -> GetUpperEnergyEdge(particle, material);
848  dedxCacheTransitionEnergy = transitionEnergy;
849 
850  // Computing dE/dx from low-energy parameterisation at
851  // transition energy
852  G4double dEdxParam = (*iter) -> GetDEDX(particle, material,
853  transitionEnergy);
854 
855  G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material,
856  particle,
857  transitionEnergy,
858  cutEnergy);
859  dEdxParam -= dEdxDeltaRays;
860 
861  // Computing dE/dx from Bethe-Bloch formula at transition
862  // energy
863  G4double transitionChargeSquare =
864  GetChargeSquareRatio(particle, material, transitionEnergy);
865 
866  G4double scaledTransitionEnergy = transitionEnergy * massRatio;
867 
868  G4double dEdxBetheBloch =
869  betheBlochModel -> ComputeDEDXPerVolume(
870  material, genericIon,
871  scaledTransitionEnergy, cutEnergy);
872  dEdxBetheBloch *= transitionChargeSquare;
873 
874  // Additionally, high order corrections are added
875  dEdxBetheBloch +=
876  corrections -> ComputeIonCorrections(particle,
877  material, transitionEnergy);
878 
879  // Computing transition factor from both dE/dx values
880  dedxCacheTransitionFactor =
881  (dEdxParam - dEdxBetheBloch)/dEdxBetheBloch
882  * transitionEnergy;
883  }
884  else {
885 
886  dedxCacheParticle = particle;
887  dedxCacheMaterial = material;
888  dedxCacheEnergyCut = cutEnergy;
889 
890  dedxCacheGenIonMassRatio =
891  genericIonPDGMass / particle -> GetPDGMass();
892 
893  dedxCacheTransitionEnergy = 0.0;
894  dedxCacheTransitionFactor = 0.0;
895  }
896  }
897 }
898 
899 // #########################################################################
900 
902  const G4MaterialCutsCouple* couple,
903  const G4DynamicParticle* dynamicParticle,
904  G4double& eloss,
905  G4double&,
906  G4double length) {
907 
908  // ############## Corrections for along step energy loss calculation ######
909  // The computed energy loss (due to electronic stopping) is overwritten
910  // by this function if an ion data parameterization is available for the
911  // current ion-material pair.
912  // No action on the energy loss (due to electronic stopping) is performed
913  // if no parameterization is available. In this case the original
914  // generic ion tables (in combination with the effective charge) are used
915  // in the along step DoIt function.
916  //
917  //
918  // (Implementation partly adapted from G4BraggIonModel/G4BetheBlochModel)
919 
920  const G4ParticleDefinition* particle = dynamicParticle -> GetDefinition();
921  const G4Material* material = couple -> GetMaterial();
922 
923  G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
924 
925  if(kineticEnergy == eloss) { return; }
926 
927  G4double cutEnergy = DBL_MAX;
928  size_t cutIndex = couple -> GetIndex();
929  cutEnergy = cutEnergies[cutIndex];
930 
931  UpdateDEDXCache(particle, material, cutEnergy);
932 
933  LossTableList::iterator iter = dedxCacheIter;
934 
935  // If parameterization for ions is available the electronic energy loss
936  // is overwritten
937  if(iter != lossTableList.end()) {
938 
939  // The energy loss is calculated using the ComputeDEDXPerVolume function
940  // and the step length (it is assumed that dE/dx does not change
941  // considerably along the step)
942  eloss =
943  length * ComputeDEDXPerVolume(material, particle,
944  kineticEnergy, cutEnergy);
945 
946 #ifdef PRINT_DEBUG
947  G4cout.precision(6);
948  G4cout << "########################################################"
949  << G4endl
950  << "# G4IonParametrisedLossModel::CorrectionsAlongStep"
951  << G4endl
952  << "# cut(MeV) = " << cutEnergy/MeV
953  << G4endl;
954 
955  G4cout << "#"
956  << std::setw(13) << std::right << "E(MeV)"
957  << std::setw(14) << "l(um)"
958  << std::setw(14) << "l*dE/dx(MeV)"
959  << std::setw(14) << "(l*dE/dx)/E"
960  << G4endl
961  << "# ------------------------------------------------------"
962  << G4endl;
963 
964  G4cout << std::setw(14) << std::right << kineticEnergy / MeV
965  << std::setw(14) << length / um
966  << std::setw(14) << eloss / MeV
967  << std::setw(14) << eloss / kineticEnergy * 100.0
968  << G4endl;
969 #endif
970 
971  // If the energy loss exceeds a certain fraction of the kinetic energy
972  // (the fraction is indicated by the parameter "energyLossLimit") then
973  // the range tables are used to derive a more accurate value of the
974  // energy loss
975  if(eloss > energyLossLimit * kineticEnergy) {
976 
977  eloss = ComputeLossForStep(couple, particle,
978  kineticEnergy,length);
979 
980 #ifdef PRINT_DEBUG
981  G4cout << "# Correction applied:"
982  << G4endl;
983 
984  G4cout << std::setw(14) << std::right << kineticEnergy / MeV
985  << std::setw(14) << length / um
986  << std::setw(14) << eloss / MeV
987  << std::setw(14) << eloss / kineticEnergy * 100.0
988  << G4endl;
989 #endif
990 
991  }
992  }
993 
994  // For all corrections below a kinetic energy between the Pre- and
995  // Post-step energy values is used
996  G4double energy = kineticEnergy - eloss * 0.5;
997  if(energy < 0.0) energy = kineticEnergy * 0.5;
998 
999  G4double chargeSquareRatio = corrections ->
1000  EffectiveChargeSquareRatio(particle,
1001  material,
1002  energy);
1003  GetModelOfFluctuations() -> SetParticleAndCharge(particle,
1004  chargeSquareRatio);
1005 
1006  // A correction is applied considering the change of the effective charge
1007  // along the step (the parameter "corrFactor" refers to the effective
1008  // charge at the beginning of the step). Note: the correction is not
1009  // applied for energy loss values deriving directly from parameterized
1010  // ion stopping power tables
1011  G4double transitionEnergy = dedxCacheTransitionEnergy;
1012 
1013  if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
1014  chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1015  material,
1016  energy);
1017 
1018  G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1019  eloss *= chargeSquareRatioCorr;
1020  }
1021  else if (iter == lossTableList.end()) {
1022 
1023  chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1024  material,
1025  energy);
1026 
1027  G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1028  eloss *= chargeSquareRatioCorr;
1029  }
1030 
1031  // Ion high order corrections are applied if the current model does not
1032  // overwrite the energy loss (i.e. when the effective charge approach is
1033  // used)
1034  if(iter == lossTableList.end()) {
1035 
1036  G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
1037  G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
1038 
1039  // Corrections are only applied in the Bethe-Bloch energy region
1040  if(scaledKineticEnergy > lowEnergyLimit)
1041  eloss += length *
1042  corrections -> IonHighOrderCorrections(particle, couple, energy);
1043  }
1044 }
1045 
1046 // #########################################################################
1047 
1048 void G4IonParametrisedLossModel::BuildRangeVector(
1049  const G4ParticleDefinition* particle,
1050  const G4MaterialCutsCouple* matCutsCouple) {
1051 
1052  G4double cutEnergy = DBL_MAX;
1053  size_t cutIndex = matCutsCouple -> GetIndex();
1054  cutEnergy = cutEnergies[cutIndex];
1055 
1056  const G4Material* material = matCutsCouple -> GetMaterial();
1057 
1058  G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
1059 
1060  G4double lowerEnergy = lowerEnergyEdgeIntegr / massRatio;
1061  G4double upperEnergy = upperEnergyEdgeIntegr / massRatio;
1062 
1063  G4double logLowerEnergyEdge = std::log(lowerEnergy);
1064  G4double logUpperEnergyEdge = std::log(upperEnergy);
1065 
1066  G4double logDeltaEnergy = (logUpperEnergyEdge - logLowerEnergyEdge) /
1067  G4double(nmbBins);
1068 
1069  G4double logDeltaIntegr = logDeltaEnergy / G4double(nmbSubBins);
1070 
1071  G4LPhysicsFreeVector* energyRangeVector =
1072  new G4LPhysicsFreeVector(nmbBins+1,
1073  lowerEnergy,
1074  upperEnergy);
1075 
1076  G4double dedxLow = ComputeDEDXPerVolume(material,
1077  particle,
1078  lowerEnergy,
1079  cutEnergy);
1080 
1081  G4double range = 2.0 * lowerEnergy / dedxLow;
1082 
1083  energyRangeVector -> PutValues(0, lowerEnergy, range);
1084 
1085  G4double logEnergy = std::log(lowerEnergy);
1086  for(size_t i = 1; i < nmbBins+1; i++) {
1087 
1088  G4double logEnergyIntegr = logEnergy;
1089 
1090  for(size_t j = 0; j < nmbSubBins; j++) {
1091 
1092  G4double binLowerBoundary = std::exp(logEnergyIntegr);
1093  logEnergyIntegr += logDeltaIntegr;
1094 
1095  G4double binUpperBoundary = std::exp(logEnergyIntegr);
1096  G4double deltaIntegr = binUpperBoundary - binLowerBoundary;
1097 
1098  G4double energyIntegr = binLowerBoundary + 0.5 * deltaIntegr;
1099 
1100  G4double dedxValue = ComputeDEDXPerVolume(material,
1101  particle,
1102  energyIntegr,
1103  cutEnergy);
1104 
1105  if(dedxValue > 0.0) range += deltaIntegr / dedxValue;
1106 
1107 #ifdef PRINT_DEBUG_DETAILS
1108  G4cout << " E = "<< energyIntegr/MeV
1109  << " MeV -> dE = " << deltaIntegr/MeV
1110  << " MeV -> dE/dx = " << dedxValue/MeV*mm
1111  << " MeV/mm -> dE/(dE/dx) = " << deltaIntegr /
1112  dedxValue / mm
1113  << " mm -> range = " << range / mm
1114  << " mm " << G4endl;
1115 #endif
1116  }
1117 
1118  logEnergy += logDeltaEnergy;
1119 
1120  G4double energy = std::exp(logEnergy);
1121 
1122  energyRangeVector -> PutValues(i, energy, range);
1123 
1124 #ifdef PRINT_DEBUG_DETAILS
1125  G4cout << "G4IonParametrisedLossModel::BuildRangeVector() bin = "
1126  << i <<", E = "
1127  << energy / MeV << " MeV, R = "
1128  << range / mm << " mm"
1129  << G4endl;
1130 #endif
1131 
1132  }
1133  energyRangeVector -> SetSpline(true);
1134 
1135  G4double lowerRangeEdge =
1136  energyRangeVector -> Value(lowerEnergy);
1137  G4double upperRangeEdge =
1138  energyRangeVector -> Value(upperEnergy);
1139 
1140  G4LPhysicsFreeVector* rangeEnergyVector
1141  = new G4LPhysicsFreeVector(nmbBins+1,
1142  lowerRangeEdge,
1143  upperRangeEdge);
1144 
1145  for(size_t i = 0; i < nmbBins+1; i++) {
1146  G4double energy = energyRangeVector -> Energy(i);
1147  rangeEnergyVector ->
1148  PutValues(i, energyRangeVector -> Value(energy), energy);
1149  }
1150 
1151  rangeEnergyVector -> SetSpline(true);
1152 
1153 #ifdef PRINT_DEBUG_TABLES
1154  G4cout << *energyLossVector
1155  << *energyRangeVector
1156  << *rangeEnergyVector << G4endl;
1157 #endif
1158 
1159  IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
1160 
1161  E[ionMatCouple] = energyRangeVector;
1162  r[ionMatCouple] = rangeEnergyVector;
1163 }
1164 
1165 // #########################################################################
1166 
1168  const G4MaterialCutsCouple* matCutsCouple,
1169  const G4ParticleDefinition* particle,
1170  G4double kineticEnergy,
1171  G4double stepLength) {
1172 
1173  G4double loss = 0.0;
1174 
1175  UpdateRangeCache(particle, matCutsCouple);
1176 
1177  G4PhysicsVector* energyRange = rangeCacheEnergyRange;
1178  G4PhysicsVector* rangeEnergy = rangeCacheRangeEnergy;
1179 
1180  if(energyRange != 0 && rangeEnergy != 0) {
1181 
1182  G4double lowerEnEdge = energyRange -> Energy( 0 );
1183  G4double lowerRangeEdge = rangeEnergy -> Energy( 0 );
1184 
1185  // Computing range for pre-step kinetic energy:
1186  G4double range = energyRange -> Value(kineticEnergy);
1187 
1188  // Energy below vector boundary:
1189  if(kineticEnergy < lowerEnEdge) {
1190 
1191  range = energyRange -> Value(lowerEnEdge);
1192  range *= std::sqrt(kineticEnergy / lowerEnEdge);
1193  }
1194 
1195 #ifdef PRINT_DEBUG
1196  G4cout << "G4IonParametrisedLossModel::ComputeLossForStep() range = "
1197  << range / mm << " mm, step = " << stepLength / mm << " mm"
1198  << G4endl;
1199 #endif
1200 
1201  // Remaining range:
1202  G4double remRange = range - stepLength;
1203 
1204  // If range is smaller than step length, the loss is set to kinetic
1205  // energy
1206  if(remRange < 0.0) loss = kineticEnergy;
1207  else if(remRange < lowerRangeEdge) {
1208 
1209  G4double ratio = remRange / lowerRangeEdge;
1210  loss = kineticEnergy - ratio * ratio * lowerEnEdge;
1211  }
1212  else {
1213 
1214  G4double energy = rangeEnergy -> Value(range - stepLength);
1215  loss = kineticEnergy - energy;
1216  }
1217  }
1218 
1219  if(loss < 0.0) loss = 0.0;
1220 
1221  return loss;
1222 }
1223 
1224 // #########################################################################
1225 
1227  const G4String& nam,
1228  G4VIonDEDXTable* table,
1229  G4VIonDEDXScalingAlgorithm* algorithm) {
1230 
1231  if(table == 0) {
1232  G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1233  << " add table: Invalid pointer."
1234  << G4endl;
1235 
1236  return false;
1237  }
1238 
1239  // Checking uniqueness of name
1240  LossTableList::iterator iter = lossTableList.begin();
1241  LossTableList::iterator iter_end = lossTableList.end();
1242 
1243  for(;iter != iter_end; iter++) {
1244  G4String tableName = (*iter) -> GetName();
1245 
1246  if(tableName == nam) {
1247  G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1248  << " add table: Name already exists."
1249  << G4endl;
1250 
1251  return false;
1252  }
1253  }
1254 
1255  G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm;
1256  if(scalingAlgorithm == 0)
1257  scalingAlgorithm = new G4VIonDEDXScalingAlgorithm;
1258 
1259  G4IonDEDXHandler* handler =
1260  new G4IonDEDXHandler(table, scalingAlgorithm, nam);
1261 
1262  lossTableList.push_front(handler);
1263 
1264  return true;
1265 }
1266 
1267 // #########################################################################
1268 
1270  const G4String& nam) {
1271 
1272  LossTableList::iterator iter = lossTableList.begin();
1273  LossTableList::iterator iter_end = lossTableList.end();
1274 
1275  for(;iter != iter_end; iter++) {
1276  G4String tableName = (*iter) -> GetName();
1277 
1278  if(tableName == nam) {
1279  delete (*iter);
1280 
1281  // Remove from table list
1282  lossTableList.erase(iter);
1283 
1284  // Range vs energy and energy vs range vectors are cleared
1285  RangeEnergyTable::iterator iterRange = r.begin();
1286  RangeEnergyTable::iterator iterRange_end = r.end();
1287 
1288  for(;iterRange != iterRange_end; iterRange++)
1289  delete iterRange -> second;
1290  r.clear();
1291 
1292  EnergyRangeTable::iterator iterEnergy = E.begin();
1293  EnergyRangeTable::iterator iterEnergy_end = E.end();
1294 
1295  for(;iterEnergy != iterEnergy_end; iterEnergy++)
1296  delete iterEnergy -> second;
1297  E.clear();
1298 
1299  return true;
1300  }
1301  }
1302 
1303  return false;
1304 }
1305 
1306 // #########################################################################
1307 
1309 
1310  RemoveDEDXTable("ICRU73");
1311  AddDEDXTable("ICRU73", new G4IonStoppingData("ion_stopping_data/icru73"));
1312 }
1313 
1314 // #########################################################################
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:448
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:107
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
static G4GenericIon * Definition()
Definition: G4GenericIon.cc:49
static G4Electron * Definition()
Definition: G4Electron.cc:49
G4bool AddDEDXTable(const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:571
int G4int
Definition: G4Types.hh:78
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
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:87
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
const G4ParticleDefinition const G4Material *G4double range
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=0)
Definition: G4VEmModel.cc:387
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
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
float electron_mass_c2
Definition: hepunit.py:274
G4double DeltaRayMeanEnergyTransferRate(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
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)
std::pair< const G4ParticleDefinition *, const G4MaterialCutsCouple * > IonMatCouple
Hep3Vector unit() const
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)
#define G4endl
Definition: G4ios.hh:61
const G4String & GetName() const
Definition: G4VEmModel.hh:753
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
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:346
G4GLOB_DLL std::ostream G4cerr
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &, G4double &, G4double)