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

Last change on this file since 1199 was 1192, checked in by garnier, 16 years ago

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