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

Last change on this file since 1014 was 1007, checked in by garnier, 17 years ago

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