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