111 particleChangeLoss(0),
113 energyLossLimit(0.01),
117 genericIonPDGMass = genericIon -> GetPDGMass();
133 lowerEnergyEdgeIntegr = 0.025 *
MeV;
139 cacheElecMassRatio = 0;
140 cacheChargeSquare = 0;
143 rangeCacheParticle = 0;
144 rangeCacheMatCutsCouple = 0;
145 rangeCacheEnergyRange = 0;
146 rangeCacheRangeEnergy = 0;
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;
163 LossTableList::iterator iterTables = lossTableList.begin();
164 LossTableList::iterator iterTables_end = lossTableList.end();
166 for(;iterTables != iterTables_end; iterTables++)
delete *iterTables;
167 lossTableList.clear();
176 return couple -> GetMaterial() -> GetIonisation() ->
177 GetMeanExcitationEnergy();
198 if(particle != cacheParticle) UpdateCache(particle);
200 G4double tau = kineticEnergy/cacheMass;
202 (1. + 2.0 * (tau + 1.) * cacheElecMassRatio +
203 cacheElecMassRatio * cacheElecMassRatio);
215 G4double chargeSquareRatio = corrections ->
216 EffectiveChargeSquareRatio(particle,
219 corrFactor = chargeSquareRatio *
220 corrections -> EffectiveChargeCorrection(particle,
245 cacheElecMassRatio = 0;
246 cacheChargeSquare = 0;
249 rangeCacheParticle = 0;
250 rangeCacheMatCutsCouple = 0;
251 rangeCacheEnergyRange = 0;
252 rangeCacheRangeEnergy = 0;
255 dedxCacheParticle = 0;
256 dedxCacheMaterial = 0;
257 dedxCacheEnergyCut = 0;
258 dedxCacheIter = lossTableList.end();
259 dedxCacheTransitionEnergy = 0.0;
260 dedxCacheTransitionFactor = 0.0;
261 dedxCacheGenIonMassRatio = 0.0;
264 LossTableList::iterator iterTables = lossTableList.begin();
265 LossTableList::iterator iterTables_end = lossTableList.end();
267 for(;iterTables != iterTables_end; iterTables++)
268 (*iterTables) -> ClearCache();
272 RangeEnergyTable::iterator iterRange = r.begin();
273 RangeEnergyTable::iterator iterRange_end = r.end();
275 for(;iterRange != iterRange_end; iterRange++)
delete iterRange ->
second;
278 EnergyRangeTable::iterator iterEnergy = E.begin();
279 EnergyRangeTable::iterator iterEnergy_end = E.end();
281 for(;iterEnergy != iterEnergy_end; iterEnergy++)
delete iterEnergy ->
second;
285 size_t size = cuts.size();
287 for(
size_t i = 0; i < size; i++) cutEnergies.push_back(cuts[i]);
292 size_t nmbCouples = coupleTable -> GetTableSize();
294 #ifdef PRINT_TABLE_BUILT
295 G4cout <<
"G4IonParametrisedLossModel::Initialise():"
296 <<
" Building dE/dx vectors:"
300 for (
size_t i = 0; i < nmbCouples; i++) {
303 coupleTable -> GetMaterialCutsCouple(i);
308 for(
G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) {
310 LossTableList::iterator iter = lossTableList.begin();
311 LossTableList::iterator iter_end = lossTableList.end();
313 for(;iter != iter_end; iter++) {
316 G4cout <<
"G4IonParametrisedLossModel::Initialise():"
317 <<
" Skipping illegal table."
322 (*iter) -> BuildDEDXTable(atomicNumberIon, material);
325 #ifdef PRINT_TABLE_BUILT
326 G4cout <<
" Atomic Number Ion = " << atomicNumberIon
327 <<
", Material = " << material ->
GetName()
328 <<
", Table = " << (*iter) ->
GetName()
338 if(! particleChangeLoss) {
347 betheBlochModel ->
Initialise(particle, cuts);
374 if(cutEnergy < tmax) {
377 G4double betaSquared = kineticEnergy *
378 (energy + cacheMass) / (energy * energy);
380 crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
381 betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
386 #ifdef PRINT_DEBUG_CS
387 G4cout <<
"########################################################"
389 <<
"# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
391 <<
"# particle =" << particle -> GetParticleName()
393 <<
"# cut(MeV) = " << cutEnergy/
MeV
398 << std::setw(14) <<
"CS(um)"
399 << std::setw(14) <<
"E_max_sec(MeV)"
401 <<
"# ------------------------------------------------------"
405 << std::setw(14) << crosssection / (um * um)
406 << std::setw(14) << tmax /
MeV
410 crosssection *= atomicNumber;
424 G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
470 UpdateDEDXCache(particle, material, cutEnergy);
472 LossTableList::iterator iter = dedxCacheIter;
474 if(iter != lossTableList.end()) {
476 G4double transitionEnergy = dedxCacheTransitionEnergy;
478 if(transitionEnergy > kineticEnergy) {
480 dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
486 dEdx -= dEdxDeltaRays;
489 G4double massRatio = dedxCacheGenIonMassRatio;
494 G4double scaledKineticEnergy = kineticEnergy * massRatio;
495 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
499 if(scaledTransitionEnergy >= lowEnergyLimit) {
502 material, genericIon,
503 scaledKineticEnergy, cutEnergy);
505 dEdx *= chargeSquare;
507 dEdx += corrections -> ComputeIonCorrections(particle,
508 material, kineticEnergy);
510 G4double factor = 1.0 + dedxCacheTransitionFactor /
521 if(particle != genericIon) {
524 massRatio = genericIonPDGMass / particle -> GetPDGMass();
527 G4double scaledKineticEnergy = kineticEnergy * massRatio;
530 if(scaledKineticEnergy < lowEnergyLimit) {
532 material, genericIon,
533 scaledKineticEnergy, cutEnergy);
535 dEdx *= chargeSquare;
539 material, genericIon,
540 lowEnergyLimit, cutEnergy);
543 material, genericIon,
544 lowEnergyLimit, cutEnergy);
546 if(particle != genericIon) {
547 G4double chargeSquareLowEnergyLimit =
549 lowEnergyLimit / massRatio);
551 dEdxLimitParam *= chargeSquareLowEnergyLimit;
552 dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
554 dEdxLimitBetheBloch +=
555 corrections -> ComputeIonCorrections(particle,
556 material, lowEnergyLimit / massRatio);
559 G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
560 * lowEnergyLimit / scaledKineticEnergy);
563 material, genericIon,
564 scaledKineticEnergy, cutEnergy);
566 dEdx *= chargeSquare;
568 if(particle != genericIon) {
569 dEdx += corrections -> ComputeIonCorrections(particle,
570 material, kineticEnergy);
578 if (dEdx < 0.0) dEdx = 0.0;
593 G4double atomicMassNumber = particle -> GetAtomicMass();
594 G4double materialDensity = material -> GetDensity();
596 G4cout <<
"# dE/dx table for " << particle -> GetParticleName()
597 <<
" in material " << material ->
GetName()
598 <<
" of density " << materialDensity /
g *
cm3
601 <<
"# Projectile mass number A1 = " << atomicMassNumber
603 <<
"# ------------------------------------------------------"
607 << std::setw(14) <<
"E/A1"
608 << std::setw(14) <<
"dE/dx"
609 << std::setw(14) <<
"1/rho*dE/dx"
613 << std::setw(14) <<
"(MeV)"
614 << std::setw(14) <<
"(MeV/cm)"
615 << std::setw(14) <<
"(MeV*cm2/mg)"
617 <<
"# ------------------------------------------------------"
620 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
621 G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
625 energyLowerBoundary = std::log(energyLowerBoundary);
626 energyUpperBoundary = std::log(energyUpperBoundary);
629 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
632 for(
int i = 0; i < numBins + 1; i++) {
635 if(logScaleEnergy) energy = std::exp(energy);
640 << std::setw(14) << energy / atomicMassNumber /
MeV
641 << std::setw(14) << dedx /
MeV *
cm
642 << std::setw(14) << dedx / materialDensity / (
MeV*
cm2/(0.001*
g))
657 LossTableList::iterator iter = lossTableList.begin();
658 LossTableList::iterator iter_end = lossTableList.end();
660 for(;iter != iter_end; iter++) {
664 lowerBoundary, upperBoundary,
665 numBins,logScaleEnergy);
674 std::vector<G4DynamicParticle*>* secondaries,
708 std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
710 if(cutKinEnergySec >= maxKinEnergySec)
return;
712 G4double kineticEnergy = particle -> GetKineticEnergy();
716 G4double betaSquared = kineticEnergy *
717 (energy + cacheMass) / (energy * energy);
726 kinEnergySec = cutKinEnergySec * maxKinEnergySec /
727 (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
731 grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
734 G4cout <<
"G4IonParametrisedLossModel::SampleSecondary Warning: "
736 << grej <<
" for e= " << kinEnergySec
745 G4double totMomentum = energy*std::sqrt(betaSquared);
747 (momentumSec * totMomentum);
748 if(cost > 1.0) cost = 1.0;
749 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
753 G4ThreeVector directionSec(sint*std::cos(phi),sint*std::sin(phi), cost) ;
761 secondaries -> push_back(delta);
764 kineticEnergy -= kinEnergySec;
765 G4ThreeVector finalP = direction*totMomentum - directionSec*momentumSec;
766 finalP = finalP.
unit();
768 particleChangeLoss -> SetProposedKineticEnergy(kineticEnergy);
769 particleChangeLoss -> SetProposedMomentumDirection(finalP);
774 void G4IonParametrisedLossModel::UpdateRangeCache(
782 if(particle == rangeCacheParticle &&
783 matCutsCouple == rangeCacheMatCutsCouple) {
786 rangeCacheParticle = particle;
787 rangeCacheMatCutsCouple = matCutsCouple;
790 LossTableList::iterator iter =
IsApplicable(particle, material);
793 if(iter != lossTableList.end()) {
796 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
797 RangeEnergyTable::iterator iterRange = r.find(ionMatCouple);
799 if(iterRange == r.end()) BuildRangeVector(particle, matCutsCouple);
801 rangeCacheEnergyRange = E[ionMatCouple];
802 rangeCacheRangeEnergy = r[ionMatCouple];
805 rangeCacheEnergyRange = 0;
806 rangeCacheRangeEnergy = 0;
813 void G4IonParametrisedLossModel::UpdateDEDXCache(
826 if(particle == dedxCacheParticle &&
827 material == dedxCacheMaterial &&
828 cutEnergy == dedxCacheEnergyCut) {
832 dedxCacheParticle = particle;
834 dedxCacheEnergyCut = cutEnergy;
836 G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
837 dedxCacheGenIonMassRatio = massRatio;
839 LossTableList::iterator iter =
IsApplicable(particle, material);
840 dedxCacheIter = iter;
843 if(iter != lossTableList.end()) {
847 (*iter) -> GetUpperEnergyEdge(particle, material);
848 dedxCacheTransitionEnergy = transitionEnergy;
852 G4double dEdxParam = (*iter) -> GetDEDX(particle, material,
859 dEdxParam -= dEdxDeltaRays;
866 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
870 material, genericIon,
871 scaledTransitionEnergy, cutEnergy);
872 dEdxBetheBloch *= transitionChargeSquare;
876 corrections -> ComputeIonCorrections(particle,
877 material, transitionEnergy);
880 dedxCacheTransitionFactor =
881 (dEdxParam - dEdxBetheBloch)/dEdxBetheBloch
886 dedxCacheParticle = particle;
888 dedxCacheEnergyCut = cutEnergy;
890 dedxCacheGenIonMassRatio =
891 genericIonPDGMass / particle -> GetPDGMass();
893 dedxCacheTransitionEnergy = 0.0;
894 dedxCacheTransitionFactor = 0.0;
921 const G4Material* material = couple -> GetMaterial();
923 G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
925 if(kineticEnergy == eloss) {
return; }
928 size_t cutIndex = couple -> GetIndex();
929 cutEnergy = cutEnergies[cutIndex];
931 UpdateDEDXCache(particle, material, cutEnergy);
933 LossTableList::iterator iter = dedxCacheIter;
937 if(iter != lossTableList.end()) {
944 kineticEnergy, cutEnergy);
948 G4cout <<
"########################################################"
950 <<
"# G4IonParametrisedLossModel::CorrectionsAlongStep"
952 <<
"# cut(MeV) = " << cutEnergy/
MeV
957 << std::setw(14) <<
"l(um)"
958 << std::setw(14) <<
"l*dE/dx(MeV)"
959 << std::setw(14) <<
"(l*dE/dx)/E"
961 <<
"# ------------------------------------------------------"
965 << std::setw(14) << length / um
966 << std::setw(14) << eloss /
MeV
967 << std::setw(14) << eloss / kineticEnergy * 100.0
975 if(eloss > energyLossLimit * kineticEnergy) {
978 kineticEnergy,length);
981 G4cout <<
"# Correction applied:"
985 << std::setw(14) << length / um
986 << std::setw(14) << eloss /
MeV
987 << std::setw(14) << eloss / kineticEnergy * 100.0
997 if(energy < 0.0) energy = kineticEnergy * 0.5;
999 G4double chargeSquareRatio = corrections ->
1000 EffectiveChargeSquareRatio(particle,
1011 G4double transitionEnergy = dedxCacheTransitionEnergy;
1013 if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
1014 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1018 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1019 eloss *= chargeSquareRatioCorr;
1021 else if (iter == lossTableList.end()) {
1023 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
1027 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
1028 eloss *= chargeSquareRatioCorr;
1034 if(iter == lossTableList.end()) {
1036 G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
1040 if(scaledKineticEnergy > lowEnergyLimit)
1042 corrections -> IonHighOrderCorrections(particle, couple, energy);
1048 void G4IonParametrisedLossModel::BuildRangeVector(
1053 size_t cutIndex = matCutsCouple -> GetIndex();
1054 cutEnergy = cutEnergies[cutIndex];
1056 const G4Material* material = matCutsCouple -> GetMaterial();
1058 G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
1060 G4double lowerEnergy = lowerEnergyEdgeIntegr / massRatio;
1061 G4double upperEnergy = upperEnergyEdgeIntegr / massRatio;
1063 G4double logLowerEnergyEdge = std::log(lowerEnergy);
1064 G4double logUpperEnergyEdge = std::log(upperEnergy);
1066 G4double logDeltaEnergy = (logUpperEnergyEdge - logLowerEnergyEdge) /
1083 energyRangeVector -> PutValues(0, lowerEnergy, range);
1085 G4double logEnergy = std::log(lowerEnergy);
1086 for(
size_t i = 1; i < nmbBins+1; i++) {
1088 G4double logEnergyIntegr = logEnergy;
1090 for(
size_t j = 0; j < nmbSubBins; j++) {
1092 G4double binLowerBoundary = std::exp(logEnergyIntegr);
1093 logEnergyIntegr += logDeltaIntegr;
1095 G4double binUpperBoundary = std::exp(logEnergyIntegr);
1096 G4double deltaIntegr = binUpperBoundary - binLowerBoundary;
1098 G4double energyIntegr = binLowerBoundary + 0.5 * deltaIntegr;
1105 if(dedxValue > 0.0) range += deltaIntegr / dedxValue;
1107 #ifdef PRINT_DEBUG_DETAILS
1109 <<
" MeV -> dE = " << deltaIntegr/
MeV
1110 <<
" MeV -> dE/dx = " << dedxValue/
MeV*
mm
1111 <<
" MeV/mm -> dE/(dE/dx) = " << deltaIntegr /
1113 <<
" mm -> range = " << range /
mm
1118 logEnergy += logDeltaEnergy;
1122 energyRangeVector -> PutValues(i, energy, range);
1124 #ifdef PRINT_DEBUG_DETAILS
1125 G4cout <<
"G4IonParametrisedLossModel::BuildRangeVector() bin = "
1127 << energy /
MeV <<
" MeV, R = "
1128 << range /
mm <<
" mm"
1133 energyRangeVector -> SetSpline(
true);
1136 energyRangeVector ->
Value(lowerEnergy);
1138 energyRangeVector ->
Value(upperEnergy);
1145 for(
size_t i = 0; i < nmbBins+1; i++) {
1147 rangeEnergyVector ->
1148 PutValues(i, energyRangeVector ->
Value(energy), energy);
1151 rangeEnergyVector -> SetSpline(
true);
1153 #ifdef PRINT_DEBUG_TABLES
1154 G4cout << *energyLossVector
1155 << *energyRangeVector
1156 << *rangeEnergyVector <<
G4endl;
1159 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
1161 E[ionMatCouple] = energyRangeVector;
1162 r[ionMatCouple] = rangeEnergyVector;
1175 UpdateRangeCache(particle, matCutsCouple);
1180 if(energyRange != 0 && rangeEnergy != 0) {
1182 G4double lowerEnEdge = energyRange -> Energy( 0 );
1183 G4double lowerRangeEdge = rangeEnergy -> Energy( 0 );
1189 if(kineticEnergy < lowerEnEdge) {
1191 range = energyRange ->
Value(lowerEnEdge);
1192 range *= std::sqrt(kineticEnergy / lowerEnEdge);
1196 G4cout <<
"G4IonParametrisedLossModel::ComputeLossForStep() range = "
1197 << range /
mm <<
" mm, step = " << stepLength /
mm <<
" mm"
1202 G4double remRange = range - stepLength;
1206 if(remRange < 0.0) loss = kineticEnergy;
1207 else if(remRange < lowerRangeEdge) {
1209 G4double ratio = remRange / lowerRangeEdge;
1210 loss = kineticEnergy - ratio * ratio * lowerEnEdge;
1215 loss = kineticEnergy -
energy;
1219 if(loss < 0.0) loss = 0.0;
1232 G4cerr <<
"G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1233 <<
" add table: Invalid pointer."
1240 LossTableList::iterator iter = lossTableList.begin();
1241 LossTableList::iterator iter_end = lossTableList.end();
1243 for(;iter != iter_end; iter++) {
1246 if(tableName == nam) {
1247 G4cerr <<
"G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1248 <<
" add table: Name already exists."
1256 if(scalingAlgorithm == 0)
1262 lossTableList.push_front(handler);
1272 LossTableList::iterator iter = lossTableList.begin();
1273 LossTableList::iterator iter_end = lossTableList.end();
1275 for(;iter != iter_end; iter++) {
1278 if(tableName == nam) {
1282 lossTableList.erase(iter);
1285 RangeEnergyTable::iterator iterRange = r.begin();
1286 RangeEnergyTable::iterator iterRange_end = r.end();
1288 for(;iterRange != iterRange_end; iterRange++)
1289 delete iterRange ->
second;
1292 EnergyRangeTable::iterator iterEnergy = E.begin();
1293 EnergyRangeTable::iterator iterEnergy_end = E.end();
1295 for(;iterEnergy != iterEnergy_end; iterEnergy++)
1296 delete iterEnergy ->
second;
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
G4double LowEnergyLimit() const
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
G4double HighEnergyLimit() const
static G4GenericIon * Definition()
static G4Electron * Definition()
G4bool AddDEDXTable(const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
G4VEmFluctuationModel * GetModelOfFluctuations()
virtual ~G4IonParametrisedLossModel()
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
void SetHighEnergyLimit(G4double)
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double)
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition const G4Material *G4double range
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=0)
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double)
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
G4double DeltaRayMeanEnergyTransferRate(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
void PrintDEDXTableHandlers(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
void DeactivateICRU73Scaling()
static G4ProductionCutsTable * GetProductionCutsTable()
G4bool RemoveDEDXTable(const G4String &name)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double)
std::pair< const G4ParticleDefinition *, const G4MaterialCutsCouple * > IonMatCouple
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double, G4double)
const G4String & GetName() const
G4IonParametrisedLossModel(const G4ParticleDefinition *particle=0, const G4String &name="ParamICRU73")
LossTableList::iterator IsApplicable(const G4ParticleDefinition *, const G4Material *)
void PrintDEDXTable(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
G4double ComputeLossForStep(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double, G4double)
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
G4GLOB_DLL std::ostream G4cerr
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &, G4double &, G4double)