Loading...
Searching...
No Matches
G4IonParametrisedLossModel Class Reference
#include <G4IonParametrisedLossModel.hh>
Inheritance diagram for G4IonParametrisedLossModel:
Protected Member Functions | |
| virtual G4double | MaxSecondaryEnergy (const G4ParticleDefinition *, G4double) |
| Protected Member Functions inherited from G4VEmModel | |
| G4ParticleChangeForLoss * | GetParticleChangeForLoss () |
| G4ParticleChangeForGamma * | GetParticleChangeForGamma () |
| virtual G4double | MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy) |
| const G4MaterialCutsCouple * | CurrentCouple () const |
| void | SetCurrentElement (const G4Element *) |
Additional Inherited Members | |
| Protected Attributes inherited from G4VEmModel | |
| G4VParticleChange * | pParticleChange |
| G4PhysicsTable * | xSectionTable |
| const std::vector< G4double > * | theDensityFactor |
| const std::vector< G4int > * | theDensityIdx |
Detailed Description
Definition at line 93 of file G4IonParametrisedLossModel.hh.
Constructor & Destructor Documentation
◆ G4IonParametrisedLossModel()
| G4IonParametrisedLossModel::G4IonParametrisedLossModel | ( | const G4ParticleDefinition * | particle = 0, |
| const G4String & | name = "ParamICRU73" |
||
| ) |
Definition at line 103 of file G4IonParametrisedLossModel.cc.
106 : G4VEmModel(nam),
107 braggIonModel(0),
108 betheBlochModel(0),
109 nmbBins(90),
110 nmbSubBins(100),
111 particleChangeLoss(0),
112 corrFactor(1.0),
113 energyLossLimit(0.01),
114 cutEnergies(0)
115{
116 genericIon = G4GenericIon::Definition();
117 genericIonPDGMass = genericIon -> GetPDGMass();
118 corrections = G4LossTableManager::Instance() -> EmCorrections();
119
120 // The upper limit of the current model is set to 100 TeV
121 SetHighEnergyLimit(100.0 * TeV);
122
123 // The Bragg ion and Bethe Bloch models are instantiated
126
127 // By default ICRU 73 stopping power tables are loaded:
131
132 // The boundaries for the range tables are set
133 lowerEnergyEdgeIntegr = 0.025 * MeV;
134 upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit();
135
136 // Cache parameters are set
137 cacheParticle = 0;
138 cacheMass = 0;
139 cacheElecMassRatio = 0;
140 cacheChargeSquare = 0;
141
142 // Cache parameters are set
143 rangeCacheParticle = 0;
144 rangeCacheMatCutsCouple = 0;
145 rangeCacheEnergyRange = 0;
146 rangeCacheRangeEnergy = 0;
147
148 // Cache parameters are set
149 dedxCacheParticle = 0;
150 dedxCacheMaterial = 0;
151 dedxCacheEnergyCut = 0;
152 dedxCacheIter = lossTableList.end();
153 dedxCacheTransitionEnergy = 0.0;
154 dedxCacheTransitionFactor = 0.0;
155 dedxCacheGenIonMassRatio = 0.0;
156}
Definition: G4BetheBlochModel.hh:73
Definition: G4BraggIonModel.hh:69
G4bool AddDEDXTable(const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
Definition: G4IonParametrisedLossModel.cc:1246
Definition: G4IonStoppingData.hh:60
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:107
Definition: G4VEmModel.hh:102
◆ ~G4IonParametrisedLossModel()
|
virtual |
Definition at line 160 of file G4IonParametrisedLossModel.cc.
160 {
161
162 // Range vs energy table objects are deleted and the container is cleared
163 RangeEnergyTable::iterator iterRange = r.begin();
164 RangeEnergyTable::iterator iterRange_end = r.end();
165
166 for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second;
167 r.clear();
168
169 // Energy vs range table objects are deleted and the container is cleared
170 EnergyRangeTable::iterator iterEnergy = E.begin();
171 EnergyRangeTable::iterator iterEnergy_end = E.end();
172
173 for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second;
174 E.clear();
175
176 // dE/dx table objects are deleted and the container is cleared
177 LossTableList::iterator iterTables = lossTableList.begin();
178 LossTableList::iterator iterTables_end = lossTableList.end();
179
180 for(;iterTables != iterTables_end; iterTables++) delete *iterTables;
181 lossTableList.clear();
182
183 // The Bragg ion and Bethe Bloch objects are deleted
184 delete betheBlochModel;
185 delete braggIonModel;
186}
Member Function Documentation
◆ AddDEDXTable()
| G4bool G4IonParametrisedLossModel::AddDEDXTable | ( | const G4String & | name, |
| G4VIonDEDXTable * | table, | ||
| G4VIonDEDXScalingAlgorithm * | algorithm = 0 |
||
| ) |
Definition at line 1246 of file G4IonParametrisedLossModel.cc.
1249 {
1250
1251 if(table == 0) {
1253 << " add table: Invalid pointer."
1254 << G4endl;
1255
1256 return false;
1257 }
1258
1259 // Checking uniqueness of name
1260 LossTableList::iterator iter = lossTableList.begin();
1261 LossTableList::iterator iter_end = lossTableList.end();
1262
1263 for(;iter != iter_end; iter++) {
1265
1266 if(tableName == nam) {
1268 << " add table: Name already exists."
1269 << G4endl;
1270
1271 return false;
1272 }
1273 }
1274
1275 G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm;
1276 if(scalingAlgorithm == 0)
1278
1279 G4IonDEDXHandler* handler =
1281
1282 lossTableList.push_front(handler);
1283
1284 return true;
1285}
G4DLLIMPORT std::ostream G4cerr
Definition: G4IonDEDXHandler.hh:82
Definition: G4String.hh:105
Referenced by DeactivateICRU73Scaling(), and G4IonParametrisedLossModel().
◆ ComputeCrossSectionPerAtom()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 370 of file G4IonParametrisedLossModel.cc.
376 {
377
378 // ############## Cross section per atom ################################
379 // Function computes ionization cross section per atom
380 //
381 // See Geant4 physics reference manual (version 9.1), section 9.1.3
382 //
383 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
384 // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
385 //
386 // (Implementation adapted from G4BraggIonModel)
387
388 G4double crosssection = 0.0;
390 G4double maxEnergy = std::min(tmax, maxKinEnergy);
391
392 if(cutEnergy < tmax) {
393
394 G4double energy = kineticEnergy + cacheMass;
395 G4double betaSquared = kineticEnergy *
396 (energy + cacheMass) / (energy * energy);
397
398 crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
399 betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
400
401 crosssection *= twopi_mc2_rcl2 * cacheChargeSquare / betaSquared;
402 }
403
404#ifdef PRINT_DEBUG_CS
406 << G4endl
407 << "# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
408 << G4endl
409 << "# particle =" << particle -> GetParticleName()
410 << G4endl
411 << "# cut(MeV) = " << cutEnergy/MeV
412 << G4endl;
413
415 << std::setw(13) << std::right << "E(MeV)"
416 << std::setw(14) << "CS(um)"
417 << std::setw(14) << "E_max_sec(MeV)"
418 << G4endl
419 << "# ------------------------------------------------------"
420 << G4endl;
421
422 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
423 << std::setw(14) << crosssection / (um * um)
424 << std::setw(14) << tmax / MeV
425 << G4endl;
426#endif
427
428 crosssection *= atomicNumber;
429
430 return crosssection;
431}
G4DLLIMPORT std::ostream G4cout
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double)
Definition: G4IonParametrisedLossModel.cc:200
Referenced by CrossSectionPerVolume().
◆ ComputeDEDXPerVolume()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 454 of file G4IonParametrisedLossModel.cc.
458 {
459
460 // ############## dE/dx ##################################################
461 // Function computes dE/dx values, where following rules are adopted:
462 // A. If the ion-material pair is covered by any native ion data
463 // parameterisation, then:
464 // * This parameterization is used for energies below a given energy
465 // limit,
466 // * whereas above the limit the Bethe-Bloch model is applied, in
467 // combination with an effective charge estimate and high order
468 // correction terms.
469 // A smoothing procedure is applied to dE/dx values computed with
470 // the second approach. The smoothing factor is based on the dE/dx
471 // values of both approaches at the transition energy (high order
472 // correction terms are included in the calculation of the transition
473 // factor).
474 // B. If the particle is a generic ion, the BraggIon and Bethe-Bloch
475 // models are used and a smoothing procedure is applied to values
476 // obtained with the second approach.
477 // C. If the ion-material is not covered by any ion data parameterization
478 // then:
479 // * The BraggIon model is used for energies below a given energy
480 // limit,
481 // * whereas above the limit the Bethe-Bloch model is applied, in
482 // combination with an effective charge estimate and high order
483 // correction terms.
484 // Also in this case, a smoothing procedure is applied to dE/dx values
485 // computed with the second model.
486
487 G4double dEdx = 0.0;
488 UpdateDEDXCache(particle, material, cutEnergy);
489
490 LossTableList::iterator iter = dedxCacheIter;
491
492 if(iter != lossTableList.end()) {
493
494 G4double transitionEnergy = dedxCacheTransitionEnergy;
495
496 if(transitionEnergy > kineticEnergy) {
497
498 dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
499
501 particle,
502 kineticEnergy,
503 cutEnergy);
504 dEdx -= dEdxDeltaRays;
505 }
506 else {
507 G4double massRatio = dedxCacheGenIonMassRatio;
508
509 G4double chargeSquare =
510 GetChargeSquareRatio(particle, material, kineticEnergy);
511
512 G4double scaledKineticEnergy = kineticEnergy * massRatio;
513 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
514
516
517 if(scaledTransitionEnergy >= lowEnergyLimit) {
518
519 dEdx = betheBlochModel -> ComputeDEDXPerVolume(
520 material, genericIon,
521 scaledKineticEnergy, cutEnergy);
522
523 dEdx *= chargeSquare;
524
525 dEdx += corrections -> ComputeIonCorrections(particle,
526 material, kineticEnergy);
527
528 G4double factor = 1.0 + dedxCacheTransitionFactor /
529 kineticEnergy;
530
531 dEdx *= factor;
532 }
533 }
534 }
535 else {
536 G4double massRatio = 1.0;
537 G4double chargeSquare = 1.0;
538
539 if(particle != genericIon) {
540
541 chargeSquare = GetChargeSquareRatio(particle, material, kineticEnergy);
542 massRatio = genericIonPDGMass / particle -> GetPDGMass();
543 }
544
545 G4double scaledKineticEnergy = kineticEnergy * massRatio;
546
548 if(scaledKineticEnergy < lowEnergyLimit) {
549 dEdx = braggIonModel -> ComputeDEDXPerVolume(
550 material, genericIon,
551 scaledKineticEnergy, cutEnergy);
552
553 dEdx *= chargeSquare;
554 }
555 else {
557 material, genericIon,
558 lowEnergyLimit, cutEnergy);
559
561 material, genericIon,
562 lowEnergyLimit, cutEnergy);
563
564 if(particle != genericIon) {
565 G4double chargeSquareLowEnergyLimit =
566 GetChargeSquareRatio(particle, material,
567 lowEnergyLimit / massRatio);
568
569 dEdxLimitParam *= chargeSquareLowEnergyLimit;
570 dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
571
572 dEdxLimitBetheBloch +=
573 corrections -> ComputeIonCorrections(particle,
574 material, lowEnergyLimit / massRatio);
575 }
576
577 G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
578 * lowEnergyLimit / scaledKineticEnergy);
579
580 dEdx = betheBlochModel -> ComputeDEDXPerVolume(
581 material, genericIon,
582 scaledKineticEnergy, cutEnergy);
583
584 dEdx *= chargeSquare;
585
586 if(particle != genericIon) {
587 dEdx += corrections -> ComputeIonCorrections(particle,
588 material, kineticEnergy);
589 }
590
591 dEdx *= factor;
592 }
593
594 }
595
596 if (dEdx < 0.0) dEdx = 0.0;
597
598 return dEdx;
599}
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double)
Definition: G4IonParametrisedLossModel.cc:228
G4double DeltaRayMeanEnergyTransferRate(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
Definition: G4IonParametrisedLossModel.cc:454
Referenced by ComputeDEDXPerVolume(), CorrectionsAlongStep(), and PrintDEDXTable().
◆ ComputeLossForStep()
| G4double G4IonParametrisedLossModel::ComputeLossForStep | ( | const G4MaterialCutsCouple * | matCutsCouple, |
| const G4ParticleDefinition * | particle, | ||
| G4double | kineticEnergy, | ||
| G4double | stepLength | ||
| ) |
Definition at line 1186 of file G4IonParametrisedLossModel.cc.
1190 {
1191
1192 G4double loss = 0.0;
1193
1194 UpdateRangeCache(particle, matCutsCouple);
1195
1196 G4PhysicsVector* energyRange = rangeCacheEnergyRange;
1197 G4PhysicsVector* rangeEnergy = rangeCacheRangeEnergy;
1198
1199 if(energyRange != 0 && rangeEnergy != 0) {
1200 G4bool b;
1201
1202 G4double lowerEnEdge = energyRange -> GetLowEdgeEnergy( 0 );
1203 G4double lowerRangeEdge = rangeEnergy -> GetLowEdgeEnergy( 0 );
1204
1205 // Computing range for pre-step kinetic energy:
1206 G4double range = energyRange -> GetValue(kineticEnergy, b);
1207
1208 // Energy below vector boundary:
1209 if(kineticEnergy < lowerEnEdge) {
1210
1211 range = energyRange -> GetValue(lowerEnEdge, b);
1212 range *= std::sqrt(kineticEnergy / lowerEnEdge);
1213 }
1214
1215#ifdef PRINT_DEBUG
1217 << range / mm << " mm, step = " << stepLength / mm << " mm"
1218 << G4endl;
1219#endif
1220
1221 // Remaining range:
1222 G4double remRange = range - stepLength;
1223
1224 // If range is smaller than step length, the loss is set to kinetic
1225 // energy
1226 if(remRange < 0.0) loss = kineticEnergy;
1227 else if(remRange < lowerRangeEdge) {
1228
1229 G4double ratio = remRange / lowerRangeEdge;
1230 loss = kineticEnergy - ratio * ratio * lowerEnEdge;
1231 }
1232 else {
1233
1234 G4double energy = rangeEnergy -> GetValue(range - stepLength, b);
1235 loss = kineticEnergy - energy;
1236 }
1237 }
1238
1239 if(loss < 0.0) loss = 0.0;
1240
1241 return loss;
1242}
Definition: G4PhysicsVector.hh:77
Referenced by CorrectionsAlongStep().
◆ CorrectionsAlongStep()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 919 of file G4IonParametrisedLossModel.cc.
924 {
925
926 // ############## Corrections for along step energy loss calculation ######
927 // The computed energy loss (due to electronic stopping) is overwritten
928 // by this function if an ion data parameterization is available for the
929 // current ion-material pair.
930 // No action on the energy loss (due to electronic stopping) is performed
931 // if no parameterization is available. In this case the original
932 // generic ion tables (in combination with the effective charge) are used
933 // in the along step DoIt function.
934 //
935 //
936 // (Implementation partly adapted from G4BraggIonModel/G4BetheBlochModel)
937
940
941 G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
942
943 if(kineticEnergy == eloss) { return; }
944
946 size_t cutIndex = couple -> GetIndex();
947 cutEnergy = cutEnergies[cutIndex];
948
949 UpdateDEDXCache(particle, material, cutEnergy);
950
951 LossTableList::iterator iter = dedxCacheIter;
952
953 // If parameterization for ions is available the electronic energy loss
954 // is overwritten
955 if(iter != lossTableList.end()) {
956
957 // The energy loss is calculated using the ComputeDEDXPerVolume function
958 // and the step length (it is assumed that dE/dx does not change
959 // considerably along the step)
960 eloss =
961 length * ComputeDEDXPerVolume(material, particle,
962 kineticEnergy, cutEnergy);
963
964#ifdef PRINT_DEBUG
965 G4cout.precision(6);
967 << G4endl
968 << "# G4IonParametrisedLossModel::CorrectionsAlongStep"
969 << G4endl
970 << "# cut(MeV) = " << cutEnergy/MeV
971 << G4endl;
972
974 << std::setw(13) << std::right << "E(MeV)"
975 << std::setw(14) << "l(um)"
976 << std::setw(14) << "l*dE/dx(MeV)"
977 << std::setw(14) << "(l*dE/dx)/E"
978 << G4endl
979 << "# ------------------------------------------------------"
980 << G4endl;
981
982 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
983 << std::setw(14) << length / um
984 << std::setw(14) << eloss / MeV
985 << std::setw(14) << eloss / kineticEnergy * 100.0
986 << G4endl;
987#endif
988
989 // If the energy loss exceeds a certain fraction of the kinetic energy
990 // (the fraction is indicated by the parameter "energyLossLimit") then
991 // the range tables are used to derive a more accurate value of the
992 // energy loss
993 if(eloss > energyLossLimit * kineticEnergy) {
994
995 eloss = ComputeLossForStep(couple, particle,
996 kineticEnergy,length);
997
998#ifdef PRINT_DEBUG
1000 << G4endl;
1001
1002 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
1003 << std::setw(14) << length / um
1004 << std::setw(14) << eloss / MeV
1005 << std::setw(14) << eloss / kineticEnergy * 100.0
1006 << G4endl;
1007#endif
1008
1009 }
1010 }
1011
1012 // For all corrections below a kinetic energy between the Pre- and
1013 // Post-step energy values is used
1014 G4double energy = kineticEnergy - eloss * 0.5;
1015 if(energy < 0.0) energy = kineticEnergy * 0.5;
1016
1017 G4double chargeSquareRatio = corrections ->
1018 EffectiveChargeSquareRatio(particle,
1019 material,
1020 energy);
1021 GetModelOfFluctuations() -> SetParticleAndCharge(particle,
1022 chargeSquareRatio);
1023
1024 // A correction is applied considering the change of the effective charge
1025 // along the step (the parameter "corrFactor" refers to the effective
1026 // charge at the beginning of the step). Note: the correction is not
1027 // applied for energy loss values deriving directly from parameterized
1028 // ion stopping power tables
1029 G4double transitionEnergy = dedxCacheTransitionEnergy;
1030
1031 if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
1032 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1033 material,
1034 energy);
1035
1036 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1037 eloss *= chargeSquareRatioCorr;
1038 }
1039 else if (iter == lossTableList.end()) {
1040
1041 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1042 material,
1043 energy);
1044
1045 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1046 eloss *= chargeSquareRatioCorr;
1047 }
1048
1049 // Ion high order corrections are applied if the current model does not
1050 // overwrite the energy loss (i.e. when the effective charge approach is
1051 // used)
1052 if(iter == lossTableList.end()) {
1053
1054 G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
1056
1057 // Corrections are only applied in the Bethe-Bloch energy region
1058 if(scaledKineticEnergy > lowEnergyLimit)
1059 eloss += length *
1060 corrections -> IonHighOrderCorrections(particle, couple, energy);
1061 }
1062}
G4double ComputeLossForStep(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double, G4double)
Definition: G4IonParametrisedLossModel.cc:1186
Definition: G4Material.hh:119
Definition: G4ParticleDefinition.hh:68
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:501
◆ CrossSectionPerVolume()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 435 of file G4IonParametrisedLossModel.cc.
440 {
441
442 G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
444 kineticEnergy,
445 nbElecPerVolume, 0,
446 cutEnergy,
447 maxEnergy);
448
449 return cross;
450}
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
Definition: G4IonParametrisedLossModel.cc:370
◆ DeactivateICRU73Scaling()
| void G4IonParametrisedLossModel::DeactivateICRU73Scaling | ( | ) |
Definition at line 1328 of file G4IonParametrisedLossModel.cc.
1328 {
1329
1332}
G4bool RemoveDEDXTable(const G4String &name)
Definition: G4IonParametrisedLossModel.cc:1289
◆ DeltaRayMeanEnergyTransferRate()
|
inline |
Referenced by ComputeDEDXPerVolume().
◆ GetChargeSquareRatio()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 228 of file G4IonParametrisedLossModel.cc.
231 { // Kinetic energy
232
233 G4double chargeSquareRatio = corrections ->
234 EffectiveChargeSquareRatio(particle,
235 material,
236 kineticEnergy);
237 corrFactor = chargeSquareRatio *
238 corrections -> EffectiveChargeCorrection(particle,
239 material,
240 kineticEnergy);
241 return corrFactor;
242}
Referenced by ComputeDEDXPerVolume().
◆ GetParticleCharge()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 246 of file G4IonParametrisedLossModel.cc.
249 { // Kinetic energy
250
252}
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double)
Definition: G4IonParametrisedLossModel.cc:246
Referenced by GetParticleCharge().
◆ Initialise()
|
virtual |
Implements G4VEmModel.
Definition at line 256 of file G4IonParametrisedLossModel.cc.
258 {
259
260 // Cached parameters are reset
261 cacheParticle = 0;
262 cacheMass = 0;
263 cacheElecMassRatio = 0;
264 cacheChargeSquare = 0;
265
266 // Cached parameters are reset
267 rangeCacheParticle = 0;
268 rangeCacheMatCutsCouple = 0;
269 rangeCacheEnergyRange = 0;
270 rangeCacheRangeEnergy = 0;
271
272 // Cached parameters are reset
273 dedxCacheParticle = 0;
274 dedxCacheMaterial = 0;
275 dedxCacheEnergyCut = 0;
276 dedxCacheIter = lossTableList.end();
277 dedxCacheTransitionEnergy = 0.0;
278 dedxCacheTransitionFactor = 0.0;
279 dedxCacheGenIonMassRatio = 0.0;
280
281 // The cache of loss tables is cleared
282 LossTableList::iterator iterTables = lossTableList.begin();
283 LossTableList::iterator iterTables_end = lossTableList.end();
284
285 for(;iterTables != iterTables_end; iterTables++)
286 (*iterTables) -> ClearCache();
287
288 // Range vs energy and energy vs range vectors from previous runs are
289 // cleared
290 RangeEnergyTable::iterator iterRange = r.begin();
291 RangeEnergyTable::iterator iterRange_end = r.end();
292
293 for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second;
294 r.clear();
295
296 EnergyRangeTable::iterator iterEnergy = E.begin();
297 EnergyRangeTable::iterator iterEnergy_end = E.end();
298
299 for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second;
300 E.clear();
301
302 // The cut energies are (re)loaded
303 size_t size = cuts.size();
304 cutEnergies.clear();
305 for(size_t i = 0; i < size; i++) cutEnergies.push_back(cuts[i]);
306
307 // All dE/dx vectors are built
310 size_t nmbCouples = coupleTable -> GetTableSize();
311
312#ifdef PRINT_TABLE_BUILT
314 << " Building dE/dx vectors:"
315 << G4endl;
316#endif
317
318 for (size_t i = 0; i < nmbCouples; i++) {
319
321 coupleTable -> GetMaterialCutsCouple(i);
322
324 // G4ProductionCuts* productionCuts = couple -> GetProductionCuts();
325
327
328 LossTableList::iterator iter = lossTableList.begin();
329 LossTableList::iterator iter_end = lossTableList.end();
330
331 for(;iter != iter_end; iter++) {
332
333 if(*iter == 0) {
335 << " Skipping illegal table."
336 << G4endl;
337 }
338
339 G4bool isApplicable =
340 (*iter) -> BuildDEDXTable(atomicNumberIon, material);
341 if(isApplicable) {
342
343#ifdef PRINT_TABLE_BUILT
347 << G4endl;
348#endif
349 break;
350 }
351 }
352 }
353 }
354
355 // The particle change object
356 if(! particleChangeLoss) {
357 particleChangeLoss = GetParticleChangeForLoss();
358 braggIonModel -> SetParticleChange(particleChangeLoss, 0);
359 betheBlochModel -> SetParticleChange(particleChangeLoss, 0);
360 }
361
362 // The G4BraggIonModel and G4BetheBlochModel instances are initialised with
363 // the same settings as the current model:
364 braggIonModel -> Initialise(particle, cuts);
365 betheBlochModel -> Initialise(particle, cuts);
366}
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4IonParametrisedLossModel.cc:256
Definition: G4MaterialCutsCouple.hh:51
Definition: G4ProductionCutsTable.hh:72
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=0)
Definition: G4VEmModel.cc:318
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:95
Referenced by Initialise().
◆ IsApplicable()
|
inline |
Referenced by PrintDEDXTableHandlers().
◆ MaxSecondaryEnergy()
|
protectedvirtual |
Reimplemented from G4VEmModel.
Definition at line 200 of file G4IonParametrisedLossModel.cc.
202 {
203
204 // ############## Maximum energy of secondaries ##########################
205 // Function computes maximum energy of secondary electrons which are
206 // released by an ion
207 //
208 // See Geant4 physics reference manual (version 9.1), section 9.1.1
209 //
210 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
211 // C.Caso et al. (Part. Data Group), Europ. Phys. Jour. C 3 1 (1998).
212 // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
213 //
214 // (Implementation adapted from G4BraggIonModel)
215
216 if(particle != cacheParticle) UpdateCache(particle);
217
218 G4double tau = kineticEnergy/cacheMass;
219 G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
220 (1. + 2.0 * (tau + 1.) * cacheElecMassRatio +
221 cacheElecMassRatio * cacheElecMassRatio);
222
223 return tmax;
224}
Referenced by ComputeCrossSectionPerAtom().
◆ MinEnergyCut()
|
virtual |
Definition at line 190 of file G4IonParametrisedLossModel.cc.
192 {
193
194 return couple -> GetMaterial() -> GetIonisation() ->
195 GetMeanExcitationEnergy();
196}
◆ PrintDEDXTable()
| void G4IonParametrisedLossModel::PrintDEDXTable | ( | const G4ParticleDefinition * | particle, |
| const G4Material * | material, | ||
| G4double | lowerBoundary, | ||
| G4double | upperBoundary, | ||
| G4int | numBins, | ||
| G4bool | logScaleEnergy | ||
| ) |
Definition at line 603 of file G4IonParametrisedLossModel.cc.
609 { // Logarithmic scaling of energy
610
611 G4double atomicMassNumber = particle -> GetAtomicMass();
612 G4double materialDensity = material -> GetDensity();
613
616 << " of density " << materialDensity / g * cm3
617 << " g/cm3"
618 << G4endl
619 << "# Projectile mass number A1 = " << atomicMassNumber
620 << G4endl
621 << "# ------------------------------------------------------"
622 << G4endl;
624 << std::setw(13) << std::right << "E"
625 << std::setw(14) << "E/A1"
626 << std::setw(14) << "dE/dx"
627 << std::setw(14) << "1/rho*dE/dx"
628 << G4endl;
630 << std::setw(13) << std::right << "(MeV)"
631 << std::setw(14) << "(MeV)"
632 << std::setw(14) << "(MeV/cm)"
633 << std::setw(14) << "(MeV*cm2/mg)"
634 << G4endl
635 << "# ------------------------------------------------------"
636 << G4endl;
637
638 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
639 G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
640
641 if(logScaleEnergy) {
642
643 energyLowerBoundary = std::log(energyLowerBoundary);
644 energyUpperBoundary = std::log(energyUpperBoundary);
645 }
646
647 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
648 G4double(nmbBins);
649
650 for(int i = 0; i < numBins + 1; i++) {
651
652 G4double energy = energyLowerBoundary + i * deltaEnergy;
653 if(logScaleEnergy) energy = std::exp(energy);
654
656 G4cout.precision(6);
657 G4cout << std::setw(14) << std::right << energy / MeV
658 << std::setw(14) << energy / atomicMassNumber / MeV
659 << std::setw(14) << dedx / MeV * cm
660 << std::setw(14) << dedx / materialDensity / (MeV*cm2/(0.001*g))
661 << G4endl;
662 }
663}
Referenced by PrintDEDXTableHandlers().
◆ PrintDEDXTableHandlers()
| void G4IonParametrisedLossModel::PrintDEDXTableHandlers | ( | const G4ParticleDefinition * | particle, |
| const G4Material * | material, | ||
| G4double | lowerBoundary, | ||
| G4double | upperBoundary, | ||
| G4int | numBins, | ||
| G4bool | logScaleEnergy | ||
| ) |
Definition at line 667 of file G4IonParametrisedLossModel.cc.
673 { // Logarithmic scaling of energy
674
675 LossTableList::iterator iter = lossTableList.begin();
676 LossTableList::iterator iter_end = lossTableList.end();
677
678 for(;iter != iter_end; iter++) {
680 if(isApplicable) {
681 (*iter) -> PrintDEDXTable(particle, material,
682 lowerBoundary, upperBoundary,
683 numBins,logScaleEnergy);
684 break;
685 }
686 }
687}
LossTableList::iterator IsApplicable(const G4ParticleDefinition *, const G4Material *)
void PrintDEDXTable(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
Definition: G4IonParametrisedLossModel.cc:603
◆ RemoveDEDXTable()
Definition at line 1289 of file G4IonParametrisedLossModel.cc.
1290 {
1291
1292 LossTableList::iterator iter = lossTableList.begin();
1293 LossTableList::iterator iter_end = lossTableList.end();
1294
1295 for(;iter != iter_end; iter++) {
1297
1298 if(tableName == nam) {
1299 delete (*iter);
1300
1301 // Remove from table list
1302 lossTableList.erase(iter);
1303
1304 // Range vs energy and energy vs range vectors are cleared
1305 RangeEnergyTable::iterator iterRange = r.begin();
1306 RangeEnergyTable::iterator iterRange_end = r.end();
1307
1308 for(;iterRange != iterRange_end; iterRange++)
1309 delete iterRange -> second;
1310 r.clear();
1311
1312 EnergyRangeTable::iterator iterEnergy = E.begin();
1313 EnergyRangeTable::iterator iterEnergy_end = E.end();
1314
1315 for(;iterEnergy != iterEnergy_end; iterEnergy++)
1316 delete iterEnergy -> second;
1317 E.clear();
1318
1319 return true;
1320 }
1321 }
1322
1323 return false;
1324}
Referenced by DeactivateICRU73Scaling().
◆ SampleSecondaries()
|
virtual |
Implements G4VEmModel.
Definition at line 691 of file G4IonParametrisedLossModel.cc.
696 {
697
698
699 // ############## Sampling of secondaries #################################
700 // The probability density function (pdf) of the kinetic energy T of a
701 // secondary electron may be written as:
702 // pdf(T) = f(T) * g(T)
703 // where
704 // f(T) = (Tmax - Tcut) / (Tmax * Tcut) * (1 / T^2)
705 // g(T) = 1 - beta^2 * T / Tmax
706 // where Tmax is the maximum kinetic energy of the secondary, Tcut
707 // is the lower energy cut and beta is the kinetic energy of the
708 // projectile.
709 //
710 // Sampling of the kinetic energy of a secondary electron:
711 // 1) T0 is sampled from f(T) using the cumulated distribution function
712 // F(T) = int_Tcut^T f(T')dT'
713 // 2) T is accepted or rejected by evaluating the rejection function g(T)
714 // at the sampled energy T0 against a randomly sampled value
715 //
716 //
717 // See Geant4 physics reference manual (version 9.1), section 9.1.4
718 //
719 //
720 // Reference pdf: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
721 //
722 // (Implementation adapted from G4BraggIonModel)
723
725 G4double maxKinEnergySec =
726 std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
727
728 if(cutKinEnergySec >= maxKinEnergySec) return;
729
730 G4double kineticEnergy = particle -> GetKineticEnergy();
732
733 G4double energy = kineticEnergy + cacheMass;
734 G4double betaSquared = kineticEnergy *
735 (energy + cacheMass) / (energy * energy);
736
737 G4double kinEnergySec;
738 G4double grej;
739
740 do {
741
742 // Sampling kinetic energy from f(T) (using F(T)):
744 kinEnergySec = cutKinEnergySec * maxKinEnergySec /
745 (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
746
747 // Deriving the value of the rejection function at the obtained kinetic
748 // energy:
749 grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
750
751 if(grej > 1.0) {
753 << "Majorant 1.0 < "
754 << grej << " for e= " << kinEnergySec
755 << G4endl;
756 }
757
759
760 G4double momentumSec =
761 std::sqrt(kinEnergySec * (kinEnergySec + 2.0 * electron_mass_c2));
762
763 G4double totMomentum = energy*std::sqrt(betaSquared);
764 G4double cost = kinEnergySec * (energy + electron_mass_c2) /
765 (momentumSec * totMomentum);
766 if(cost > 1.0) cost = 1.0;
767 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
768
770
771 G4ThreeVector directionSec(sint*std::cos(phi),sint*std::sin(phi), cost) ;
772 directionSec.rotateUz(direction);
773
774 // create G4DynamicParticle object for delta ray
776 directionSec,
777 kinEnergySec);
778
779 secondaries -> push_back(delta);
780
781 // Change kinematics of primary particle
782 kineticEnergy -= kinEnergySec;
783 G4ThreeVector finalP = direction*totMomentum - directionSec*momentumSec;
784 finalP = finalP.unit();
785
786 particleChangeLoss -> SetProposedKineticEnergy(kineticEnergy);
787 particleChangeLoss -> SetProposedMomentumDirection(finalP);
788}
Definition: ThreeVector.h:41
Hep3Vector unit() const
Definition: G4DynamicParticle.hh:74
const G4ThreeVector & GetMomentumDirection() const
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:399
◆ SetEnergyLossLimit()
|
inline |
The documentation for this class was generated from the following files:
