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