112 particleChangeLoss(0),
114 energyLossLimit(0.01),
168 LossTableList::iterator iterTables_end =
lossTableList.end();
170 for(;iterTables != iterTables_end; ++iterTables) {
delete *iterTables; }
174 RangeEnergyTable::iterator itr =
r.begin();
175 RangeEnergyTable::iterator itr_end =
r.end();
176 for(;itr != itr_end; ++itr) {
delete itr->second; }
180 EnergyRangeTable::iterator ite =
E.begin();
181 EnergyRangeTable::iterator ite_end =
E.end();
182 for(;ite != ite_end; ++ite) {
delete ite->second; }
193 return couple -> GetMaterial() -> GetIonisation() ->
194 GetMeanExcitationEnergy();
218 G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
233 EffectiveChargeSquareRatio(particle,
282 LossTableList::iterator iterTables_end =
lossTableList.end();
284 for(;iterTables != iterTables_end; iterTables++)
285 (*iterTables) -> ClearCache();
289 RangeEnergyTable::iterator iterRange =
r.begin();
290 RangeEnergyTable::iterator iterRange_end =
r.end();
292 for(;iterRange != iterRange_end; iterRange++) {
293 delete iterRange->second;
297 EnergyRangeTable::iterator iterEnergy =
E.begin();
298 EnergyRangeTable::iterator iterEnergy_end =
E.end();
300 for(;iterEnergy != iterEnergy_end; iterEnergy++) {
301 delete iterEnergy->second;
306 size_t size = cuts.size();
308 for(
size_t i = 0; i < size; i++)
cutEnergies.push_back(cuts[i]);
313 size_t nmbCouples = coupleTable -> GetTableSize();
315 #ifdef PRINT_TABLE_BUILT
316 G4cout <<
"G4IonParametrisedLossModel::Initialise():"
317 <<
" Building dE/dx vectors:"
321 for (
size_t i = 0; i < nmbCouples; i++) {
324 coupleTable -> GetMaterialCutsCouple(i);
326 const G4Material* material = couple -> GetMaterial();
329 for(
G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) {
334 for(;iter != iter_end; iter++) {
337 G4cout <<
"G4IonParametrisedLossModel::Initialise():"
338 <<
" Skipping illegal table."
343 (*iter) -> BuildDEDXTable(atomicNumberIon, material);
346 #ifdef PRINT_TABLE_BUILT
347 G4cout <<
" Atomic Number Ion = " << atomicNumberIon
348 <<
", Material = " << material ->
GetName()
349 <<
", Table = " << (*iter) ->
GetName()
395 if(cutEnergy < tmax) {
398 G4double betaSquared = kineticEnergy *
399 (energy +
cacheMass) / (energy * energy);
401 crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
402 betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
407 #ifdef PRINT_DEBUG_CS
408 G4cout <<
"########################################################"
410 <<
"# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
412 <<
"# particle =" << particle -> GetParticleName()
414 <<
"# cut(MeV) = " << cutEnergy/
MeV
419 << std::setw(14) <<
"CS(um)"
420 << std::setw(14) <<
"E_max_sec(MeV)"
422 <<
"# ------------------------------------------------------"
426 << std::setw(14) << crosssection / (um * um)
427 << std::setw(14) << tmax /
MeV
431 crosssection *= atomicNumber;
445 G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
499 if(transitionEnergy > kineticEnergy) {
501 dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
507 dEdx -= dEdxDeltaRays;
515 G4double scaledKineticEnergy = kineticEnergy * massRatio;
516 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
520 if(scaledTransitionEnergy >= lowEnergyLimit) {
524 scaledKineticEnergy, cutEnergy);
526 dEdx *= chargeSquare;
528 dEdx +=
corrections -> ComputeIonCorrections(particle,
529 material, kineticEnergy);
548 G4double scaledKineticEnergy = kineticEnergy * massRatio;
551 if(scaledKineticEnergy < lowEnergyLimit) {
554 scaledKineticEnergy, cutEnergy);
556 dEdx *= chargeSquare;
561 lowEnergyLimit, cutEnergy);
565 lowEnergyLimit, cutEnergy);
568 G4double chargeSquareLowEnergyLimit =
570 lowEnergyLimit / massRatio);
572 dEdxLimitParam *= chargeSquareLowEnergyLimit;
573 dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
575 dEdxLimitBetheBloch +=
577 material, lowEnergyLimit / massRatio);
580 G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
581 * lowEnergyLimit / scaledKineticEnergy);
585 scaledKineticEnergy, cutEnergy);
587 dEdx *= chargeSquare;
590 dEdx +=
corrections -> ComputeIonCorrections(particle,
591 material, kineticEnergy);
599 if (dEdx < 0.0) dEdx = 0.0;
614 G4double atomicMassNumber = particle -> GetAtomicMass();
615 G4double materialDensity = material -> GetDensity();
617 G4cout <<
"# dE/dx table for " << particle -> GetParticleName()
618 <<
" in material " << material ->
GetName()
619 <<
" of density " << materialDensity /
g *
cm3
622 <<
"# Projectile mass number A1 = " << atomicMassNumber
624 <<
"# ------------------------------------------------------"
628 << std::setw(14) <<
"E/A1"
629 << std::setw(14) <<
"dE/dx"
630 << std::setw(14) <<
"1/rho*dE/dx"
634 << std::setw(14) <<
"(MeV)"
635 << std::setw(14) <<
"(MeV/cm)"
636 << std::setw(14) <<
"(MeV*cm2/mg)"
638 <<
"# ------------------------------------------------------"
641 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
642 G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
646 energyLowerBoundary = std::log(energyLowerBoundary);
647 energyUpperBoundary = std::log(energyUpperBoundary);
650 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
653 for(
int i = 0; i < numBins + 1; i++) {
656 if(logScaleEnergy) energy = std::exp(energy);
661 << std::setw(14) << energy / atomicMassNumber /
MeV
662 << std::setw(14) << dedx /
MeV *
cm
663 << std::setw(14) << dedx / materialDensity / (
MeV*
cm2/(0.001*
g))
681 for(;iter != iter_end; iter++) {
685 lowerBoundary, upperBoundary,
686 numBins,logScaleEnergy);
695 std::vector<G4DynamicParticle*>* secondaries,
729 std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
731 if(cutKinEnergySec >= maxKinEnergySec)
return;
733 G4double kineticEnergy = particle -> GetKineticEnergy();
746 kinEnergySec = cutKinEnergySec * maxKinEnergySec /
747 (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
751 grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
754 G4cout <<
"G4IonParametrisedLossModel::SampleSecondary Warning: "
756 << grej <<
" for e= " << kinEnergySec
773 secondaries->push_back(delta);
777 G4double totalMomentum = std::sqrt(kineticEnergy*(energy + cacheMass));
780 finalP = finalP.unit();
782 kineticEnergy -= kinEnergySec;
805 const G4Material* material = matCutsCouple -> GetMaterial();
806 LossTableList::iterator iter =
IsApplicable(particle, material);
812 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
813 RangeEnergyTable::iterator iterRange =
r.find(ionMatCouple);
855 LossTableList::iterator iter =
IsApplicable(particle, material);
863 (*iter) -> GetUpperEnergyEdge(particle, material);
868 G4double dEdxParam = (*iter) -> GetDEDX(particle, material,
875 dEdxParam -= dEdxDeltaRays;
882 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
887 scaledTransitionEnergy, cutEnergy);
888 dEdxBetheBloch *= transitionChargeSquare;
893 material, transitionEnergy);
897 (dEdxParam - dEdxBetheBloch)/dEdxBetheBloch
937 const G4Material* material = couple -> GetMaterial();
939 G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
941 if(kineticEnergy == eloss) {
return; }
944 size_t cutIndex = couple -> GetIndex();
960 kineticEnergy, cutEnergy);
964 G4cout <<
"########################################################"
966 <<
"# G4IonParametrisedLossModel::CorrectionsAlongStep"
968 <<
"# cut(MeV) = " << cutEnergy/
MeV
973 << std::setw(14) <<
"l(um)"
974 << std::setw(14) <<
"l*dE/dx(MeV)"
975 << std::setw(14) <<
"(l*dE/dx)/E"
977 <<
"# ------------------------------------------------------"
981 << std::setw(14) << length / um
982 << std::setw(14) << eloss /
MeV
983 << std::setw(14) << eloss / kineticEnergy * 100.0
994 kineticEnergy,length);
997 G4cout <<
"# Correction applied:"
1001 << std::setw(14) << length / um
1002 << std::setw(14) << eloss /
MeV
1003 << std::setw(14) << eloss / kineticEnergy * 100.0
1013 if(energy < 0.0) energy = kineticEnergy * 0.5;
1016 EffectiveChargeSquareRatio(particle,
1029 if(iter !=
lossTableList.end() && transitionEnergy < kineticEnergy) {
1030 chargeSquareRatio *=
corrections -> EffectiveChargeCorrection(particle,
1035 eloss *= chargeSquareRatioCorr;
1039 chargeSquareRatio *=
corrections -> EffectiveChargeCorrection(particle,
1044 eloss *= chargeSquareRatioCorr;
1056 if(scaledKineticEnergy > lowEnergyLimit)
1058 corrections -> IonHighOrderCorrections(particle, couple, energy);
1069 size_t cutIndex = matCutsCouple -> GetIndex();
1072 const G4Material* material = matCutsCouple -> GetMaterial();
1079 G4double logLowerEnergyEdge = std::log(lowerEnergy);
1080 G4double logUpperEnergyEdge = std::log(upperEnergy);
1082 G4double logDeltaEnergy = (logUpperEnergyEdge - logLowerEnergyEdge) /
1097 G4double range = 2.0 * lowerEnergy / dedxLow;
1099 energyRangeVector -> PutValues(0, lowerEnergy, range);
1101 G4double logEnergy = std::log(lowerEnergy);
1102 for(
size_t i = 1; i <
nmbBins+1; i++) {
1104 G4double logEnergyIntegr = logEnergy;
1108 G4double binLowerBoundary = std::exp(logEnergyIntegr);
1109 logEnergyIntegr += logDeltaIntegr;
1111 G4double binUpperBoundary = std::exp(logEnergyIntegr);
1112 G4double deltaIntegr = binUpperBoundary - binLowerBoundary;
1114 G4double energyIntegr = binLowerBoundary + 0.5 * deltaIntegr;
1121 if(dedxValue > 0.0) range += deltaIntegr / dedxValue;
1123 #ifdef PRINT_DEBUG_DETAILS
1125 <<
" MeV -> dE = " << deltaIntegr/
MeV
1126 <<
" MeV -> dE/dx = " << dedxValue/
MeV*
mm
1127 <<
" MeV/mm -> dE/(dE/dx) = " << deltaIntegr /
1129 <<
" mm -> range = " << range /
mm
1134 logEnergy += logDeltaEnergy;
1138 energyRangeVector -> PutValues(i, energy, range);
1140 #ifdef PRINT_DEBUG_DETAILS
1141 G4cout <<
"G4IonParametrisedLossModel::BuildRangeVector() bin = "
1143 << energy /
MeV <<
" MeV, R = "
1144 << range /
mm <<
" mm"
1149 energyRangeVector -> SetSpline(
true);
1152 energyRangeVector ->
Value(lowerEnergy);
1154 energyRangeVector ->
Value(upperEnergy);
1161 for(
size_t i = 0; i < nmbBins+1; i++) {
1163 rangeEnergyVector ->
1164 PutValues(i, energyRangeVector ->
Value(energy), energy);
1167 rangeEnergyVector -> SetSpline(
true);
1169 #ifdef PRINT_DEBUG_TABLES
1170 G4cout << *energyLossVector
1171 << *energyRangeVector
1172 << *rangeEnergyVector <<
G4endl;
1175 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
1177 E[ionMatCouple] = energyRangeVector;
1178 r[ionMatCouple] = rangeEnergyVector;
1196 if(energyRange != 0 && rangeEnergy != 0) {
1198 G4double lowerEnEdge = energyRange -> Energy( 0 );
1199 G4double lowerRangeEdge = rangeEnergy -> Energy( 0 );
1205 if(kineticEnergy < lowerEnEdge) {
1207 range = energyRange ->
Value(lowerEnEdge);
1208 range *= std::sqrt(kineticEnergy / lowerEnEdge);
1212 G4cout <<
"G4IonParametrisedLossModel::ComputeLossForStep() range = "
1213 << range /
mm <<
" mm, step = " << stepLength /
mm <<
" mm"
1218 G4double remRange = range - stepLength;
1222 if(remRange < 0.0) loss = kineticEnergy;
1223 else if(remRange < lowerRangeEdge) {
1225 G4double ratio = remRange / lowerRangeEdge;
1226 loss = kineticEnergy - ratio * ratio * lowerEnEdge;
1231 loss = kineticEnergy -
energy;
1235 if(loss < 0.0) loss = 0.0;
1248 G4cerr <<
"G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1249 <<
" add table: Invalid pointer."
1259 for(;iter != iter_end; iter++) {
1262 if(tableName == nam) {
1263 G4cerr <<
"G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1264 <<
" add table: Name already exists."
1272 if(scalingAlgorithm == 0)
1291 for(;iter != iter_end; iter++) {
1294 if(tableName == nam) {
1301 RangeEnergyTable::iterator iterRange =
r.begin();
1302 RangeEnergyTable::iterator iterRange_end =
r.end();
1304 for(;iterRange != iterRange_end; iterRange++)
1305 delete iterRange ->
second;
1308 EnergyRangeTable::iterator iterEnergy =
E.begin();
1309 EnergyRangeTable::iterator iterEnergy_end =
E.end();
1311 for(;iterEnergy != iterEnergy_end; iterEnergy++)
1312 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
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)