| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // ===========================================
28 // GEANT4 class source file
29 //
30 // Class: G4IonParametrisedLoss
31 //
32 // Base class: G4VEmModel (utils)
33 //
34 // Author: Anton Lechner (Anton.
35 //
36 // First implementation: 10. 11. 2008
37 //
38 // Modifications: 03. 02. 2009 - Bug fix itera
39 // 11. 03. 2009 - Introduced ne
40 // and modified
41 // (tables are n
42 // Minor bug fix
43 // 11. 05. 2009 - Introduced sc
44 // G4IonDEDXScal
45 // 12. 11. 2009 - Moved from or
46 // class (G4IonS
47 // of reading st
48 // in G4LEDATA (
49 // Simultanesoul
50 // ICRU 73 is in
51 // - Removed nucle
52 // AlongStep sin
53 // - Added functio
54 // of heavy ions
55 // - Minor fix in
56 // - Minor fix in
57 // 23. 11. 2009 - Changed energ
58 // to improve ac
59 // 24. 11. 2009 - Bug fix: Rang
60 // materials app
61 // regions (adde
62 // modified Buil
63 // functions acc
64 // - Removed GetRa
65 // 04. 11. 2010 - Moved virtual
66 //
67 //
68 // Class description:
69 // Model for computing the energy loss of i
70 // parameterisation of dE/dx tables (by def
71 // ion-material combinations and/or project
72 // by this model, the G4BraggIonModel and G
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::G4IonParametrisedL
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()
118
119 // The Bragg ion and Bethe Bloch models are
120 braggIonModel = new G4BraggIonModel();
121 betheBlochModel = new G4BetheBlochModel();
122
123 // The boundaries for the range tables are s
124 lowerEnergyEdgeIntegr = 0.025 * MeV;
125 upperEnergyEdgeIntegr = betheBlochModel -> H
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::~G4IonParametrised
155
156 // dE/dx table objects are deleted and the c
157 LossTableList::iterator iterTables = lossTab
158 LossTableList::iterator iterTables_end = los
159
160 for(;iterTables != iterTables_end; ++iterTab
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->se
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->se
173 E.clear();
174
175 }
176
177 // ###########################################
178
179 G4double G4IonParametrisedLossModel::MinEnergy
180 const G
181 const G
182
183 return couple->GetMaterial()->GetIonisation(
184 }
185
186 // ###########################################
187
188 G4double G4IonParametrisedLossModel::MaxSecond
189 const G4ParticleD
190 G4double kineticE
191
192 // ############## Maximum energy of secondar
193 // Function computes maximum energy of secon
194 // released by an ion
195 //
196 // See Geant4 physics reference manual (vers
197 //
198 // Ref.: W.M. Yao et al, Jour. of Phys. G 33
199 // C.Caso et al. (Part. Data Group), E
200 // B. Rossi, High energy particles, Ne
201 //
202 // (Implementation adapted from G4BraggIonMo
203
204 if(particle != cacheParticle) UpdateCache(pa
205
206 G4double tau = kineticEnergy/cacheMass;
207 G4double tmax = 2.0 * electron_mass_c2 * tau
208 (1. + 2.0 * (tau + 1.) * cac
209 cacheElecMassRatio * cacheEl
210
211 return tmax;
212 }
213
214 // ###########################################
215
216 G4double G4IonParametrisedLossModel::GetCharge
217 const G4ParticleD
218 const G4Material* material,
219 G4double kineticE
220
221 G4double chargeSquareRatio = corrections ->
222 Effective
223
224
225 corrFactor = chargeSquareRatio *
226 corrections -> Effectiv
227
228
229 return corrFactor;
230 }
231
232 // ###########################################
233
234 G4double G4IonParametrisedLossModel::GetPartic
235 const G4ParticleD
236 const G4Material* material,
237 G4double kineticE
238
239 return corrections -> GetParticleCharge(part
240 }
241
242 // ###########################################
243
244 void G4IonParametrisedLossModel::Initialise(
245 const G
246 const G
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
270 if(!isInitialised) {
271 G4bool icru90 = G4EmParameters::Instance()
272 isInitialised = true;
273 AddDEDXTable("ICRU73",
274 new G4IonStoppingData("ion_stopping_data/
275 new G4IonDEDXScalingICRU73());
276 }
277 // The cache of loss tables is cleared
278 LossTableList::iterator iterTables = lossTab
279 LossTableList::iterator iterTables_end = los
280
281 for(;iterTables != iterTables_end; ++iterTab
282 (*iterTables) -> ClearCache();
283 }
284
285 // Range vs energy and energy vs range vecto
286 // cleared
287 RangeEnergyTable::iterator iterRange = r.beg
288 RangeEnergyTable::iterator iterRange_end = r
289
290 for(;iterRange != iterRange_end; ++iterRange
291 delete iterRange->second;
292 }
293 r.clear();
294
295 EnergyRangeTable::iterator iterEnergy = E.be
296 EnergyRangeTable::iterator iterEnergy_end =
297
298 for(;iterEnergy != iterEnergy_end; ++iterEne
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::Ge
309 G4int nmbCouples = (G4int)coupleTable->GetTa
310
311 #ifdef PRINT_TABLE_BUILT
312 G4cout << "G4IonParametrisedLossModel::Ini
313 << " Building dE/dx vectors:"
314 << G4endl;
315 #endif
316
317 for (G4int i = 0; i < nmbCouples; ++i) {
318
319 const G4MaterialCutsCouple* couple = coupl
320 const G4Material* material = couple->GetMa
321
322 for(G4int atomicNumberIon = 3; atomicNumbe
323
324 LossTableList::iterator iter = lossTabl
325 LossTableList::iterator iter_end = loss
326
327 for(;iter != iter_end; ++iter) {
328
329 if(*iter == 0) {
330 G4cout << "G4IonParametrisedLoss
331 << " Skipping illegal tab
332 << G4endl;
333 }
334
335 if((*iter)->BuildDEDXTable(atomicNumberIon
336
337 #ifdef PRINT_TABLE_BUILT
338 G4cout << " Atomic Number Ion =
339 << ", Material = " << mate
340 << ", Table = " << (*iter)
341 << G4endl;
342 #endif
343 break;
344 }
345 }
346 }
347 }
348
349 // The particle change object
350 if(! particleChangeLoss) {
351 particleChangeLoss = GetParticleChangeForL
352 braggIonModel->SetParticleChange(particleC
353 betheBlochModel->SetParticleChange(particl
354 }
355
356 // The G4BraggIonModel and G4BetheBlochModel
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::ComputeCr
365 const G4Par
366 G4double ki
367 G4double atomicNumber,
368 G4double,
369 G4double cu
370 G4double ma
371
372 // ############## Cross section per atom ##
373 // Function computes ionization cross sectio
374 //
375 // See Geant4 physics reference manual (vers
376 //
377 // Ref.: W.M. Yao et al, Jour. of Phys. G 33
378 // B. Rossi, High energy particles, Ne
379 //
380 // (Implementation adapted from G4BraggIonMo
381
382 G4double crosssection = 0.0;
383 G4double tmax = MaxSecondaryEnergy(part
384 G4double maxEnergy = std::min(tmax, maxKinEn
385
386 if(cutEnergy < tmax) {
387
388 G4double energy = kineticEnergy + cacheM
389 G4double betaSquared = kineticEnergy *
390 (energy +
391
392 crosssection = 1.0 / cutEnergy - 1.0 / ma
393 betaSquared * std::lo
394
395 crosssection *= twopi_mc2_rcl2 * cacheCha
396 }
397
398 #ifdef PRINT_DEBUG_CS
399 G4cout << "#################################
400 << G4endl
401 << "# G4IonParametrisedLossModel::Com
402 << G4endl
403 << "# particle =" << particle -> GetP
404 << G4endl
405 << "# cut(MeV) = " << cutEnergy/MeV
406 << G4endl;
407
408 G4cout << "#"
409 << std::setw(13) << std::right << "E(
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 << kin
417 << std::setw(14) << crosssection / (u
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::CrossSect
430 const G4Material* material,
431 const G4Par
432 G4double ki
433 G4double cu
434 G4double ma
435
436 G4double nbElecPerVolume = material -> GetTo
437 G4double cross = ComputeCrossSectionPerAtom(
438
439
440
441
442
443 return cross;
444 }
445
446 // ###########################################
447
448 G4double G4IonParametrisedLossModel::ComputeDE
449 const G4
450 const G4ParticleDefinition
451 G4double kineticEnergy,
452 G4double cutEnergy) {
453
454 // ############## dE/dx ####################
455 // Function computes dE/dx values, where fol
456 // A. If the ion-material pair is covered
457 // parameterisation, then:
458 // * This parameterization is used for
459 // limit,
460 // * whereas above the limit the Bethe-
461 // combination with an effective char
462 // correction terms.
463 // A smoothing procedure is applied to
464 // the second approach. The smoothing f
465 // values of both approaches at the tra
466 // correction terms are included in the
467 // factor).
468 // B. If the particle is a generic ion, th
469 // models are used and a smoothing proc
470 // obtained with the second approach.
471 // C. If the ion-material is not covered b
472 // then:
473 // * The BraggIon model is used for ene
474 // limit,
475 // * whereas above the limit the Bethe-
476 // combination with an effective char
477 // correction terms.
478 // Also in this case, a smoothing proce
479 // computed with the second model.
480
481 G4double dEdx = 0.0;
482 UpdateDEDXCache(particle, material, cutEnerg
483
484 LossTableList::iterator iter = dedxCacheIter
485
486 if(iter != lossTableList.end()) {
487
488 G4double transitionEnergy = dedxCacheTran
489
490 if(transitionEnergy > kineticEnergy) {
491
492 dEdx = (*iter) -> GetDEDX(particle, ma
493
494 G4double dEdxDeltaRays = DeltaRayMeanE
495 par
496 kin
497 cut
498 dEdx -= dEdxDeltaRays;
499 }
500 else {
501 G4double massRatio = dedxCacheGenIonMa
502
503 G4double chargeSquare =
504 GetChargeSquareRatio(pa
505
506 G4double scaledKineticEnergy = kinetic
507 G4double scaledTransitionEnergy = tran
508
509 G4double lowEnergyLimit = betheBlochMo
510
511 if(scaledTransitionEnergy >= lowEnergy
512
513 dEdx = betheBlochModel -> ComputeDE
514 material
515 scaledKineticEnergy, cutEner
516
517 dEdx *= chargeSquare;
518
519 dEdx += corrections -> ComputeIonCo
520
521
522 G4double factor = 1.0 + dedxCacheTr
523
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(pa
536 massRatio = genericIonPDGMass / partic
537 }
538
539 G4double scaledKineticEnergy = kineticEne
540
541 G4double lowEnergyLimit = betheBlochModel
542 if(scaledKineticEnergy < lowEnergyLimit)
543 dEdx = braggIonModel -> ComputeDEDXPer
544 material
545 scaledKineticEnergy, cutEner
546
547 dEdx *= chargeSquare;
548 }
549 else {
550 G4double dEdxLimitParam = braggIonMode
551 material
552 lowEnergyLimit, cutEnergy);
553
554 G4double dEdxLimitBetheBloch = betheBl
555 material
556 lowEnergyLimit, cutEnergy);
557
558 if(particle != genericIon) {
559 G4double chargeSquareLowEnergyLimit
560 GetChargeSquareRatio(particle,
561 lowEnergyL
562
563 dEdxLimitParam *= chargeSquareLowEn
564 dEdxLimitBetheBloch *= chargeSquare
565
566 dEdxLimitBetheBloch +=
567 corrections -> ComputeIonC
568 material
569 }
570
571 G4double factor = (1.0 + (dEdxLimitPar
572 * lowEnergyLimi
573
574 dEdx = betheBlochModel -> ComputeDEDXP
575 material
576 scaledKineticEnergy, cutEner
577
578 dEdx *= chargeSquare;
579
580 if(particle != genericIon) {
581 dEdx += corrections -> ComputeIonCo
582 material
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::PrintDEDXTabl
598 const G4ParticleDefinition*
599 const G4Material* material,
600 G4double lowerBoundary,
601 G4double upperBoundary,
602 G4int numBins,
603 G4bool logScaleEnergy) {
604
605 G4double atomicMassNumber = particle -> GetA
606 G4double materialDensity = material -> GetDe
607
608 G4cout << "# dE/dx table for " << particle -
609 << " in material " << material -> Get
610 << " of density " << materialDensity
611 << " g/cm3"
612 << G4endl
613 << "# Projectile mass number A1 = " <
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 << "(M
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
633 G4double energyUpperBoundary = upperBoundary
634
635 if(logScaleEnergy) {
636
637 energyLowerBoundary = std::log(energyLowe
638 energyUpperBoundary = std::log(energyUppe
639 }
640
641 G4double deltaEnergy = (energyUpperBoundary
642
643
644 for(int i = 0; i < numBins + 1; i++) {
645
646 G4double energy = energyLowerBoundary +
647 if(logScaleEnergy) energy = G4Exp(energy
648
649 G4double dedx = ComputeDEDXPerVolume(mat
650 G4cout.precision(6);
651 G4cout << std::setw(14) << std::right <<
652 << std::setw(14) << energy / atom
653 << std::setw(14) << dedx / MeV *
654 << std::setw(14) << dedx / materi
655 << G4endl;
656 }
657 }
658
659 // ###########################################
660
661 void G4IonParametrisedLossModel::PrintDEDXTabl
662 const G4ParticleDefinition*
663 const G4Material* material,
664 G4double lowerBoundary,
665 G4double upperBoundary,
666 G4int numBins,
667 G4bool logScaleEnergy) {
668
669 LossTableList::iterator iter = lossTableList
670 LossTableList::iterator iter_end = lossTable
671
672 for(;iter != iter_end; iter++) {
673 G4bool isApplicable = (*iter) -> IsAppl
674 if(isApplicable) {
675 (*iter) -> PrintDEDXTable(particle, material
676 lowerBoundar
677 numBins,logS
678 break;
679 }
680 }
681 }
682
683 // ###########################################
684
685 void G4IonParametrisedLossModel::SampleSeconda
686 std::vector<G4Dyn
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
693 // secondary electron may be written as:
694 // pdf(T) = f(T) * g(T)
695 // where
696 // f(T) = (Tmax - Tcut) / (Tmax * Tcut) *
697 // g(T) = 1 - beta^2 * T / Tmax
698 // where Tmax is the maximum kinetic energy
699 // is the lower energy cut and beta is the k
700 // projectile.
701 //
702 // Sampling of the kinetic energy of a secon
703 // 1) T0 is sampled from f(T) using the cum
704 // F(T) = int_Tcut^T f(T')dT'
705 // 2) T is accepted or rejected by evaluati
706 // at the sampled energy T0 against a ra
707 //
708 //
709 // See Geant4 physics reference manual (vers
710 //
711 //
712 // Reference pdf: W.M. Yao et al, Jour. of P
713 //
714 // (Implementation adapted from G4BraggIonMo
715
716 G4double rossiMaxKinEnergySec = MaxSecondary
717 G4double maxKinEnergySec =
718 std::min(rossiMaxKinEn
719
720 if(cutKinEnergySec >= maxKinEnergySec) retur
721
722 G4double kineticEnergy = particle -> GetKine
723
724 G4double energy = kineticEnergy + cacheMass
725 G4double betaSquared = kineticEnergy * (ener
726 / (energy * energy);
727
728 G4double kinEnergySec;
729 G4double grej;
730
731 do {
732
733 // Sampling kinetic energy from f(T) (usin
734 G4double xi = G4UniformRand();
735 kinEnergySec = cutKinEnergySec * maxKinEne
736 (maxKinEnergySec * (1.
737
738 // Deriving the value of the rejection fun
739 // energy:
740 grej = 1.0 - betaSquared * kinEnergySec /
741
742 if(grej > 1.0) {
743 G4cout << "G4IonParametrisedLossModel:
744 << "Majorant 1.0 < "
745 << grej << " for e= " << kinEne
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 = G4Ele
755
756 G4DynamicParticle* delta = new G4DynamicPart
757 GetAngularDistribution()->SampleDirection(
758
759
760
761 secondaries->push_back(delta);
762
763 // Change kinematics of primary particle
764 G4ThreeVector direction = particle ->GetMome
765 G4double totalMomentum = std::sqrt(kineticEn
766
767 G4ThreeVector finalP = totalMomentum*directi
768 finalP = finalP.unit();
769
770 kineticEnergy -= kinEnergySec;
771
772 particleChangeLoss->SetProposedKineticEnergy
773 particleChangeLoss->SetProposedMomentumDirec
774 }
775
776 // ###########################################
777
778 void G4IonParametrisedLossModel::UpdateRangeCa
779 const G4ParticleDefinitio
780 const G4MaterialCutsCoupl
781
782 // ############## Caching ##################
783 // If the ion-material-cut combination is co
784 // parameterisation (for low energies), rang
785
786 if(particle == rangeCacheParticle &&
787 matCutsCouple == rangeCacheMatCutsCouple)
788 }
789 else{
790 rangeCacheParticle = particle;
791 rangeCacheMatCutsCouple = matCutsCouple;
792
793 const G4Material* material = matCutsCoupl
794 LossTableList::iterator iter = IsApplicab
795
796 // If any table is applicable, the transi
797 if(iter != lossTableList.end()) {
798
799 // Build range-energy and energy-range
800 IonMatCouple ionMatCouple = std::make_
801 RangeEnergyTable::iterator iterRange =
802
803 if(iterRange == r.end()) BuildRangeVec
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::UpdateDEDXCac
818 const G4ParticleDefinitio
819 const G4Material* materia
820 G4double cutEnergy) {
821
822 // ############## Caching ##################
823 // If the ion-material combination is covere
824 // parameterisation (for low energies), a tr
825 // which is applied to Bethe-Bloch results a
826 // guarantee a smooth transition.
827 // This factor only needs to be calculated f
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 /
841 dedxCacheGenIonMassRatio = massRatio;
842
843 LossTableList::iterator iter = IsApplicab
844 dedxCacheIter = iter;
845
846 // If any table is applicable, the transi
847 if(iter != lossTableList.end()) {
848
849 // Retrieving the transition energy fr
850 G4double transitionEnergy =
851 (*iter) -> GetUpperEnergyEdge
852 dedxCacheTransitionEnergy = transition
853
854 // Computing dE/dx from low-energy par
855 // transition energy
856 G4double dEdxParam = (*iter) -> GetDED
857 tra
858
859 G4double dEdxDeltaRays = DeltaRayMeanE
860 par
861 tra
862 cut
863 dEdxParam -= dEdxDeltaRays;
864
865 // Computing dE/dx from Bethe-Bloch fo
866 // energy
867 G4double transitionChargeSquare =
868 GetChargeSquareRatio(particle, m
869
870 G4double scaledTransitionEnergy = tran
871
872 G4double dEdxBetheBloch =
873 betheBlochModel ->
874 materi
875 scaledTransitionEnergy, cu
876 dEdxBetheBloch *= transitionChargeSqua
877
878 // Additionally, high order correction
879 dEdxBetheBloch +=
880 corrections -> ComputeIonCorrectio
881
882
883 // Computing transition factor from bo
884 dedxCacheTransitionFactor =
885 (dEdxParam - dEdxBetheBlo
886 * transitionEnerg
887 }
888 else {
889
890 dedxCacheParticle = particle;
891 dedxCacheMaterial = material;
892 dedxCacheEnergyCut = cutEnergy;
893
894 dedxCacheGenIonMassRatio =
895 genericIonPDGMass
896
897 dedxCacheTransitionEnergy = 0.0;
898 dedxCacheTransitionFactor = 0.0;
899 }
900 }
901 }
902
903 // ###########################################
904
905 void G4IonParametrisedLossModel::CorrectionsAl
906 const G4MaterialC
907 const G4DynamicParticle* dynamicPar
908 const G4double& l
909 G4double& eloss) {
910
911 // ############## Corrections for along step
912 // The computed energy loss (due to electron
913 // by this function if an ion data parameter
914 // current ion-material pair.
915 // No action on the energy loss (due to elec
916 // if no parameterization is available. In t
917 // generic ion tables (in combination with t
918 // in the along step DoIt function.
919 //
920 //
921 // (Implementation partly adapted from G4Bra
922
923 const G4ParticleDefinition* particle = dynam
924 const G4Material* material = couple -> GetMa
925
926 G4double kineticEnergy = dynamicParticle ->
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, cutEnerg
935
936 LossTableList::iterator iter = dedxCacheIter
937
938 // If parameterization for ions is available
939 // is overwritten
940 if(iter != lossTableList.end()) {
941
942 // The energy loss is calculated using th
943 // and the step length (it is assumed tha
944 // considerably along the step)
945 eloss =
946 length * ComputeDEDXPerVolume(material
947 kineticE
948
949 #ifdef PRINT_DEBUG
950 G4cout.precision(6);
951 G4cout << "#################################
952 << G4endl
953 << "# G4IonParametrisedLossModel::Cor
954 << G4endl
955 << "# cut(MeV) = " << cutEnergy/MeV
956 << G4endl;
957
958 G4cout << "#"
959 << std::setw(13) << std::right << "E(
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 << kin
968 << std::setw(14) << length / um
969 << std::setw(14) << eloss / MeV
970 << std::setw(14) << eloss / kineticEn
971 << G4endl;
972 #endif
973
974 // If the energy loss exceeds a certain f
975 // (the fraction is indicated by the para
976 // the range tables are used to derive a
977 // energy loss
978 if(eloss > energyLossLimit * kineticEnerg
979
980 eloss = ComputeLossForStep(couple, par
981 kineticEner
982
983 #ifdef PRINT_DEBUG
984 G4cout << "# Correction applied:"
985 << G4endl;
986
987 G4cout << std::setw(14) << std::right << kin
988 << std::setw(14) << length / um
989 << std::setw(14) << eloss / MeV
990 << std::setw(14) << eloss / kineticEn
991 << G4endl;
992 #endif
993
994 }
995 }
996
997 // For all corrections below a kinetic energ
998 // Post-step energy values is used
999 G4double energy = kineticEnergy - eloss * 0.
1000 if(energy < 0.0) energy = kineticEnergy * 0
1001
1002 G4double chargeSquareRatio = corrections ->
1003 Effectiv
1004
1005
1006 GetModelOfFluctuations() -> SetParticleAndC
1007
1008
1009 // A correction is applied considering the
1010 // along the step (the parameter "corrFacto
1011 // charge at the beginning of the step). No
1012 // applied for energy loss values deriving
1013 // ion stopping power tables
1014 G4double transitionEnergy = dedxCacheTransi
1015
1016 if(iter != lossTableList.end() && transitio
1017 chargeSquareRatio *= corrections -> Effe
1018
1019
1020
1021 G4double chargeSquareRatioCorr = chargeS
1022 eloss *= chargeSquareRatioCorr;
1023 }
1024 else if (iter == lossTableList.end()) {
1025
1026 chargeSquareRatio *= corrections -> Effe
1027
1028
1029
1030 G4double chargeSquareRatioCorr = chargeS
1031 eloss *= chargeSquareRatioCorr;
1032 }
1033
1034 // Ion high order corrections are applied i
1035 // overwrite the energy loss (i.e. when the
1036 // used)
1037 if(iter == lossTableList.end()) {
1038
1039 G4double scaledKineticEnergy = kineticEn
1040 G4double lowEnergyLimit = betheBlochMode
1041
1042 // Corrections are only applied in the B
1043 if(scaledKineticEnergy > lowEnergyLimit)
1044 eloss += length *
1045 corrections -> IonHighOrderCorrec
1046 }
1047 }
1048
1049 // ##########################################
1050
1051 void G4IonParametrisedLossModel::BuildRangeVe
1052 const G4ParticleDefiniti
1053 const G4MaterialCutsCoup
1054
1055 G4double cutEnergy = DBL_MAX;
1056 std::size_t cutIndex = matCutsCouple -> Get
1057 cutEnergy = cutEnergies[cutIndex];
1058
1059 const G4Material* material = matCutsCouple
1060
1061 G4double massRatio = genericIonPDGMass / pa
1062
1063 G4double lowerEnergy = lowerEnergyEdgeInteg
1064 G4double upperEnergy = upperEnergyEdgeInteg
1065
1066 G4double logLowerEnergyEdge = std::log(lowe
1067 G4double logUpperEnergyEdge = std::log(uppe
1068
1069 G4double logDeltaEnergy = (logUpperEnergyEd
1070
1071
1072 G4double logDeltaIntegr = logDeltaEnergy /
1073
1074 G4PhysicsFreeVector* energyRangeVector =
1075 new G4PhysicsFreeVector(nmbBins+1,
1076 lowerEnergy,
1077 upperEnergy,
1078 /*spline=*/true);
1079
1080 G4double dedxLow = ComputeDEDXPerVolume(mat
1081 par
1082 low
1083 cut
1084
1085 G4double range = 2.0 * lowerEnergy / dedxLo
1086
1087 energyRangeVector -> PutValues(0, lowerEner
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;
1095
1096 G4double binLowerBoundary = G4Exp(l
1097 logEnergyIntegr += logDeltaIntegr;
1098
1099 G4double binUpperBoundary = G4Exp(l
1100 G4double deltaIntegr = binUpperBoun
1101
1102 G4double energyIntegr = binLowerBou
1103
1104 G4double dedxValue = ComputeDEDXPer
1105
1106
1107 cutEnergy);
1108
1109 if(dedxValue > 0.0) range += deltaI
1110
1111 #ifdef PRINT_DEBUG_DETAILS
1112 G4cout << " E = "<< energyInt
1113 << " MeV -> dE = " << de
1114 << " MeV -> dE/dx = " <<
1115 << " MeV/mm -> dE/(dE/dx
1116
1117 << " mm -> range = " <<
1118 << " mm " << G4endl;
1119 #endif
1120 }
1121
1122 logEnergy += logDeltaEnergy;
1123
1124 G4double energy = G4Exp(logEnergy);
1125
1126 energyRangeVector -> PutValues(i, energ
1127
1128 #ifdef PRINT_DEBUG_DETAILS
1129 G4cout << "G4IonParametrisedLossModel::
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 -> Valu
1142 G4double upperRangeEdge =
1143 energyRangeVector -> Valu
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 ->
1153 rangeEnergyVector ->
1154 PutValues(i, energyRangeVector -
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(
1166
1167 E[ionMatCouple] = energyRangeVector;
1168 r[ionMatCouple] = rangeEnergyVector;
1169 }
1170
1171 // ##########################################
1172
1173 G4double G4IonParametrisedLossModel::ComputeL
1174 const G4MaterialCutsCoup
1175 const G4ParticleDefiniti
1176 G4double kineticEnergy,
1177 G4double stepLength) {
1178
1179 G4double loss = 0.0;
1180
1181 UpdateRangeCache(particle, matCutsCouple);
1182
1183 G4PhysicsVector* energyRange = rangeCacheEn
1184 G4PhysicsVector* rangeEnergy = rangeCacheRa
1185
1186 if(energyRange != 0 && rangeEnergy != 0) {
1187
1188 G4double lowerEnEdge = energyRange -> En
1189 G4double lowerRangeEdge = rangeEnergy ->
1190
1191 // Computing range for pre-step kinetic
1192 G4double range = energyRange -> Value(ki
1193
1194 // Energy below vector boundary:
1195 if(kineticEnergy < lowerEnEdge) {
1196
1197 range = energyRange -> Value(lowerEn
1198 range *= std::sqrt(kineticEnergy / lo
1199 }
1200
1201 #ifdef PRINT_DEBUG
1202 G4cout << "G4IonParametrisedLossModel::C
1203 << range / mm << " mm, step = " <
1204 << G4endl;
1205 #endif
1206
1207 // Remaining range:
1208 G4double remRange = range - stepLength;
1209
1210 // If range is smaller than step length,
1211 // energy
1212 if(remRange < 0.0) loss = kineticEnergy;
1213 else if(remRange < lowerRangeEdge) {
1214
1215 G4double ratio = remRange / lowerRang
1216 loss = kineticEnergy - ratio * ratio
1217 }
1218 else {
1219
1220 G4double energy = rangeEnergy -> Valu
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::AddDEDXTab
1233 const G4Strin
1234 G4VIonDEDXTab
1235 G4VIonDEDXSca
1236
1237 if(table == 0) {
1238 G4cout << "G4IonParametrisedLossModel::A
1239 << " add table: Invalid pointer."
1240 << G4endl;
1241
1242 return false;
1243 }
1244
1245 // Checking uniqueness of name
1246 LossTableList::iterator iter = lossTableLis
1247 LossTableList::iterator iter_end = lossTabl
1248
1249 for(;iter != iter_end; ++iter) {
1250 const G4String& tableName = (*iter)->Get
1251
1252 if(tableName == nam) {
1253 G4cout << "G4IonParametrisedLossModel
1254 << " add table: Name already e
1255 << G4endl;
1256
1257 return false;
1258 }
1259 }
1260
1261 G4VIonDEDXScalingAlgorithm* scalingAlgorith
1262 if(scalingAlgorithm == 0)
1263 scalingAlgorithm = new G4VIonDEDXScaling
1264
1265 G4IonDEDXHandler* handler =
1266 new G4IonDEDXHandler(ta
1267
1268 lossTableList.push_front(handler);
1269
1270 return true;
1271 }
1272
1273 // ##########################################
1274
1275 G4bool G4IonParametrisedLossModel::RemoveDEDX
1276 const G4String& nam) {
1277
1278 LossTableList::iterator iter = lossTableLis
1279 LossTableList::iterator iter_end = lossTabl
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 rang
1291 RangeEnergyTable::iterator iterRange
1292 RangeEnergyTable::iterator iterRange_
1293
1294 for(;iterRange != iterRange_end; iter
1295 d
1296 r.clear();
1297
1298 EnergyRangeTable::iterator iterEnergy
1299 EnergyRangeTable::iterator iterEnergy
1300
1301 for(;iterEnergy != iterEnergy_end; it
1302 d
1303 E.clear();
1304
1305 return true;
1306 }
1307 }
1308
1309 return false;
1310 }
1311