source: trunk/source/processes/electromagnetic/lowenergy/src/G4IonParametrisedLossModel.cc@ 1316

Last change on this file since 1316 was 1228, checked in by garnier, 16 years ago

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