113 particleChangeLoss(0),
115 energyLossLimit(0.01),
169 LossTableList::iterator iterTables_end =
lossTableList.end();
171 for(;iterTables != iterTables_end; ++iterTables) {
delete *iterTables; }
175 RangeEnergyTable::iterator itr =
r.begin();
176 RangeEnergyTable::iterator itr_end =
r.end();
177 for(;itr != itr_end; ++itr) {
delete itr->second; }
181 EnergyRangeTable::iterator ite =
E.begin();
182 EnergyRangeTable::iterator ite_end =
E.end();
183 for(;ite != ite_end; ++ite) {
delete ite->second; }
194 return couple -> GetMaterial() -> GetIonisation() ->
195 GetMeanExcitationEnergy();
219 G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
234 EffectiveChargeSquareRatio(particle,
283 LossTableList::iterator iterTables_end =
lossTableList.end();
285 for(;iterTables != iterTables_end; iterTables++)
286 (*iterTables) -> ClearCache();
290 RangeEnergyTable::iterator iterRange =
r.begin();
291 RangeEnergyTable::iterator iterRange_end =
r.end();
293 for(;iterRange != iterRange_end; iterRange++) {
294 delete iterRange->second;
298 EnergyRangeTable::iterator iterEnergy =
E.begin();
299 EnergyRangeTable::iterator iterEnergy_end =
E.end();
301 for(;iterEnergy != iterEnergy_end; iterEnergy++) {
302 delete iterEnergy->second;
307 size_t size = cuts.size();
309 for(
size_t i = 0; i < size; i++)
cutEnergies.push_back(cuts[i]);
314 size_t nmbCouples = coupleTable -> GetTableSize();
316 #ifdef PRINT_TABLE_BUILT
317 G4cout <<
"G4IonParametrisedLossModel::Initialise():"
318 <<
" Building dE/dx vectors:"
322 for (
size_t i = 0; i < nmbCouples; i++) {
325 coupleTable -> GetMaterialCutsCouple(i);
327 const G4Material* material = couple -> GetMaterial();
330 for(
G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) {
335 for(;iter != iter_end; iter++) {
338 G4cout <<
"G4IonParametrisedLossModel::Initialise():"
339 <<
" Skipping illegal table."
344 (*iter) -> BuildDEDXTable(atomicNumberIon, material);
347 #ifdef PRINT_TABLE_BUILT
348 G4cout <<
" Atomic Number Ion = " << atomicNumberIon
349 <<
", Material = " << material ->
GetName()
350 <<
", Table = " << (*iter) ->
GetName()
396 if(cutEnergy < tmax) {
399 G4double betaSquared = kineticEnergy *
400 (energy +
cacheMass) / (energy * energy);
402 crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
403 betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
408 #ifdef PRINT_DEBUG_CS
409 G4cout <<
"########################################################"
411 <<
"# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
413 <<
"# particle =" << particle -> GetParticleName()
415 <<
"# cut(MeV) = " << cutEnergy/
MeV
420 << std::setw(14) <<
"CS(um)"
421 << std::setw(14) <<
"E_max_sec(MeV)"
423 <<
"# ------------------------------------------------------"
427 << std::setw(14) << crosssection / (
um *
um)
428 << std::setw(14) << tmax /
MeV
432 crosssection *= atomicNumber;
446 G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
500 if(transitionEnergy > kineticEnergy) {
502 dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
508 dEdx -= dEdxDeltaRays;
516 G4double scaledKineticEnergy = kineticEnergy * massRatio;
517 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
521 if(scaledTransitionEnergy >= lowEnergyLimit) {
525 scaledKineticEnergy, cutEnergy);
527 dEdx *= chargeSquare;
529 dEdx +=
corrections -> ComputeIonCorrections(particle,
530 material, kineticEnergy);
549 G4double scaledKineticEnergy = kineticEnergy * massRatio;
552 if(scaledKineticEnergy < lowEnergyLimit) {
555 scaledKineticEnergy, cutEnergy);
557 dEdx *= chargeSquare;
562 lowEnergyLimit, cutEnergy);
566 lowEnergyLimit, cutEnergy);
569 G4double chargeSquareLowEnergyLimit =
571 lowEnergyLimit / massRatio);
573 dEdxLimitParam *= chargeSquareLowEnergyLimit;
574 dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
576 dEdxLimitBetheBloch +=
578 material, lowEnergyLimit / massRatio);
581 G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
582 * lowEnergyLimit / scaledKineticEnergy);
586 scaledKineticEnergy, cutEnergy);
588 dEdx *= chargeSquare;
591 dEdx +=
corrections -> ComputeIonCorrections(particle,
592 material, kineticEnergy);
600 if (dEdx < 0.0) dEdx = 0.0;
615 G4double atomicMassNumber = particle -> GetAtomicMass();
616 G4double materialDensity = material -> GetDensity();
618 G4cout <<
"# dE/dx table for " << particle -> GetParticleName()
619 <<
" in material " << material ->
GetName()
620 <<
" of density " << materialDensity /
g *
cm3
623 <<
"# Projectile mass number A1 = " << atomicMassNumber
625 <<
"# ------------------------------------------------------"
629 << std::setw(14) <<
"E/A1"
630 << std::setw(14) <<
"dE/dx"
631 << std::setw(14) <<
"1/rho*dE/dx"
635 << std::setw(14) <<
"(MeV)"
636 << std::setw(14) <<
"(MeV/cm)"
637 << std::setw(14) <<
"(MeV*cm2/mg)"
639 <<
"# ------------------------------------------------------"
642 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
643 G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
647 energyLowerBoundary = std::log(energyLowerBoundary);
648 energyUpperBoundary = std::log(energyUpperBoundary);
651 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
654 for(
int i = 0; i < numBins + 1; i++) {
657 if(logScaleEnergy) energy =
G4Exp(energy);
662 << std::setw(14) << energy / atomicMassNumber /
MeV
663 << std::setw(14) << dedx /
MeV *
cm
664 << std::setw(14) << dedx / materialDensity / (
MeV*
cm2/(0.001*
g))
682 for(;iter != iter_end; iter++) {
686 lowerBoundary, upperBoundary,
687 numBins,logScaleEnergy);
696 std::vector<G4DynamicParticle*>* secondaries,
730 std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
732 if(cutKinEnergySec >= maxKinEnergySec)
return;
734 G4double kineticEnergy = particle -> GetKineticEnergy();
747 kinEnergySec = cutKinEnergySec * maxKinEnergySec /
748 (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
752 grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
755 G4cout <<
"G4IonParametrisedLossModel::SampleSecondary Warning: "
757 << grej <<
" for e= " << kinEnergySec
774 secondaries->push_back(delta);
778 G4double totalMomentum = std::sqrt(kineticEnergy*(energy + cacheMass));
781 finalP = finalP.unit();
783 kineticEnergy -= kinEnergySec;
806 const G4Material* material = matCutsCouple -> GetMaterial();
807 LossTableList::iterator iter =
IsApplicable(particle, material);
813 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
814 RangeEnergyTable::iterator iterRange =
r.find(ionMatCouple);
856 LossTableList::iterator iter =
IsApplicable(particle, material);
864 (*iter) -> GetUpperEnergyEdge(particle, material);
869 G4double dEdxParam = (*iter) -> GetDEDX(particle, material,
876 dEdxParam -= dEdxDeltaRays;
883 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
888 scaledTransitionEnergy, cutEnergy);
889 dEdxBetheBloch *= transitionChargeSquare;
894 material, transitionEnergy);
898 (dEdxParam - dEdxBetheBloch)/dEdxBetheBloch
938 const G4Material* material = couple -> GetMaterial();
940 G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
942 if(kineticEnergy == eloss) {
return; }
945 size_t cutIndex = couple -> GetIndex();
961 kineticEnergy, cutEnergy);
965 G4cout <<
"########################################################"
967 <<
"# G4IonParametrisedLossModel::CorrectionsAlongStep"
969 <<
"# cut(MeV) = " << cutEnergy/
MeV
974 << std::setw(14) <<
"l(um)"
975 << std::setw(14) <<
"l*dE/dx(MeV)"
976 << std::setw(14) <<
"(l*dE/dx)/E"
978 <<
"# ------------------------------------------------------"
982 << std::setw(14) << length /
um
983 << std::setw(14) << eloss /
MeV
984 << std::setw(14) << eloss / kineticEnergy * 100.0
995 kineticEnergy,length);
998 G4cout <<
"# Correction applied:"
1002 << std::setw(14) << length /
um
1003 << std::setw(14) << eloss /
MeV
1004 << std::setw(14) << eloss / kineticEnergy * 100.0
1014 if(energy < 0.0) energy = kineticEnergy * 0.5;
1017 EffectiveChargeSquareRatio(particle,
1030 if(iter !=
lossTableList.end() && transitionEnergy < kineticEnergy) {
1031 chargeSquareRatio *=
corrections -> EffectiveChargeCorrection(particle,
1036 eloss *= chargeSquareRatioCorr;
1040 chargeSquareRatio *=
corrections -> EffectiveChargeCorrection(particle,
1045 eloss *= chargeSquareRatioCorr;
1057 if(scaledKineticEnergy > lowEnergyLimit)
1059 corrections -> IonHighOrderCorrections(particle, couple, energy);
1070 size_t cutIndex = matCutsCouple -> GetIndex();
1073 const G4Material* material = matCutsCouple -> GetMaterial();
1080 G4double logLowerEnergyEdge = std::log(lowerEnergy);
1081 G4double logUpperEnergyEdge = std::log(upperEnergy);
1083 G4double logDeltaEnergy = (logUpperEnergyEdge - logLowerEnergyEdge) /
1098 G4double range = 2.0 * lowerEnergy / dedxLow;
1100 energyRangeVector -> PutValues(0, lowerEnergy, range);
1102 G4double logEnergy = std::log(lowerEnergy);
1103 for(
size_t i = 1; i <
nmbBins+1; i++) {
1105 G4double logEnergyIntegr = logEnergy;
1110 logEnergyIntegr += logDeltaIntegr;
1113 G4double deltaIntegr = binUpperBoundary - binLowerBoundary;
1115 G4double energyIntegr = binLowerBoundary + 0.5 * deltaIntegr;
1122 if(dedxValue > 0.0) range += deltaIntegr / dedxValue;
1124 #ifdef PRINT_DEBUG_DETAILS
1126 <<
" MeV -> dE = " << deltaIntegr/
MeV
1127 <<
" MeV -> dE/dx = " << dedxValue/
MeV*
mm
1128 <<
" MeV/mm -> dE/(dE/dx) = " << deltaIntegr /
1130 <<
" mm -> range = " << range /
mm
1135 logEnergy += logDeltaEnergy;
1139 energyRangeVector -> PutValues(i, energy, range);
1141 #ifdef PRINT_DEBUG_DETAILS
1142 G4cout <<
"G4IonParametrisedLossModel::BuildRangeVector() bin = "
1144 << energy /
MeV <<
" MeV, R = "
1145 << range /
mm <<
" mm"
1150 energyRangeVector -> SetSpline(
true);
1153 energyRangeVector ->
Value(lowerEnergy);
1155 energyRangeVector ->
Value(upperEnergy);
1162 for(
size_t i = 0; i < nmbBins+1; i++) {
1164 rangeEnergyVector ->
1165 PutValues(i, energyRangeVector ->
Value(energy), energy);
1168 rangeEnergyVector -> SetSpline(
true);
1170 #ifdef PRINT_DEBUG_TABLES
1171 G4cout << *energyLossVector
1172 << *energyRangeVector
1173 << *rangeEnergyVector <<
G4endl;
1176 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
1178 E[ionMatCouple] = energyRangeVector;
1179 r[ionMatCouple] = rangeEnergyVector;
1197 if(energyRange != 0 && rangeEnergy != 0) {
1199 G4double lowerEnEdge = energyRange -> Energy( 0 );
1200 G4double lowerRangeEdge = rangeEnergy -> Energy( 0 );
1206 if(kineticEnergy < lowerEnEdge) {
1208 range = energyRange ->
Value(lowerEnEdge);
1209 range *= std::sqrt(kineticEnergy / lowerEnEdge);
1213 G4cout <<
"G4IonParametrisedLossModel::ComputeLossForStep() range = "
1214 << range /
mm <<
" mm, step = " << stepLength /
mm <<
" mm"
1219 G4double remRange = range - stepLength;
1223 if(remRange < 0.0) loss = kineticEnergy;
1224 else if(remRange < lowerRangeEdge) {
1226 G4double ratio = remRange / lowerRangeEdge;
1227 loss = kineticEnergy - ratio * ratio * lowerEnEdge;
1232 loss = kineticEnergy -
energy;
1236 if(loss < 0.0) loss = 0.0;
1249 G4cerr <<
"G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1250 <<
" add table: Invalid pointer."
1260 for(;iter != iter_end; iter++) {
1263 if(tableName == nam) {
1264 G4cerr <<
"G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1265 <<
" add table: Name already exists."
1273 if(scalingAlgorithm == 0)
1292 for(;iter != iter_end; iter++) {
1295 if(tableName == nam) {
1302 RangeEnergyTable::iterator iterRange =
r.begin();
1303 RangeEnergyTable::iterator iterRange_end =
r.end();
1305 for(;iterRange != iterRange_end; iterRange++)
1306 delete iterRange ->
second;
1309 EnergyRangeTable::iterator iterEnergy =
E.begin();
1310 EnergyRangeTable::iterator iterEnergy_end =
E.end();
1312 for(;iterEnergy != iterEnergy_end; iterEnergy++)
1313 delete iterEnergy ->
second;
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
G4double LowEnergyLimit() const
G4double dedxCacheTransitionEnergy
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
G4double upperEnergyEdgeIntegr
void BuildRangeVector(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
static G4GenericIon * Definition()
G4VEmAngularDistribution * GetAngularDistribution()
G4EmCorrections * corrections
G4bool AddDEDXTable(const G4String &name, G4VIonDEDXTable *table, G4VIonDEDXScalingAlgorithm *algorithm=0)
G4BetheBlochModel * betheBlochModel
G4VEmFluctuationModel * GetModelOfFluctuations()
virtual ~G4IonParametrisedLossModel()
G4double dedxCacheTransitionFactor
G4double cacheElecMassRatio
const G4ParticleDefinition * cacheParticle
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
void UpdateRangeCache(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
void SetHighEnergyLimit(G4double)
const G4MaterialCutsCouple * rangeCacheMatCutsCouple
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
LossTableList::iterator dedxCacheIter
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double)
G4GLOB_DLL std::ostream G4cout
G4ParticleChangeForLoss * particleChangeLoss
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=0)
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
G4double dedxCacheGenIonMassRatio
const G4Material * dedxCacheMaterial
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double)
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static const double second
G4double cacheChargeSquare
LossTableList lossTableList
void SetProposedMomentumDirection(const G4ThreeVector &dir)
G4double DeltaRayMeanEnergyTransferRate(const G4Material *, const G4ParticleDefinition *, G4double, G4double)
G4double genericIonPDGMass
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4PhysicsVector * rangeCacheRangeEnergy
void PrintDEDXTableHandlers(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
void DeactivateICRU73Scaling()
static G4ProductionCutsTable * GetProductionCutsTable()
G4bool RemoveDEDXTable(const G4String &name)
void UpdateCache(const G4ParticleDefinition *)
G4double lowerEnergyEdgeIntegr
static const G4double factor
G4PhysicsVector * rangeCacheEnergyRange
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double)
G4double energy(const ThreeVector &p, const G4double m)
std::pair< const G4ParticleDefinition *, const G4MaterialCutsCouple * > IonMatCouple
void SetAngularDistribution(G4VEmAngularDistribution *)
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)
static G4Electron * Electron()
const G4String & GetName() const
G4double dedxCacheEnergyCut
G4IonParametrisedLossModel(const G4ParticleDefinition *particle=0, const G4String &name="ParamICRU73")
LossTableList::iterator IsApplicable(const G4ParticleDefinition *, const G4Material *)
void UpdateDEDXCache(const G4ParticleDefinition *, const G4Material *, G4double cutEnergy)
const G4ParticleDefinition * dedxCacheParticle
void PrintDEDXTable(const G4ParticleDefinition *, const G4Material *, G4double, G4double, G4int, G4bool)
G4double ComputeLossForStep(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double, G4double)
const G4ParticleDefinition * rangeCacheParticle
G4BraggIonModel * braggIonModel
G4ThreeVector GetMomentum() const
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
const G4Material * GetMaterial() const
G4GLOB_DLL std::ostream G4cerr
G4int SelectRandomAtomNumber(const G4Material *)
G4ParticleDefinition * genericIon
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &, G4double &, G4double)