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

Last change on this file since 1066 was 1055, checked in by garnier, 17 years ago

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