Loading...
Searching...
No Matches
G4IonParametrisedLossModel Class Reference
#include <G4IonParametrisedLossModel.hh>
Inheritance diagram for G4IonParametrisedLossModel:
Protected Member Functions | |
| G4double | MaxSecondaryEnergy (const G4ParticleDefinition *, G4double) override |
| Protected Member Functions inherited from G4VEmModel | |
| G4ParticleChangeForLoss * | GetParticleChangeForLoss () |
| G4ParticleChangeForGamma * | GetParticleChangeForGamma () |
| const G4MaterialCutsCouple * | CurrentCouple () const |
| void | SetCurrentElement (const G4Element *) |
Additional Inherited Members | |
| Protected Attributes inherited from G4VEmModel | |
| G4ElementData * | fElementData = nullptr |
| G4VParticleChange * | pParticleChange = nullptr |
| G4PhysicsTable * | xSectionTable = nullptr |
| const G4Material * | pBaseMaterial = nullptr |
| const std::vector< G4double > * | theDensityFactor = nullptr |
| const std::vector< G4int > * | theDensityIdx = nullptr |
| G4double | inveplus |
| G4double | pFactor = 1.0 |
| std::size_t | currentCoupleIndex = 0 |
| std::size_t | basedCoupleIndex = 0 |
| G4bool | lossFlucFlag = true |
Detailed Description
Definition at line 91 of file G4IonParametrisedLossModel.hh.
Constructor & Destructor Documentation
◆ G4IonParametrisedLossModel()
|
explicit |
Definition at line 101 of file G4IonParametrisedLossModel.cc.
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();
118
119 // The Bragg ion and Bethe Bloch models are instantiated
120 braggIonModel = new G4BraggIonModel();
121 betheBlochModel = new G4BetheBlochModel();
122
123 // The boundaries for the range tables are set
124 lowerEnergyEdgeIntegr = 0.025 * MeV;
125 upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit();
126
127 // Cache parameters are set
128 cacheParticle = 0;
129 cacheMass = 0;
130 cacheElecMassRatio = 0;
131 cacheChargeSquare = 0;
132
133 // Cache parameters are set
134 rangeCacheParticle = 0;
135 rangeCacheMatCutsCouple = 0;
136 rangeCacheEnergyRange = 0;
137 rangeCacheRangeEnergy = 0;
138
139 // Cache parameters are set
140 dedxCacheParticle = 0;
141 dedxCacheMaterial = 0;
142 dedxCacheEnergyCut = 0;
143 dedxCacheIter = lossTableList.end();
144 dedxCacheTransitionEnergy = 0.0;
145 dedxCacheTransitionFactor = 0.0;
146 dedxCacheGenIonMassRatio = 0.0;
147
148 // default generator
150}
static G4LossTableManager * Instance()
Definition G4LossTableManager.cc:88
G4EmCorrections * EmCorrections()
Definition G4LossTableManager.hh:383
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition G4VEmModel.hh:606
Referenced by SetEnergyLossLimit().
◆ ~G4IonParametrisedLossModel()
|
virtual |
Definition at line 154 of file G4IonParametrisedLossModel.cc.
154 {
155
156 // dE/dx table objects are deleted and the container is cleared
157 LossTableList::iterator iterTables = lossTableList.begin();
158 LossTableList::iterator iterTables_end = lossTableList.end();
159
160 for(;iterTables != iterTables_end; ++iterTables) { delete *iterTables; }
161 lossTableList.clear();
162
163 // range table
164 RangeEnergyTable::iterator itr = r.begin();
165 RangeEnergyTable::iterator itr_end = r.end();
166 for(;itr != itr_end; ++itr) { delete itr->second; }
167 r.clear();
168
169 // inverse range
170 EnergyRangeTable::iterator ite = E.begin();
171 EnergyRangeTable::iterator ite_end = E.end();
172 for(;ite != ite_end; ++ite) { delete ite->second; }
173 E.clear();
174
175}
Member Function Documentation
◆ AddDEDXTable()
| G4bool G4IonParametrisedLossModel::AddDEDXTable | ( | const G4String & | name, |
| G4VIonDEDXTable * | table, | ||
| G4VIonDEDXScalingAlgorithm * | algorithm = nullptr ) |
Definition at line 1232 of file G4IonParametrisedLossModel.cc.
1235 {
1236
1237 if(table == 0) {
1239 << " add table: Invalid pointer."
1240 << G4endl;
1241
1242 return false;
1243 }
1244
1245 // Checking uniqueness of name
1246 LossTableList::iterator iter = lossTableList.begin();
1247 LossTableList::iterator iter_end = lossTableList.end();
1248
1249 for(;iter != iter_end; ++iter) {
1250 const G4String& tableName = (*iter)->GetName();
1251
1252 if(tableName == nam) {
1254 << " add table: Name already exists."
1255 << G4endl;
1256
1257 return false;
1258 }
1259 }
1260
1261 G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm;
1262 if(scalingAlgorithm == 0)
1263 scalingAlgorithm = new G4VIonDEDXScalingAlgorithm;
1264
1265 G4IonDEDXHandler* handler =
1266 new G4IonDEDXHandler(table, scalingAlgorithm, nam);
1267
1268 lossTableList.push_front(handler);
1269
1270 return true;
1271}
G4GLOB_DLL std::ostream G4cout
Referenced by Initialise().
◆ ComputeCrossSectionPerAtom()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 364 of file G4IonParametrisedLossModel.cc.
370 {
371
372 // ############## Cross section per atom ################################
373 // Function computes ionization cross section per atom
374 //
375 // See Geant4 physics reference manual (version 9.1), section 9.1.3
376 //
377 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
378 // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
379 //
380 // (Implementation adapted from G4BraggIonModel)
381
382 G4double crosssection = 0.0;
384 G4double maxEnergy = std::min(tmax, maxKinEnergy);
385
386 if(cutEnergy < tmax) {
387
389 G4double betaSquared = kineticEnergy *
390 (energy + cacheMass) / (energy * energy);
391
392 crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
393 betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
394
395 crosssection *= twopi_mc2_rcl2 * cacheChargeSquare / betaSquared;
396 }
397
398#ifdef PRINT_DEBUG_CS
400 << G4endl
401 << "# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
402 << G4endl
403 << "# particle =" << particle -> GetParticleName()
404 << G4endl
405 << "# cut(MeV) = " << cutEnergy/MeV
406 << G4endl;
407
409 << std::setw(13) << std::right << "E(MeV)"
410 << std::setw(14) << "CS(um)"
411 << std::setw(14) << "E_max_sec(MeV)"
412 << G4endl
413 << "# ------------------------------------------------------"
414 << G4endl;
415
416 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
417 << std::setw(14) << crosssection / (um * um)
418 << std::setw(14) << tmax / MeV
419 << G4endl;
420#endif
421
422 crosssection *= atomicNumber;
423
424 return crosssection;
425}
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double) override
Definition G4IonParametrisedLossModel.cc:188
G4double energy(const ThreeVector &p, const G4double m)
Definition G4INCLKinematicsUtils.cc:185
Referenced by CrossSectionPerVolume().
◆ ComputeDEDXPerVolume()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 448 of file G4IonParametrisedLossModel.cc.
452 {
453
454 // ############## dE/dx ##################################################
455 // Function computes dE/dx values, where following rules are adopted:
456 // A. If the ion-material pair is covered by any native ion data
457 // parameterisation, then:
458 // * This parameterization is used for energies below a given energy
459 // limit,
460 // * whereas above the limit the Bethe-Bloch model is applied, in
461 // combination with an effective charge estimate and high order
462 // correction terms.
463 // A smoothing procedure is applied to dE/dx values computed with
464 // the second approach. The smoothing factor is based on the dE/dx
465 // values of both approaches at the transition energy (high order
466 // correction terms are included in the calculation of the transition
467 // factor).
468 // B. If the particle is a generic ion, the BraggIon and Bethe-Bloch
469 // models are used and a smoothing procedure is applied to values
470 // obtained with the second approach.
471 // C. If the ion-material is not covered by any ion data parameterization
472 // then:
473 // * The BraggIon model is used for energies below a given energy
474 // limit,
475 // * whereas above the limit the Bethe-Bloch model is applied, in
476 // combination with an effective charge estimate and high order
477 // correction terms.
478 // Also in this case, a smoothing procedure is applied to dE/dx values
479 // computed with the second model.
480
481 G4double dEdx = 0.0;
482 UpdateDEDXCache(particle, material, cutEnergy);
483
484 LossTableList::iterator iter = dedxCacheIter;
485
486 if(iter != lossTableList.end()) {
487
488 G4double transitionEnergy = dedxCacheTransitionEnergy;
489
490 if(transitionEnergy > kineticEnergy) {
491
492 dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
493
495 particle,
496 kineticEnergy,
497 cutEnergy);
498 dEdx -= dEdxDeltaRays;
499 }
500 else {
501 G4double massRatio = dedxCacheGenIonMassRatio;
502
503 G4double chargeSquare =
504 GetChargeSquareRatio(particle, material, kineticEnergy);
505
506 G4double scaledKineticEnergy = kineticEnergy * massRatio;
507 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
508
510
511 if(scaledTransitionEnergy >= lowEnergyLimit) {
512
513 dEdx = betheBlochModel -> ComputeDEDXPerVolume(
514 material, genericIon,
515 scaledKineticEnergy, cutEnergy);
516
517 dEdx *= chargeSquare;
518
519 dEdx += corrections -> ComputeIonCorrections(particle,
520 material, kineticEnergy);
521
522 G4double factor = 1.0 + dedxCacheTransitionFactor /
523 kineticEnergy;
524
525 dEdx *= factor;
526 }
527 }
528 }
529 else {
530 G4double massRatio = 1.0;
531 G4double chargeSquare = 1.0;
532
533 if(particle != genericIon) {
534
535 chargeSquare = GetChargeSquareRatio(particle, material, kineticEnergy);
536 massRatio = genericIonPDGMass / particle -> GetPDGMass();
537 }
538
539 G4double scaledKineticEnergy = kineticEnergy * massRatio;
540
542 if(scaledKineticEnergy < lowEnergyLimit) {
543 dEdx = braggIonModel -> ComputeDEDXPerVolume(
544 material, genericIon,
545 scaledKineticEnergy, cutEnergy);
546
547 dEdx *= chargeSquare;
548 }
549 else {
551 material, genericIon,
552 lowEnergyLimit, cutEnergy);
553
555 material, genericIon,
556 lowEnergyLimit, cutEnergy);
557
558 if(particle != genericIon) {
559 G4double chargeSquareLowEnergyLimit =
560 GetChargeSquareRatio(particle, material,
561 lowEnergyLimit / massRatio);
562
563 dEdxLimitParam *= chargeSquareLowEnergyLimit;
564 dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
565
566 dEdxLimitBetheBloch +=
567 corrections -> ComputeIonCorrections(particle,
568 material, lowEnergyLimit / massRatio);
569 }
570
571 G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
572 * lowEnergyLimit / scaledKineticEnergy);
573
574 dEdx = betheBlochModel -> ComputeDEDXPerVolume(
575 material, genericIon,
576 scaledKineticEnergy, cutEnergy);
577
578 dEdx *= chargeSquare;
579
580 if(particle != genericIon) {
581 dEdx += corrections -> ComputeIonCorrections(particle,
582 material, kineticEnergy);
583 }
584
585 dEdx *= factor;
586 }
587
588 }
589
590 if (dEdx < 0.0) dEdx = 0.0;
591
592 return dEdx;
593}
G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double) override
Definition G4IonParametrisedLossModel.cc:216
G4double DeltaRayMeanEnergyTransferRate(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double) override
Definition G4IonParametrisedLossModel.cc:448
Referenced by ComputeDEDXPerVolume(), CorrectionsAlongStep(), and PrintDEDXTable().
◆ ComputeLossForStep()
| G4double G4IonParametrisedLossModel::ComputeLossForStep | ( | const G4MaterialCutsCouple * | matCutsCouple, |
| const G4ParticleDefinition * | particle, | ||
| G4double | kineticEnergy, | ||
| G4double | stepLength ) |
Definition at line 1173 of file G4IonParametrisedLossModel.cc.
1177 {
1178
1179 G4double loss = 0.0;
1180
1181 UpdateRangeCache(particle, matCutsCouple);
1182
1183 G4PhysicsVector* energyRange = rangeCacheEnergyRange;
1184 G4PhysicsVector* rangeEnergy = rangeCacheRangeEnergy;
1185
1186 if(energyRange != 0 && rangeEnergy != 0) {
1187
1188 G4double lowerEnEdge = energyRange -> Energy( 0 );
1189 G4double lowerRangeEdge = rangeEnergy -> Energy( 0 );
1190
1191 // Computing range for pre-step kinetic energy:
1193
1194 // Energy below vector boundary:
1195 if(kineticEnergy < lowerEnEdge) {
1196
1197 range = energyRange -> Value(lowerEnEdge);
1198 range *= std::sqrt(kineticEnergy / lowerEnEdge);
1199 }
1200
1201#ifdef PRINT_DEBUG
1203 << range / mm << " mm, step = " << stepLength / mm << " mm"
1204 << G4endl;
1205#endif
1206
1207 // Remaining range:
1208 G4double remRange = range - stepLength;
1209
1210 // If range is smaller than step length, the loss is set to kinetic
1211 // energy
1212 if(remRange < 0.0) loss = kineticEnergy;
1213 else if(remRange < lowerRangeEdge) {
1214
1215 G4double ratio = remRange / lowerRangeEdge;
1216 loss = kineticEnergy - ratio * ratio * lowerEnEdge;
1217 }
1218 else {
1219
1221 loss = kineticEnergy - energy;
1222 }
1223 }
1224
1225 if(loss < 0.0) loss = 0.0;
1226
1227 return loss;
1228}
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition G4VEmModel.cc:347
Referenced by CorrectionsAlongStep().
◆ CorrectionsAlongStep()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 905 of file G4IonParametrisedLossModel.cc.
909 {
910
911 // ############## Corrections for along step energy loss calculation ######
912 // The computed energy loss (due to electronic stopping) is overwritten
913 // by this function if an ion data parameterization is available for the
914 // current ion-material pair.
915 // No action on the energy loss (due to electronic stopping) is performed
916 // if no parameterization is available. In this case the original
917 // generic ion tables (in combination with the effective charge) are used
918 // in the along step DoIt function.
919 //
920 //
921 // (Implementation partly adapted from G4BraggIonModel/G4BetheBlochModel)
922
923 const G4ParticleDefinition* particle = dynamicParticle -> GetDefinition();
924 const G4Material* material = couple -> GetMaterial();
925
926 G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
927
928 if(kineticEnergy == eloss) { return; }
929
931 std::size_t cutIndex = couple -> GetIndex();
932 cutEnergy = cutEnergies[cutIndex];
933
934 UpdateDEDXCache(particle, material, cutEnergy);
935
936 LossTableList::iterator iter = dedxCacheIter;
937
938 // If parameterization for ions is available the electronic energy loss
939 // is overwritten
940 if(iter != lossTableList.end()) {
941
942 // The energy loss is calculated using the ComputeDEDXPerVolume function
943 // and the step length (it is assumed that dE/dx does not change
944 // considerably along the step)
945 eloss =
946 length * ComputeDEDXPerVolume(material, particle,
947 kineticEnergy, cutEnergy);
948
949#ifdef PRINT_DEBUG
950 G4cout.precision(6);
952 << G4endl
953 << "# G4IonParametrisedLossModel::CorrectionsAlongStep"
954 << G4endl
955 << "# cut(MeV) = " << cutEnergy/MeV
956 << G4endl;
957
959 << std::setw(13) << std::right << "E(MeV)"
960 << std::setw(14) << "l(um)"
961 << std::setw(14) << "l*dE/dx(MeV)"
962 << std::setw(14) << "(l*dE/dx)/E"
963 << G4endl
964 << "# ------------------------------------------------------"
965 << G4endl;
966
967 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
968 << std::setw(14) << length / um
969 << std::setw(14) << eloss / MeV
970 << std::setw(14) << eloss / kineticEnergy * 100.0
971 << G4endl;
972#endif
973
974 // If the energy loss exceeds a certain fraction of the kinetic energy
975 // (the fraction is indicated by the parameter "energyLossLimit") then
976 // the range tables are used to derive a more accurate value of the
977 // energy loss
978 if(eloss > energyLossLimit * kineticEnergy) {
979
980 eloss = ComputeLossForStep(couple, particle,
981 kineticEnergy,length);
982
983#ifdef PRINT_DEBUG
985 << G4endl;
986
987 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
988 << std::setw(14) << length / um
989 << std::setw(14) << eloss / MeV
990 << std::setw(14) << eloss / kineticEnergy * 100.0
991 << G4endl;
992#endif
993
994 }
995 }
996
997 // For all corrections below a kinetic energy between the Pre- and
998 // Post-step energy values is used
1001
1002 G4double chargeSquareRatio = corrections ->
1003 EffectiveChargeSquareRatio(particle,
1004 material,
1005 energy);
1006 GetModelOfFluctuations() -> SetParticleAndCharge(particle,
1007 chargeSquareRatio);
1008
1009 // A correction is applied considering the change of the effective charge
1010 // along the step (the parameter "corrFactor" refers to the effective
1011 // charge at the beginning of the step). Note: the correction is not
1012 // applied for energy loss values deriving directly from parameterized
1013 // ion stopping power tables
1014 G4double transitionEnergy = dedxCacheTransitionEnergy;
1015
1016 if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
1017 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1018 material,
1019 energy);
1020
1021 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1022 eloss *= chargeSquareRatioCorr;
1023 }
1024 else if (iter == lossTableList.end()) {
1025
1026 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1027 material,
1028 energy);
1029
1030 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1031 eloss *= chargeSquareRatioCorr;
1032 }
1033
1034 // Ion high order corrections are applied if the current model does not
1035 // overwrite the energy loss (i.e. when the effective charge approach is
1036 // used)
1037 if(iter == lossTableList.end()) {
1038
1039 G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
1041
1042 // Corrections are only applied in the Bethe-Bloch energy region
1043 if(scaledKineticEnergy > lowEnergyLimit)
1044 eloss += length *
1045 corrections -> IonHighOrderCorrections(particle, couple, energy);
1046 }
1047}
G4double ComputeLossForStep(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double, G4double)
Definition G4IonParametrisedLossModel.cc:1173
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition G4VEmModel.hh:592
◆ CrossSectionPerVolume()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 429 of file G4IonParametrisedLossModel.cc.
434 {
435
436 G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
438 kineticEnergy,
439 nbElecPerVolume, 0,
440 cutEnergy,
441 maxEnergy);
442
443 return cross;
444}
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double) override
Definition G4IonParametrisedLossModel.cc:364
◆ DeltaRayMeanEnergyTransferRate()
|
inline |
Referenced by ComputeDEDXPerVolume().
◆ GetChargeSquareRatio()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 216 of file G4IonParametrisedLossModel.cc.
219 { // Kinetic energy
220
221 G4double chargeSquareRatio = corrections ->
222 EffectiveChargeSquareRatio(particle,
223 material,
224 kineticEnergy);
225 corrFactor = chargeSquareRatio *
226 corrections -> EffectiveChargeCorrection(particle,
227 material,
228 kineticEnergy);
229 return corrFactor;
230}
Referenced by ComputeDEDXPerVolume().
◆ GetParticleCharge()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 234 of file G4IonParametrisedLossModel.cc.
237 { // Kinetic energy
238
240}
G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double) override
Definition G4IonParametrisedLossModel.cc:234
Referenced by GetParticleCharge().
◆ Initialise()
|
overridevirtual |
Implements G4VEmModel.
Definition at line 244 of file G4IonParametrisedLossModel.cc.
246 {
247
248 // Cached parameters are reset
249 cacheParticle = 0;
250 cacheMass = 0;
251 cacheElecMassRatio = 0;
252 cacheChargeSquare = 0;
253
254 // Cached parameters are reset
255 rangeCacheParticle = 0;
256 rangeCacheMatCutsCouple = 0;
257 rangeCacheEnergyRange = 0;
258 rangeCacheRangeEnergy = 0;
259
260 // Cached parameters are reset
261 dedxCacheParticle = 0;
262 dedxCacheMaterial = 0;
263 dedxCacheEnergyCut = 0;
264 dedxCacheIter = lossTableList.end();
265 dedxCacheTransitionEnergy = 0.0;
266 dedxCacheTransitionFactor = 0.0;
267 dedxCacheGenIonMassRatio = 0.0;
268
269 // By default ICRU 73 stopping power tables are loaded:
270 if(!isInitialised) {
272 isInitialised = true;
274 new G4IonStoppingData("ion_stopping_data/icru",icru90),
275 new G4IonDEDXScalingICRU73());
276 }
277 // The cache of loss tables is cleared
278 LossTableList::iterator iterTables = lossTableList.begin();
279 LossTableList::iterator iterTables_end = lossTableList.end();
280
281 for(;iterTables != iterTables_end; ++iterTables) {
282 (*iterTables) -> ClearCache();
283 }
284
285 // Range vs energy and energy vs range vectors from previous runs are
286 // cleared
287 RangeEnergyTable::iterator iterRange = r.begin();
288 RangeEnergyTable::iterator iterRange_end = r.end();
289
290 for(;iterRange != iterRange_end; ++iterRange) {
291 delete iterRange->second;
292 }
293 r.clear();
294
295 EnergyRangeTable::iterator iterEnergy = E.begin();
296 EnergyRangeTable::iterator iterEnergy_end = E.end();
297
298 for(;iterEnergy != iterEnergy_end; ++iterEnergy) {
299 delete iterEnergy->second;
300 }
301 E.clear();
302
303 // The cut energies
304 cutEnergies = cuts;
305
306 // All dE/dx vectors are built
307 const G4ProductionCutsTable* coupleTable=
310
311#ifdef PRINT_TABLE_BUILT
313 << " Building dE/dx vectors:"
314 << G4endl;
315#endif
316
318
321
323
324 LossTableList::iterator iter = lossTableList.begin();
325 LossTableList::iterator iter_end = lossTableList.end();
326
327 for(;iter != iter_end; ++iter) {
328
329 if(*iter == 0) {
331 << " Skipping illegal table."
332 << G4endl;
333 }
334
335 if((*iter)->BuildDEDXTable(atomicNumberIon, material)) {
336
337#ifdef PRINT_TABLE_BUILT
341 << G4endl;
342#endif
343 break;
344 }
345 }
346 }
347 }
348
349 // The particle change object
350 if(! particleChangeLoss) {
351 particleChangeLoss = GetParticleChangeForLoss();
352 braggIonModel->SetParticleChange(particleChangeLoss, 0);
353 betheBlochModel->SetParticleChange(particleChangeLoss, 0);
354 }
355
356 // The G4BraggIonModel and G4BetheBlochModel instances are initialised with
357 // the same settings as the current model:
358 braggIonModel -> Initialise(particle, cuts);
359 betheBlochModel -> Initialise(particle, cuts);
360}
G4bool AddDEDXTable(const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=nullptr)
Definition G4IonParametrisedLossModel.cc:1232
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition G4IonParametrisedLossModel.cc:244
const G4Material * GetMaterial() const
Definition G4MaterialCutsCouple.hh:166
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition G4ProductionCutsTable.hh:256
std::size_t GetTableSize() const
Definition G4ProductionCutsTable.hh:249
static G4ProductionCutsTable * GetProductionCutsTable()
Definition G4ProductionCutsTable.cc:58
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition G4VEmModel.cc:106
Referenced by Initialise().
◆ IsApplicable()
|
inline |
Referenced by PrintDEDXTableHandlers().
◆ MaxSecondaryEnergy()
|
overrideprotectedvirtual |
Reimplemented from G4VEmModel.
Definition at line 188 of file G4IonParametrisedLossModel.cc.
190 {
191
192 // ############## Maximum energy of secondaries ##########################
193 // Function computes maximum energy of secondary electrons which are
194 // released by an ion
195 //
196 // See Geant4 physics reference manual (version 9.1), section 9.1.1
197 //
198 // Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
199 // C.Caso et al. (Part. Data Group), Europ. Phys. Jour. C 3 1 (1998).
200 // B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
201 //
202 // (Implementation adapted from G4BraggIonModel)
203
204 if(particle != cacheParticle) UpdateCache(particle);
205
206 G4double tau = kineticEnergy/cacheMass;
207 G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
208 (1. + 2.0 * (tau + 1.) * cacheElecMassRatio +
209 cacheElecMassRatio * cacheElecMassRatio);
210
211 return tmax;
212}
Referenced by ComputeCrossSectionPerAtom().
◆ MinEnergyCut()
|
overridevirtual |
Reimplemented from G4VEmModel.
Definition at line 179 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 597 of file G4IonParametrisedLossModel.cc.
603 { // Logarithmic scaling of energy
604
605 G4double atomicMassNumber = particle -> GetAtomicMass();
606 G4double materialDensity = material -> GetDensity();
607
610 << " of density " << materialDensity / g * cm3
611 << " g/cm3"
612 << G4endl
613 << "# Projectile mass number A1 = " << atomicMassNumber
614 << G4endl
615 << "# ------------------------------------------------------"
616 << G4endl;
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;
624 << std::setw(13) << std::right << "(MeV)"
625 << std::setw(14) << "(MeV)"
626 << std::setw(14) << "(MeV/cm)"
627 << std::setw(14) << "(MeV*cm2/mg)"
628 << G4endl
629 << "# ------------------------------------------------------"
630 << G4endl;
631
632 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
633 G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
634
635 if(logScaleEnergy) {
636
637 energyLowerBoundary = std::log(energyLowerBoundary);
638 energyUpperBoundary = std::log(energyUpperBoundary);
639 }
640
641 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
642 G4double(nmbBins);
643
644 for(int i = 0; i < numBins + 1; i++) {
645
648
650 G4cout.precision(6);
652 << std::setw(14) << energy / atomicMassNumber / MeV
653 << std::setw(14) << dedx / MeV * cm
654 << std::setw(14) << dedx / materialDensity / (MeV*cm2/(0.001*g))
655 << G4endl;
656 }
657}
Referenced by PrintDEDXTableHandlers().
◆ PrintDEDXTableHandlers()
| void G4IonParametrisedLossModel::PrintDEDXTableHandlers | ( | const G4ParticleDefinition * | particle, |
| const G4Material * | material, | ||
| G4double | lowerBoundary, | ||
| G4double | upperBoundary, | ||
| G4int | numBins, | ||
| G4bool | logScaleEnergy ) |
Definition at line 661 of file G4IonParametrisedLossModel.cc.
667 { // Logarithmic scaling of energy
668
669 LossTableList::iterator iter = lossTableList.begin();
670 LossTableList::iterator iter_end = lossTableList.end();
671
672 for(;iter != iter_end; iter++) {
674 if(isApplicable) {
675 (*iter) -> PrintDEDXTable(particle, material,
676 lowerBoundary, upperBoundary,
677 numBins,logScaleEnergy);
678 break;
679 }
680 }
681}
LossTableList::iterator IsApplicable(const G4ParticleDefinition *, const G4Material *)
void PrintDEDXTable(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
Definition G4IonParametrisedLossModel.cc:597
◆ RemoveDEDXTable()
Definition at line 1275 of file G4IonParametrisedLossModel.cc.
1276 {
1277
1278 LossTableList::iterator iter = lossTableList.begin();
1279 LossTableList::iterator iter_end = lossTableList.end();
1280
1281 for(;iter != iter_end; iter++) {
1282 G4String tableName = (*iter) -> GetName();
1283
1284 if(tableName == nam) {
1285 delete (*iter);
1286
1287 // Remove from table list
1288 lossTableList.erase(iter);
1289
1290 // Range vs energy and energy vs range vectors are cleared
1291 RangeEnergyTable::iterator iterRange = r.begin();
1292 RangeEnergyTable::iterator iterRange_end = r.end();
1293
1294 for(;iterRange != iterRange_end; iterRange++)
1295 delete iterRange -> second;
1296 r.clear();
1297
1298 EnergyRangeTable::iterator iterEnergy = E.begin();
1299 EnergyRangeTable::iterator iterEnergy_end = E.end();
1300
1301 for(;iterEnergy != iterEnergy_end; iterEnergy++)
1302 delete iterEnergy -> second;
1303 E.clear();
1304
1305 return true;
1306 }
1307 }
1308
1309 return false;
1310}
◆ SampleSecondaries()
|
overridevirtual |
Implements G4VEmModel.
Definition at line 685 of file G4IonParametrisedLossModel.cc.
690 {
691 // ############## Sampling of secondaries #################################
692 // The probability density function (pdf) of the kinetic energy T of a
693 // secondary electron may be written as:
694 // pdf(T) = f(T) * g(T)
695 // where
696 // f(T) = (Tmax - Tcut) / (Tmax * Tcut) * (1 / T^2)
697 // g(T) = 1 - beta^2 * T / Tmax
698 // where Tmax is the maximum kinetic energy of the secondary, Tcut
699 // is the lower energy cut and beta is the kinetic energy of the
700 // projectile.
701 //
702 // Sampling of the kinetic energy of a secondary electron:
703 // 1) T0 is sampled from f(T) using the cumulated distribution function
704 // F(T) = int_Tcut^T f(T')dT'
705 // 2) T is accepted or rejected by evaluating the rejection function g(T)
706 // at the sampled energy T0 against a randomly sampled value
707 //
708 //
709 // See Geant4 physics reference manual (version 9.1), section 9.1.4
710 //
711 //
712 // Reference pdf: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
713 //
714 // (Implementation adapted from G4BraggIonModel)
715
717 G4double maxKinEnergySec =
718 std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
719
720 if(cutKinEnergySec >= maxKinEnergySec) return;
721
722 G4double kineticEnergy = particle -> GetKineticEnergy();
723
726 / (energy * energy);
727
728 G4double kinEnergySec;
729 G4double grej;
730
731 do {
732
733 // Sampling kinetic energy from f(T) (using F(T)):
735 kinEnergySec = cutKinEnergySec * maxKinEnergySec /
736 (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
737
738 // Deriving the value of the rejection function at the obtained kinetic
739 // energy:
740 grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
741
742 if(grej > 1.0) {
744 << "Majorant 1.0 < "
745 << grej << " for e= " << kinEnergySec
746 << G4endl;
747 }
748
750
753
755
756 G4DynamicParticle* delta = new G4DynamicParticle(electron,
757 GetAngularDistribution()->SampleDirection(particle, kinEnergySec,
758 Z, mat),
759 kinEnergySec);
760
761 secondaries->push_back(delta);
762
763 // Change kinematics of primary particle
765 G4double totalMomentum = std::sqrt(kineticEnergy*(energy + cacheMass));
766
768 finalP = finalP.unit();
769
770 kineticEnergy -= kinEnergySec;
771
772 particleChangeLoss->SetProposedKineticEnergy(kineticEnergy);
773 particleChangeLoss->SetProposedMomentumDirection(finalP);
774}
Hep3Vector unit() const
const G4ThreeVector & GetMomentumDirection() const
G4ThreeVector GetMomentum() const
G4VEmAngularDistribution * GetAngularDistribution()
Definition G4VEmModel.hh:599
G4int SelectRandomAtomNumber(const G4Material *) const
Definition G4VEmModel.cc:251
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition G4VEmModel.hh:500
◆ SetEnergyLossLimit()
|
inline |
The documentation for this class was generated from the following files:
