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 | |
| G4ElementData * | fElementData |
| G4VParticleChange * | pParticleChange |
| G4PhysicsTable * | xSectionTable |
| const G4Material * | pBaseMaterial |
| const std::vector< G4double > * | theDensityFactor |
| const std::vector< G4int > * | theDensityIdx |
| size_t | idxTable |
| G4bool | lossFlucFlag |
| G4double | inveplus |
| G4double | pFactor |
Detailed Description
Definition at line 92 of file G4IonParametrisedLossModel.hh.
Constructor & Destructor Documentation
◆ G4IonParametrisedLossModel()
| G4IonParametrisedLossModel::G4IonParametrisedLossModel | ( | const G4ParticleDefinition * | particle = 0, |
| const G4String & | name = "ParamICRU73" |
||
| ) |
Definition at line 105 of file G4IonParametrisedLossModel.cc.
108 : G4VEmModel(nam),
109 braggIonModel(0),
110 betheBlochModel(0),
111 nmbBins(90),
112 nmbSubBins(100),
113 particleChangeLoss(0),
114 corrFactor(1.0),
115 energyLossLimit(0.01),
116 cutEnergies(0),
117 isInitialised(false)
118{
119 genericIon = G4GenericIon::Definition();
120 genericIonPDGMass = genericIon->GetPDGMass();
122
123 // The Bragg ion and Bethe Bloch models are instantiated
126
127 // The boundaries for the range tables are set
128 lowerEnergyEdgeIntegr = 0.025 * MeV;
129 upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit();
130
131 // Cache parameters are set
132 cacheParticle = 0;
133 cacheMass = 0;
134 cacheElecMassRatio = 0;
135 cacheChargeSquare = 0;
136
137 // Cache parameters are set
138 rangeCacheParticle = 0;
139 rangeCacheMatCutsCouple = 0;
140 rangeCacheEnergyRange = 0;
141 rangeCacheRangeEnergy = 0;
142
143 // Cache parameters are set
144 dedxCacheParticle = 0;
145 dedxCacheMaterial = 0;
146 dedxCacheEnergyCut = 0;
147 dedxCacheIter = lossTableList.end();
148 dedxCacheTransitionEnergy = 0.0;
149 dedxCacheTransitionFactor = 0.0;
150 dedxCacheGenIonMassRatio = 0.0;
151
152 // default generator
154}
Definition: G4BetheBlochModel.hh:63
Definition: G4BraggIonModel.hh:69
Definition: G4DeltaAngle.hh:56
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:87
G4EmCorrections * EmCorrections()
Definition: G4LossTableManager.hh:397
Definition: G4VEmModel.hh:103
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:618
◆ ~G4IonParametrisedLossModel()
|
virtual |
Definition at line 158 of file G4IonParametrisedLossModel.cc.
158 {
159
160 // dE/dx table objects are deleted and the container is cleared
161 LossTableList::iterator iterTables = lossTableList.begin();
162 LossTableList::iterator iterTables_end = lossTableList.end();
163
164 for(;iterTables != iterTables_end; ++iterTables) { delete *iterTables; }
165 lossTableList.clear();
166
167 // range table
168 RangeEnergyTable::iterator itr = r.begin();
169 RangeEnergyTable::iterator itr_end = r.end();
170 for(;itr != itr_end; ++itr) { delete itr->second; }
171 r.clear();
172
173 // inverse range
174 EnergyRangeTable::iterator ite = E.begin();
175 EnergyRangeTable::iterator ite_end = E.end();
176 for(;ite != ite_end; ++ite) { delete ite->second; }
177 E.clear();
178
179}
Member Function Documentation
◆ AddDEDXTable()
| G4bool G4IonParametrisedLossModel::AddDEDXTable | ( | const G4String & | name, |
| G4VIonDEDXTable * | table, | ||
| G4VIonDEDXScalingAlgorithm * | algorithm = 0 |
||
| ) |
Definition at line 1237 of file G4IonParametrisedLossModel.cc.
1240 {
1241
1242 if(table == 0) {
1244 << " add table: Invalid pointer."
1245 << G4endl;
1246
1247 return false;
1248 }
1249
1250 // Checking uniqueness of name
1251 LossTableList::iterator iter = lossTableList.begin();
1252 LossTableList::iterator iter_end = lossTableList.end();
1253
1254 for(;iter != iter_end; ++iter) {
1256
1257 if(tableName == nam) {
1259 << " add table: Name already exists."
1260 << G4endl;
1261
1262 return false;
1263 }
1264 }
1265
1266 G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm;
1267 if(scalingAlgorithm == 0)
1269
1270 G4IonDEDXHandler* handler =
1272
1273 lossTableList.push_front(handler);
1274
1275 return true;
1276}
G4GLOB_DLL std::ostream G4cout
Definition: G4IonDEDXHandler.hh:82
Definition: G4String.hh:53
Referenced by Initialise().
◆ ComputeCrossSectionPerAtom()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 368 of file G4IonParametrisedLossModel.cc.
374 {
375
376 // ############## Cross section per atom ################################
377 // Function computes ionization cross section per atom
378 //
379 // See Geant4 physics reference manual (version 9.1), section 9.1.3
380 //
381 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
382 // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
383 //
384 // (Implementation adapted from G4BraggIonModel)
385
386 G4double crosssection = 0.0;
388 G4double maxEnergy = std::min(tmax, maxKinEnergy);
389
390 if(cutEnergy < tmax) {
391
393 G4double betaSquared = kineticEnergy *
394 (energy + cacheMass) / (energy * energy);
395
396 crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
397 betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
398
399 crosssection *= twopi_mc2_rcl2 * cacheChargeSquare / betaSquared;
400 }
401
402#ifdef PRINT_DEBUG_CS
404 << G4endl
405 << "# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
406 << G4endl
407 << "# particle =" << particle -> GetParticleName()
408 << G4endl
409 << "# cut(MeV) = " << cutEnergy/MeV
410 << G4endl;
411
413 << std::setw(13) << std::right << "E(MeV)"
414 << std::setw(14) << "CS(um)"
415 << std::setw(14) << "E_max_sec(MeV)"
416 << G4endl
417 << "# ------------------------------------------------------"
418 << G4endl;
419
420 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
421 << std::setw(14) << crosssection / (um * um)
422 << std::setw(14) << tmax / MeV
423 << G4endl;
424#endif
425
426 crosssection *= atomicNumber;
427
428 return crosssection;
429}
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double)
Definition: G4IonParametrisedLossModel.cc:192
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:158
Referenced by CrossSectionPerVolume().
◆ ComputeDEDXPerVolume()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 452 of file G4IonParametrisedLossModel.cc.
456 {
457
458 // ############## dE/dx ##################################################
459 // Function computes dE/dx values, where following rules are adopted:
460 // A. If the ion-material pair is covered by any native ion data
461 // parameterisation, then:
462 // * This parameterization is used for energies below a given energy
463 // limit,
464 // * whereas above the limit the Bethe-Bloch model is applied, in
465 // combination with an effective charge estimate and high order
466 // correction terms.
467 // A smoothing procedure is applied to dE/dx values computed with
468 // the second approach. The smoothing factor is based on the dE/dx
469 // values of both approaches at the transition energy (high order
470 // correction terms are included in the calculation of the transition
471 // factor).
472 // B. If the particle is a generic ion, the BraggIon and Bethe-Bloch
473 // models are used and a smoothing procedure is applied to values
474 // obtained with the second approach.
475 // C. If the ion-material is not covered by any ion data parameterization
476 // then:
477 // * The BraggIon model is used for energies below a given energy
478 // limit,
479 // * whereas above the limit the Bethe-Bloch model is applied, in
480 // combination with an effective charge estimate and high order
481 // correction terms.
482 // Also in this case, a smoothing procedure is applied to dE/dx values
483 // computed with the second model.
484
485 G4double dEdx = 0.0;
486 UpdateDEDXCache(particle, material, cutEnergy);
487
488 LossTableList::iterator iter = dedxCacheIter;
489
490 if(iter != lossTableList.end()) {
491
492 G4double transitionEnergy = dedxCacheTransitionEnergy;
493
494 if(transitionEnergy > kineticEnergy) {
495
496 dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
497
499 particle,
500 kineticEnergy,
501 cutEnergy);
502 dEdx -= dEdxDeltaRays;
503 }
504 else {
505 G4double massRatio = dedxCacheGenIonMassRatio;
506
507 G4double chargeSquare =
508 GetChargeSquareRatio(particle, material, kineticEnergy);
509
510 G4double scaledKineticEnergy = kineticEnergy * massRatio;
511 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
512
514
515 if(scaledTransitionEnergy >= lowEnergyLimit) {
516
517 dEdx = betheBlochModel -> ComputeDEDXPerVolume(
518 material, genericIon,
519 scaledKineticEnergy, cutEnergy);
520
521 dEdx *= chargeSquare;
522
523 dEdx += corrections -> ComputeIonCorrections(particle,
524 material, kineticEnergy);
525
526 G4double factor = 1.0 + dedxCacheTransitionFactor /
527 kineticEnergy;
528
529 dEdx *= factor;
530 }
531 }
532 }
533 else {
534 G4double massRatio = 1.0;
535 G4double chargeSquare = 1.0;
536
537 if(particle != genericIon) {
538
539 chargeSquare = GetChargeSquareRatio(particle, material, kineticEnergy);
540 massRatio = genericIonPDGMass / particle -> GetPDGMass();
541 }
542
543 G4double scaledKineticEnergy = kineticEnergy * massRatio;
544
546 if(scaledKineticEnergy < lowEnergyLimit) {
547 dEdx = braggIonModel -> ComputeDEDXPerVolume(
548 material, genericIon,
549 scaledKineticEnergy, cutEnergy);
550
551 dEdx *= chargeSquare;
552 }
553 else {
555 material, genericIon,
556 lowEnergyLimit, cutEnergy);
557
559 material, genericIon,
560 lowEnergyLimit, cutEnergy);
561
562 if(particle != genericIon) {
563 G4double chargeSquareLowEnergyLimit =
564 GetChargeSquareRatio(particle, material,
565 lowEnergyLimit / massRatio);
566
567 dEdxLimitParam *= chargeSquareLowEnergyLimit;
568 dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
569
570 dEdxLimitBetheBloch +=
571 corrections -> ComputeIonCorrections(particle,
572 material, lowEnergyLimit / massRatio);
573 }
574
575 G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
576 * lowEnergyLimit / scaledKineticEnergy);
577
578 dEdx = betheBlochModel -> ComputeDEDXPerVolume(
579 material, genericIon,
580 scaledKineticEnergy, cutEnergy);
581
582 dEdx *= chargeSquare;
583
584 if(particle != genericIon) {
585 dEdx += corrections -> ComputeIonCorrections(particle,
586 material, kineticEnergy);
587 }
588
589 dEdx *= factor;
590 }
591
592 }
593
594 if (dEdx < 0.0) dEdx = 0.0;
595
596 return dEdx;
597}
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double)
Definition: G4IonParametrisedLossModel.cc:220
G4double DeltaRayMeanEnergyTransferRate(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
Definition: G4IonParametrisedLossModel.cc:452
Referenced by ComputeDEDXPerVolume(), CorrectionsAlongStep(), and PrintDEDXTable().
◆ ComputeLossForStep()
| G4double G4IonParametrisedLossModel::ComputeLossForStep | ( | const G4MaterialCutsCouple * | matCutsCouple, |
| const G4ParticleDefinition * | particle, | ||
| G4double | kineticEnergy, | ||
| G4double | stepLength | ||
| ) |
Definition at line 1178 of file G4IonParametrisedLossModel.cc.
1182 {
1183
1184 G4double loss = 0.0;
1185
1186 UpdateRangeCache(particle, matCutsCouple);
1187
1188 G4PhysicsVector* energyRange = rangeCacheEnergyRange;
1189 G4PhysicsVector* rangeEnergy = rangeCacheRangeEnergy;
1190
1191 if(energyRange != 0 && rangeEnergy != 0) {
1192
1193 G4double lowerEnEdge = energyRange -> Energy( 0 );
1194 G4double lowerRangeEdge = rangeEnergy -> Energy( 0 );
1195
1196 // Computing range for pre-step kinetic energy:
1198
1199 // Energy below vector boundary:
1200 if(kineticEnergy < lowerEnEdge) {
1201
1202 range = energyRange -> Value(lowerEnEdge);
1203 range *= std::sqrt(kineticEnergy / lowerEnEdge);
1204 }
1205
1206#ifdef PRINT_DEBUG
1208 << range / mm << " mm, step = " << stepLength / mm << " mm"
1209 << G4endl;
1210#endif
1211
1212 // Remaining range:
1213 G4double remRange = range - stepLength;
1214
1215 // If range is smaller than step length, the loss is set to kinetic
1216 // energy
1217 if(remRange < 0.0) loss = kineticEnergy;
1218 else if(remRange < lowerRangeEdge) {
1219
1220 G4double ratio = remRange / lowerRangeEdge;
1221 loss = kineticEnergy - ratio * ratio * lowerEnEdge;
1222 }
1223 else {
1224
1226 loss = kineticEnergy - energy;
1227 }
1228 }
1229
1230 if(loss < 0.0) loss = 0.0;
1231
1232 return loss;
1233}
Definition: G4PhysicsVector.hh:61
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:415
Referenced by CorrectionsAlongStep().
◆ CorrectionsAlongStep()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 912 of file G4IonParametrisedLossModel.cc.
917 {
918
919 // ############## Corrections for along step energy loss calculation ######
920 // The computed energy loss (due to electronic stopping) is overwritten
921 // by this function if an ion data parameterization is available for the
922 // current ion-material pair.
923 // No action on the energy loss (due to electronic stopping) is performed
924 // if no parameterization is available. In this case the original
925 // generic ion tables (in combination with the effective charge) are used
926 // in the along step DoIt function.
927 //
928 //
929 // (Implementation partly adapted from G4BraggIonModel/G4BetheBlochModel)
930
933
934 G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
935
936 if(kineticEnergy == eloss) { return; }
937
939 size_t cutIndex = couple -> GetIndex();
940 cutEnergy = cutEnergies[cutIndex];
941
942 UpdateDEDXCache(particle, material, cutEnergy);
943
944 LossTableList::iterator iter = dedxCacheIter;
945
946 // If parameterization for ions is available the electronic energy loss
947 // is overwritten
948 if(iter != lossTableList.end()) {
949
950 // The energy loss is calculated using the ComputeDEDXPerVolume function
951 // and the step length (it is assumed that dE/dx does not change
952 // considerably along the step)
953 eloss =
954 length * ComputeDEDXPerVolume(material, particle,
955 kineticEnergy, cutEnergy);
956
957#ifdef PRINT_DEBUG
958 G4cout.precision(6);
960 << G4endl
961 << "# G4IonParametrisedLossModel::CorrectionsAlongStep"
962 << G4endl
963 << "# cut(MeV) = " << cutEnergy/MeV
964 << G4endl;
965
967 << std::setw(13) << std::right << "E(MeV)"
968 << std::setw(14) << "l(um)"
969 << std::setw(14) << "l*dE/dx(MeV)"
970 << std::setw(14) << "(l*dE/dx)/E"
971 << G4endl
972 << "# ------------------------------------------------------"
973 << G4endl;
974
975 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
976 << std::setw(14) << length / um
977 << std::setw(14) << eloss / MeV
978 << std::setw(14) << eloss / kineticEnergy * 100.0
979 << G4endl;
980#endif
981
982 // If the energy loss exceeds a certain fraction of the kinetic energy
983 // (the fraction is indicated by the parameter "energyLossLimit") then
984 // the range tables are used to derive a more accurate value of the
985 // energy loss
986 if(eloss > energyLossLimit * kineticEnergy) {
987
988 eloss = ComputeLossForStep(couple, particle,
989 kineticEnergy,length);
990
991#ifdef PRINT_DEBUG
993 << G4endl;
994
995 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
996 << std::setw(14) << length / um
997 << std::setw(14) << eloss / MeV
998 << std::setw(14) << eloss / kineticEnergy * 100.0
999 << G4endl;
1000#endif
1001
1002 }
1003 }
1004
1005 // For all corrections below a kinetic energy between the Pre- and
1006 // Post-step energy values is used
1009
1010 G4double chargeSquareRatio = corrections ->
1011 EffectiveChargeSquareRatio(particle,
1012 material,
1013 energy);
1014 GetModelOfFluctuations() -> SetParticleAndCharge(particle,
1015 chargeSquareRatio);
1016
1017 // A correction is applied considering the change of the effective charge
1018 // along the step (the parameter "corrFactor" refers to the effective
1019 // charge at the beginning of the step). Note: the correction is not
1020 // applied for energy loss values deriving directly from parameterized
1021 // ion stopping power tables
1022 G4double transitionEnergy = dedxCacheTransitionEnergy;
1023
1024 if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
1025 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1026 material,
1027 energy);
1028
1029 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1030 eloss *= chargeSquareRatioCorr;
1031 }
1032 else if (iter == lossTableList.end()) {
1033
1034 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1035 material,
1036 energy);
1037
1038 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1039 eloss *= chargeSquareRatioCorr;
1040 }
1041
1042 // Ion high order corrections are applied if the current model does not
1043 // overwrite the energy loss (i.e. when the effective charge approach is
1044 // used)
1045 if(iter == lossTableList.end()) {
1046
1047 G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
1049
1050 // Corrections are only applied in the Bethe-Bloch energy region
1051 if(scaledKineticEnergy > lowEnergyLimit)
1052 eloss += length *
1053 corrections -> IonHighOrderCorrections(particle, couple, energy);
1054 }
1055}
G4double ComputeLossForStep(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double, G4double)
Definition: G4IonParametrisedLossModel.cc:1178
Definition: G4Material.hh:118
Definition: G4ParticleDefinition.hh:60
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:604
◆ CrossSectionPerVolume()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 433 of file G4IonParametrisedLossModel.cc.
438 {
439
440 G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
442 kineticEnergy,
443 nbElecPerVolume, 0,
444 cutEnergy,
445 maxEnergy);
446
447 return cross;
448}
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
Definition: G4IonParametrisedLossModel.cc:368
◆ DeltaRayMeanEnergyTransferRate()
|
inline |
Referenced by ComputeDEDXPerVolume().
◆ GetChargeSquareRatio()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 220 of file G4IonParametrisedLossModel.cc.
223 { // Kinetic energy
224
225 G4double chargeSquareRatio = corrections ->
226 EffectiveChargeSquareRatio(particle,
227 material,
228 kineticEnergy);
229 corrFactor = chargeSquareRatio *
230 corrections -> EffectiveChargeCorrection(particle,
231 material,
232 kineticEnergy);
233 return corrFactor;
234}
Referenced by ComputeDEDXPerVolume().
◆ GetParticleCharge()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 238 of file G4IonParametrisedLossModel.cc.
241 { // Kinetic energy
242
244}
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double)
Definition: G4IonParametrisedLossModel.cc:238
Referenced by GetParticleCharge().
◆ Initialise()
|
virtual |
Implements G4VEmModel.
Definition at line 248 of file G4IonParametrisedLossModel.cc.
250 {
251
252 // Cached parameters are reset
253 cacheParticle = 0;
254 cacheMass = 0;
255 cacheElecMassRatio = 0;
256 cacheChargeSquare = 0;
257
258 // Cached parameters are reset
259 rangeCacheParticle = 0;
260 rangeCacheMatCutsCouple = 0;
261 rangeCacheEnergyRange = 0;
262 rangeCacheRangeEnergy = 0;
263
264 // Cached parameters are reset
265 dedxCacheParticle = 0;
266 dedxCacheMaterial = 0;
267 dedxCacheEnergyCut = 0;
268 dedxCacheIter = lossTableList.end();
269 dedxCacheTransitionEnergy = 0.0;
270 dedxCacheTransitionFactor = 0.0;
271 dedxCacheGenIonMassRatio = 0.0;
272
273 // By default ICRU 73 stopping power tables are loaded:
274 if(!isInitialised) {
276 isInitialised = true;
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
289 // Range vs energy and energy vs range vectors from previous runs are
290 // cleared
291 RangeEnergyTable::iterator iterRange = r.begin();
292 RangeEnergyTable::iterator iterRange_end = r.end();
293
294 for(;iterRange != iterRange_end; ++iterRange) {
295 delete iterRange->second;
296 }
297 r.clear();
298
299 EnergyRangeTable::iterator iterEnergy = E.begin();
300 EnergyRangeTable::iterator iterEnergy_end = E.end();
301
302 for(;iterEnergy != iterEnergy_end; iterEnergy++) {
303 delete iterEnergy->second;
304 }
305 E.clear();
306
307 // The cut energies
308 cutEnergies = cuts;
309
310 // All dE/dx vectors are built
314
315#ifdef PRINT_TABLE_BUILT
317 << " Building dE/dx vectors:"
318 << G4endl;
319#endif
320
321 for (size_t i = 0; i < nmbCouples; ++i) {
322
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 if((*iter)->BuildDEDXTable(atomicNumberIon, material)) {
340
341#ifdef PRINT_TABLE_BUILT
345 << G4endl;
346#endif
347 break;
348 }
349 }
350 }
351 }
352
353 // The particle change object
354 if(! particleChangeLoss) {
355 particleChangeLoss = GetParticleChangeForLoss();
356 braggIonModel->SetParticleChange(particleChangeLoss, 0);
357 betheBlochModel->SetParticleChange(particleChangeLoss, 0);
358 }
359
360 // The G4BraggIonModel and G4BetheBlochModel instances are initialised with
361 // the same settings as the current model:
362 braggIonModel -> Initialise(particle, cuts);
363 betheBlochModel -> Initialise(particle, cuts);
364}
G4bool AddDEDXTable(const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
Definition: G4IonParametrisedLossModel.cc:1237
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4IonParametrisedLossModel.cc:248
Definition: G4IonStoppingData.hh:59
Definition: G4MaterialCutsCouple.hh:53
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
Definition: G4ProductionCutsTable.hh:57
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:246
std::size_t GetTableSize() const
Definition: G4ProductionCutsTable.hh:239
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:58
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
Definition: G4VEmModel.cc:456
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
Referenced by Initialise().
◆ IsApplicable()
|
inline |
Referenced by PrintDEDXTableHandlers().
◆ MaxSecondaryEnergy()
|
protectedvirtual |
Reimplemented from G4VEmModel.
Definition at line 192 of file G4IonParametrisedLossModel.cc.
194 {
195
196 // ############## Maximum energy of secondaries ##########################
197 // Function computes maximum energy of secondary electrons which are
198 // released by an ion
199 //
200 // See Geant4 physics reference manual (version 9.1), section 9.1.1
201 //
202 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
203 // C.Caso et al. (Part. Data Group), Europ. Phys. Jour. C 3 1 (1998).
204 // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
205 //
206 // (Implementation adapted from G4BraggIonModel)
207
208 if(particle != cacheParticle) UpdateCache(particle);
209
210 G4double tau = kineticEnergy/cacheMass;
211 G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
212 (1. + 2.0 * (tau + 1.) * cacheElecMassRatio +
213 cacheElecMassRatio * cacheElecMassRatio);
214
215 return tmax;
216}
Referenced by ComputeCrossSectionPerAtom().
◆ MinEnergyCut()
|
virtual |
Reimplemented from G4VEmModel.
Definition at line 183 of file G4IonParametrisedLossModel.cc.
◆ PrintDEDXTable()
| void G4IonParametrisedLossModel::PrintDEDXTable | ( | const G4ParticleDefinition * | particle, |
| const G4Material * | material, | ||
| G4double | lowerBoundary, | ||
| G4double | upperBoundary, | ||
| G4int | numBins, | ||
| G4bool | logScaleEnergy | ||
| ) |
Definition at line 601 of file G4IonParametrisedLossModel.cc.
607 { // Logarithmic scaling of energy
608
609 G4double atomicMassNumber = particle -> GetAtomicMass();
610 G4double materialDensity = material -> GetDensity();
611
614 << " of density " << materialDensity / g * cm3
615 << " g/cm3"
616 << G4endl
617 << "# Projectile mass number A1 = " << atomicMassNumber
618 << G4endl
619 << "# ------------------------------------------------------"
620 << G4endl;
622 << std::setw(13) << std::right << "E"
623 << std::setw(14) << "E/A1"
624 << std::setw(14) << "dE/dx"
625 << std::setw(14) << "1/rho*dE/dx"
626 << G4endl;
628 << std::setw(13) << std::right << "(MeV)"
629 << std::setw(14) << "(MeV)"
630 << std::setw(14) << "(MeV/cm)"
631 << std::setw(14) << "(MeV*cm2/mg)"
632 << G4endl
633 << "# ------------------------------------------------------"
634 << G4endl;
635
636 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
637 G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
638
639 if(logScaleEnergy) {
640
641 energyLowerBoundary = std::log(energyLowerBoundary);
642 energyUpperBoundary = std::log(energyUpperBoundary);
643 }
644
645 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
646 G4double(nmbBins);
647
648 for(int i = 0; i < numBins + 1; i++) {
649
652
654 G4cout.precision(6);
656 << std::setw(14) << energy / atomicMassNumber / MeV
657 << std::setw(14) << dedx / MeV * cm
658 << std::setw(14) << dedx / materialDensity / (MeV*cm2/(0.001*g))
659 << G4endl;
660 }
661}
Referenced by PrintDEDXTableHandlers().
◆ PrintDEDXTableHandlers()
| void G4IonParametrisedLossModel::PrintDEDXTableHandlers | ( | const G4ParticleDefinition * | particle, |
| const G4Material * | material, | ||
| G4double | lowerBoundary, | ||
| G4double | upperBoundary, | ||
| G4int | numBins, | ||
| G4bool | logScaleEnergy | ||
| ) |
Definition at line 665 of file G4IonParametrisedLossModel.cc.
671 { // Logarithmic scaling of energy
672
673 LossTableList::iterator iter = lossTableList.begin();
674 LossTableList::iterator iter_end = lossTableList.end();
675
676 for(;iter != iter_end; iter++) {
678 if(isApplicable) {
679 (*iter) -> PrintDEDXTable(particle, material,
680 lowerBoundary, upperBoundary,
681 numBins,logScaleEnergy);
682 break;
683 }
684 }
685}
LossTableList::iterator IsApplicable(const G4ParticleDefinition *, const G4Material *)
void PrintDEDXTable(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
Definition: G4IonParametrisedLossModel.cc:601
◆ RemoveDEDXTable()
Definition at line 1280 of file G4IonParametrisedLossModel.cc.
1281 {
1282
1283 LossTableList::iterator iter = lossTableList.begin();
1284 LossTableList::iterator iter_end = lossTableList.end();
1285
1286 for(;iter != iter_end; iter++) {
1288
1289 if(tableName == nam) {
1290 delete (*iter);
1291
1292 // Remove from table list
1293 lossTableList.erase(iter);
1294
1295 // Range vs energy and energy vs range vectors are cleared
1296 RangeEnergyTable::iterator iterRange = r.begin();
1297 RangeEnergyTable::iterator iterRange_end = r.end();
1298
1299 for(;iterRange != iterRange_end; iterRange++)
1300 delete iterRange -> second;
1301 r.clear();
1302
1303 EnergyRangeTable::iterator iterEnergy = E.begin();
1304 EnergyRangeTable::iterator iterEnergy_end = E.end();
1305
1306 for(;iterEnergy != iterEnergy_end; iterEnergy++)
1307 delete iterEnergy -> second;
1308 E.clear();
1309
1310 return true;
1311 }
1312 }
1313
1314 return false;
1315}
◆ SampleSecondaries()
|
virtual |
Implements G4VEmModel.
Definition at line 689 of file G4IonParametrisedLossModel.cc.
694 {
695
696
697 // ############## Sampling of secondaries #################################
698 // The probability density function (pdf) of the kinetic energy T of a
699 // secondary electron may be written as:
700 // pdf(T) = f(T) * g(T)
701 // where
702 // f(T) = (Tmax - Tcut) / (Tmax * Tcut) * (1 / T^2)
703 // g(T) = 1 - beta^2 * T / Tmax
704 // where Tmax is the maximum kinetic energy of the secondary, Tcut
705 // is the lower energy cut and beta is the kinetic energy of the
706 // projectile.
707 //
708 // Sampling of the kinetic energy of a secondary electron:
709 // 1) T0 is sampled from f(T) using the cumulated distribution function
710 // F(T) = int_Tcut^T f(T')dT'
711 // 2) T is accepted or rejected by evaluating the rejection function g(T)
712 // at the sampled energy T0 against a randomly sampled value
713 //
714 //
715 // See Geant4 physics reference manual (version 9.1), section 9.1.4
716 //
717 //
718 // Reference pdf: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
719 //
720 // (Implementation adapted from G4BraggIonModel)
721
723 G4double maxKinEnergySec =
724 std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
725
726 if(cutKinEnergySec >= maxKinEnergySec) return;
727
728 G4double kineticEnergy = particle -> GetKineticEnergy();
729
732 / (energy * energy);
733
734 G4double kinEnergySec;
735 G4double grej;
736
737 do {
738
739 // Sampling kinetic energy from f(T) (using F(T)):
741 kinEnergySec = cutKinEnergySec * maxKinEnergySec /
742 (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
743
744 // Deriving the value of the rejection function at the obtained kinetic
745 // energy:
746 grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
747
748 if(grej > 1.0) {
750 << "Majorant 1.0 < "
751 << grej << " for e= " << kinEnergySec
752 << G4endl;
753 }
754
756
759
761
763 GetAngularDistribution()->SampleDirection(particle, kinEnergySec,
764 Z, mat),
765 kinEnergySec);
766
767
768 secondaries->push_back(delta);
769
770 // Change kinematics of primary particle
772 G4double totalMomentum = std::sqrt(kineticEnergy*(energy + cacheMass));
773
775 finalP = finalP.unit();
776
777 kineticEnergy -= kinEnergySec;
778
779 particleChangeLoss->SetProposedKineticEnergy(kineticEnergy);
780 particleChangeLoss->SetProposedMomentumDirection(finalP);
781}
Definition: ThreeVector.h:36
Hep3Vector unit() const
Definition: G4DynamicParticle.hh:65
const G4ThreeVector & GetMomentumDirection() const
G4ThreeVector GetMomentum() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForLoss.hh:135
void SetProposedMomentumDirection(const G4ThreeVector &dir)
Definition: G4ParticleChangeForLoss.hh:185
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:611
G4int SelectRandomAtomNumber(const G4Material *)
Definition: G4VEmModel.cc:315
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:510
◆ SetEnergyLossLimit()
|
inline |
The documentation for this class was generated from the following files:
