00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 #include "G4IonParametrisedLossModel.hh"
00082 #include "G4PhysicalConstants.hh"
00083 #include "G4SystemOfUnits.hh"
00084 #include "G4LPhysicsFreeVector.hh"
00085 #include "G4IonStoppingData.hh"
00086 #include "G4VIonDEDXTable.hh"
00087 #include "G4VIonDEDXScalingAlgorithm.hh"
00088 #include "G4IonDEDXScalingICRU73.hh"
00089 #include "G4BraggIonModel.hh"
00090 #include "G4BetheBlochModel.hh"
00091 #include "G4ProductionCutsTable.hh"
00092 #include "G4ParticleChangeForLoss.hh"
00093 #include "G4LossTableManager.hh"
00094 #include "G4GenericIon.hh"
00095 #include "G4Electron.hh"
00096 #include "Randomize.hh"
00097
00098
00099
00100
00101
00102
00103 G4IonParametrisedLossModel::G4IonParametrisedLossModel(
00104 const G4ParticleDefinition*,
00105 const G4String& nam)
00106 : G4VEmModel(nam),
00107 braggIonModel(0),
00108 betheBlochModel(0),
00109 nmbBins(90),
00110 nmbSubBins(100),
00111 particleChangeLoss(0),
00112 corrFactor(1.0),
00113 energyLossLimit(0.01),
00114 cutEnergies(0)
00115 {
00116 genericIon = G4GenericIon::Definition();
00117 genericIonPDGMass = genericIon -> GetPDGMass();
00118 corrections = G4LossTableManager::Instance() -> EmCorrections();
00119
00120
00121 SetHighEnergyLimit(100.0 * TeV);
00122
00123
00124 braggIonModel = new G4BraggIonModel();
00125 betheBlochModel = new G4BetheBlochModel();
00126
00127
00128 AddDEDXTable("ICRU73",
00129 new G4IonStoppingData("ion_stopping_data/icru73"),
00130 new G4IonDEDXScalingICRU73());
00131
00132
00133 lowerEnergyEdgeIntegr = 0.025 * MeV;
00134 upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit();
00135
00136
00137 cacheParticle = 0;
00138 cacheMass = 0;
00139 cacheElecMassRatio = 0;
00140 cacheChargeSquare = 0;
00141
00142
00143 rangeCacheParticle = 0;
00144 rangeCacheMatCutsCouple = 0;
00145 rangeCacheEnergyRange = 0;
00146 rangeCacheRangeEnergy = 0;
00147
00148
00149 dedxCacheParticle = 0;
00150 dedxCacheMaterial = 0;
00151 dedxCacheEnergyCut = 0;
00152 dedxCacheIter = lossTableList.end();
00153 dedxCacheTransitionEnergy = 0.0;
00154 dedxCacheTransitionFactor = 0.0;
00155 dedxCacheGenIonMassRatio = 0.0;
00156 }
00157
00158
00159
00160 G4IonParametrisedLossModel::~G4IonParametrisedLossModel() {
00161
00162
00163 RangeEnergyTable::iterator iterRange = r.begin();
00164 RangeEnergyTable::iterator iterRange_end = r.end();
00165
00166 for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second;
00167 r.clear();
00168
00169
00170 EnergyRangeTable::iterator iterEnergy = E.begin();
00171 EnergyRangeTable::iterator iterEnergy_end = E.end();
00172
00173 for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second;
00174 E.clear();
00175
00176
00177 LossTableList::iterator iterTables = lossTableList.begin();
00178 LossTableList::iterator iterTables_end = lossTableList.end();
00179
00180 for(;iterTables != iterTables_end; iterTables++) delete *iterTables;
00181 lossTableList.clear();
00182
00183
00184 delete betheBlochModel;
00185 delete braggIonModel;
00186 }
00187
00188
00189
00190 G4double G4IonParametrisedLossModel::MinEnergyCut(
00191 const G4ParticleDefinition*,
00192 const G4MaterialCutsCouple* couple) {
00193
00194 return couple -> GetMaterial() -> GetIonisation() ->
00195 GetMeanExcitationEnergy();
00196 }
00197
00198
00199
00200 G4double G4IonParametrisedLossModel::MaxSecondaryEnergy(
00201 const G4ParticleDefinition* particle,
00202 G4double kineticEnergy) {
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 if(particle != cacheParticle) UpdateCache(particle);
00217
00218 G4double tau = kineticEnergy/cacheMass;
00219 G4double tmax = 2.0 * electron_mass_c2 * tau * (tau + 2.) /
00220 (1. + 2.0 * (tau + 1.) * cacheElecMassRatio +
00221 cacheElecMassRatio * cacheElecMassRatio);
00222
00223 return tmax;
00224 }
00225
00226
00227
00228 G4double G4IonParametrisedLossModel::GetChargeSquareRatio(
00229 const G4ParticleDefinition* particle,
00230 const G4Material* material,
00231 G4double kineticEnergy) {
00232
00233 G4double chargeSquareRatio = corrections ->
00234 EffectiveChargeSquareRatio(particle,
00235 material,
00236 kineticEnergy);
00237 corrFactor = chargeSquareRatio *
00238 corrections -> EffectiveChargeCorrection(particle,
00239 material,
00240 kineticEnergy);
00241 return corrFactor;
00242 }
00243
00244
00245
00246 G4double G4IonParametrisedLossModel::GetParticleCharge(
00247 const G4ParticleDefinition* particle,
00248 const G4Material* material,
00249 G4double kineticEnergy) {
00250
00251 return corrections -> GetParticleCharge(particle, material, kineticEnergy);
00252 }
00253
00254
00255
00256 void G4IonParametrisedLossModel::Initialise(
00257 const G4ParticleDefinition* particle,
00258 const G4DataVector& cuts) {
00259
00260
00261 cacheParticle = 0;
00262 cacheMass = 0;
00263 cacheElecMassRatio = 0;
00264 cacheChargeSquare = 0;
00265
00266
00267 rangeCacheParticle = 0;
00268 rangeCacheMatCutsCouple = 0;
00269 rangeCacheEnergyRange = 0;
00270 rangeCacheRangeEnergy = 0;
00271
00272
00273 dedxCacheParticle = 0;
00274 dedxCacheMaterial = 0;
00275 dedxCacheEnergyCut = 0;
00276 dedxCacheIter = lossTableList.end();
00277 dedxCacheTransitionEnergy = 0.0;
00278 dedxCacheTransitionFactor = 0.0;
00279 dedxCacheGenIonMassRatio = 0.0;
00280
00281
00282 LossTableList::iterator iterTables = lossTableList.begin();
00283 LossTableList::iterator iterTables_end = lossTableList.end();
00284
00285 for(;iterTables != iterTables_end; iterTables++)
00286 (*iterTables) -> ClearCache();
00287
00288
00289
00290 RangeEnergyTable::iterator iterRange = r.begin();
00291 RangeEnergyTable::iterator iterRange_end = r.end();
00292
00293 for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second;
00294 r.clear();
00295
00296 EnergyRangeTable::iterator iterEnergy = E.begin();
00297 EnergyRangeTable::iterator iterEnergy_end = E.end();
00298
00299 for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second;
00300 E.clear();
00301
00302
00303 size_t size = cuts.size();
00304 cutEnergies.clear();
00305 for(size_t i = 0; i < size; i++) cutEnergies.push_back(cuts[i]);
00306
00307
00308 const G4ProductionCutsTable* coupleTable=
00309 G4ProductionCutsTable::GetProductionCutsTable();
00310 size_t nmbCouples = coupleTable -> GetTableSize();
00311
00312 #ifdef PRINT_TABLE_BUILT
00313 G4cout << "G4IonParametrisedLossModel::Initialise():"
00314 << " Building dE/dx vectors:"
00315 << G4endl;
00316 #endif
00317
00318 for (size_t i = 0; i < nmbCouples; i++) {
00319
00320 const G4MaterialCutsCouple* couple =
00321 coupleTable -> GetMaterialCutsCouple(i);
00322
00323 const G4Material* material = couple -> GetMaterial();
00324
00325
00326 for(G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) {
00327
00328 LossTableList::iterator iter = lossTableList.begin();
00329 LossTableList::iterator iter_end = lossTableList.end();
00330
00331 for(;iter != iter_end; iter++) {
00332
00333 if(*iter == 0) {
00334 G4cout << "G4IonParametrisedLossModel::Initialise():"
00335 << " Skipping illegal table."
00336 << G4endl;
00337 }
00338
00339 G4bool isApplicable =
00340 (*iter) -> BuildDEDXTable(atomicNumberIon, material);
00341 if(isApplicable) {
00342
00343 #ifdef PRINT_TABLE_BUILT
00344 G4cout << " Atomic Number Ion = " << atomicNumberIon
00345 << ", Material = " << material -> GetName()
00346 << ", Table = " << (*iter) -> GetName()
00347 << G4endl;
00348 #endif
00349 break;
00350 }
00351 }
00352 }
00353 }
00354
00355
00356 if(! particleChangeLoss) {
00357 particleChangeLoss = GetParticleChangeForLoss();
00358 braggIonModel -> SetParticleChange(particleChangeLoss, 0);
00359 betheBlochModel -> SetParticleChange(particleChangeLoss, 0);
00360 }
00361
00362
00363
00364 braggIonModel -> Initialise(particle, cuts);
00365 betheBlochModel -> Initialise(particle, cuts);
00366 }
00367
00368
00369
00370 G4double G4IonParametrisedLossModel::ComputeCrossSectionPerAtom(
00371 const G4ParticleDefinition* particle,
00372 G4double kineticEnergy,
00373 G4double atomicNumber,
00374 G4double,
00375 G4double cutEnergy,
00376 G4double maxKinEnergy) {
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 G4double crosssection = 0.0;
00389 G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
00390 G4double maxEnergy = std::min(tmax, maxKinEnergy);
00391
00392 if(cutEnergy < tmax) {
00393
00394 G4double energy = kineticEnergy + cacheMass;
00395 G4double betaSquared = kineticEnergy *
00396 (energy + cacheMass) / (energy * energy);
00397
00398 crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
00399 betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
00400
00401 crosssection *= twopi_mc2_rcl2 * cacheChargeSquare / betaSquared;
00402 }
00403
00404 #ifdef PRINT_DEBUG_CS
00405 G4cout << "########################################################"
00406 << G4endl
00407 << "# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
00408 << G4endl
00409 << "# particle =" << particle -> GetParticleName()
00410 << G4endl
00411 << "# cut(MeV) = " << cutEnergy/MeV
00412 << G4endl;
00413
00414 G4cout << "#"
00415 << std::setw(13) << std::right << "E(MeV)"
00416 << std::setw(14) << "CS(um)"
00417 << std::setw(14) << "E_max_sec(MeV)"
00418 << G4endl
00419 << "# ------------------------------------------------------"
00420 << G4endl;
00421
00422 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
00423 << std::setw(14) << crosssection / (um * um)
00424 << std::setw(14) << tmax / MeV
00425 << G4endl;
00426 #endif
00427
00428 crosssection *= atomicNumber;
00429
00430 return crosssection;
00431 }
00432
00433
00434
00435 G4double G4IonParametrisedLossModel::CrossSectionPerVolume(
00436 const G4Material* material,
00437 const G4ParticleDefinition* particle,
00438 G4double kineticEnergy,
00439 G4double cutEnergy,
00440 G4double maxEnergy) {
00441
00442 G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
00443 G4double cross = ComputeCrossSectionPerAtom(particle,
00444 kineticEnergy,
00445 nbElecPerVolume, 0,
00446 cutEnergy,
00447 maxEnergy);
00448
00449 return cross;
00450 }
00451
00452
00453
00454 G4double G4IonParametrisedLossModel::ComputeDEDXPerVolume(
00455 const G4Material* material,
00456 const G4ParticleDefinition* particle,
00457 G4double kineticEnergy,
00458 G4double cutEnergy) {
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487 G4double dEdx = 0.0;
00488 UpdateDEDXCache(particle, material, cutEnergy);
00489
00490 LossTableList::iterator iter = dedxCacheIter;
00491
00492 if(iter != lossTableList.end()) {
00493
00494 G4double transitionEnergy = dedxCacheTransitionEnergy;
00495
00496 if(transitionEnergy > kineticEnergy) {
00497
00498 dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
00499
00500 G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material,
00501 particle,
00502 kineticEnergy,
00503 cutEnergy);
00504 dEdx -= dEdxDeltaRays;
00505 }
00506 else {
00507 G4double massRatio = dedxCacheGenIonMassRatio;
00508
00509 G4double chargeSquare =
00510 GetChargeSquareRatio(particle, material, kineticEnergy);
00511
00512 G4double scaledKineticEnergy = kineticEnergy * massRatio;
00513 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
00514
00515 G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
00516
00517 if(scaledTransitionEnergy >= lowEnergyLimit) {
00518
00519 dEdx = betheBlochModel -> ComputeDEDXPerVolume(
00520 material, genericIon,
00521 scaledKineticEnergy, cutEnergy);
00522
00523 dEdx *= chargeSquare;
00524
00525 dEdx += corrections -> ComputeIonCorrections(particle,
00526 material, kineticEnergy);
00527
00528 G4double factor = 1.0 + dedxCacheTransitionFactor /
00529 kineticEnergy;
00530
00531 dEdx *= factor;
00532 }
00533 }
00534 }
00535 else {
00536 G4double massRatio = 1.0;
00537 G4double chargeSquare = 1.0;
00538
00539 if(particle != genericIon) {
00540
00541 chargeSquare = GetChargeSquareRatio(particle, material, kineticEnergy);
00542 massRatio = genericIonPDGMass / particle -> GetPDGMass();
00543 }
00544
00545 G4double scaledKineticEnergy = kineticEnergy * massRatio;
00546
00547 G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
00548 if(scaledKineticEnergy < lowEnergyLimit) {
00549 dEdx = braggIonModel -> ComputeDEDXPerVolume(
00550 material, genericIon,
00551 scaledKineticEnergy, cutEnergy);
00552
00553 dEdx *= chargeSquare;
00554 }
00555 else {
00556 G4double dEdxLimitParam = braggIonModel -> ComputeDEDXPerVolume(
00557 material, genericIon,
00558 lowEnergyLimit, cutEnergy);
00559
00560 G4double dEdxLimitBetheBloch = betheBlochModel -> ComputeDEDXPerVolume(
00561 material, genericIon,
00562 lowEnergyLimit, cutEnergy);
00563
00564 if(particle != genericIon) {
00565 G4double chargeSquareLowEnergyLimit =
00566 GetChargeSquareRatio(particle, material,
00567 lowEnergyLimit / massRatio);
00568
00569 dEdxLimitParam *= chargeSquareLowEnergyLimit;
00570 dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
00571
00572 dEdxLimitBetheBloch +=
00573 corrections -> ComputeIonCorrections(particle,
00574 material, lowEnergyLimit / massRatio);
00575 }
00576
00577 G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
00578 * lowEnergyLimit / scaledKineticEnergy);
00579
00580 dEdx = betheBlochModel -> ComputeDEDXPerVolume(
00581 material, genericIon,
00582 scaledKineticEnergy, cutEnergy);
00583
00584 dEdx *= chargeSquare;
00585
00586 if(particle != genericIon) {
00587 dEdx += corrections -> ComputeIonCorrections(particle,
00588 material, kineticEnergy);
00589 }
00590
00591 dEdx *= factor;
00592 }
00593
00594 }
00595
00596 if (dEdx < 0.0) dEdx = 0.0;
00597
00598 return dEdx;
00599 }
00600
00601
00602
00603 void G4IonParametrisedLossModel::PrintDEDXTable(
00604 const G4ParticleDefinition* particle,
00605 const G4Material* material,
00606 G4double lowerBoundary,
00607 G4double upperBoundary,
00608 G4int numBins,
00609 G4bool logScaleEnergy) {
00610
00611 G4double atomicMassNumber = particle -> GetAtomicMass();
00612 G4double materialDensity = material -> GetDensity();
00613
00614 G4cout << "# dE/dx table for " << particle -> GetParticleName()
00615 << " in material " << material -> GetName()
00616 << " of density " << materialDensity / g * cm3
00617 << " g/cm3"
00618 << G4endl
00619 << "# Projectile mass number A1 = " << atomicMassNumber
00620 << G4endl
00621 << "# ------------------------------------------------------"
00622 << G4endl;
00623 G4cout << "#"
00624 << std::setw(13) << std::right << "E"
00625 << std::setw(14) << "E/A1"
00626 << std::setw(14) << "dE/dx"
00627 << std::setw(14) << "1/rho*dE/dx"
00628 << G4endl;
00629 G4cout << "#"
00630 << std::setw(13) << std::right << "(MeV)"
00631 << std::setw(14) << "(MeV)"
00632 << std::setw(14) << "(MeV/cm)"
00633 << std::setw(14) << "(MeV*cm2/mg)"
00634 << G4endl
00635 << "# ------------------------------------------------------"
00636 << G4endl;
00637
00638 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
00639 G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
00640
00641 if(logScaleEnergy) {
00642
00643 energyLowerBoundary = std::log(energyLowerBoundary);
00644 energyUpperBoundary = std::log(energyUpperBoundary);
00645 }
00646
00647 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
00648 G4double(nmbBins);
00649
00650 for(int i = 0; i < numBins + 1; i++) {
00651
00652 G4double energy = energyLowerBoundary + i * deltaEnergy;
00653 if(logScaleEnergy) energy = std::exp(energy);
00654
00655 G4double dedx = ComputeDEDXPerVolume(material, particle, energy, DBL_MAX);
00656 G4cout.precision(6);
00657 G4cout << std::setw(14) << std::right << energy / MeV
00658 << std::setw(14) << energy / atomicMassNumber / MeV
00659 << std::setw(14) << dedx / MeV * cm
00660 << std::setw(14) << dedx / materialDensity / (MeV*cm2/(0.001*g))
00661 << G4endl;
00662 }
00663 }
00664
00665
00666
00667 void G4IonParametrisedLossModel::PrintDEDXTableHandlers(
00668 const G4ParticleDefinition* particle,
00669 const G4Material* material,
00670 G4double lowerBoundary,
00671 G4double upperBoundary,
00672 G4int numBins,
00673 G4bool logScaleEnergy) {
00674
00675 LossTableList::iterator iter = lossTableList.begin();
00676 LossTableList::iterator iter_end = lossTableList.end();
00677
00678 for(;iter != iter_end; iter++) {
00679 G4bool isApplicable = (*iter) -> IsApplicable(particle, material);
00680 if(isApplicable) {
00681 (*iter) -> PrintDEDXTable(particle, material,
00682 lowerBoundary, upperBoundary,
00683 numBins,logScaleEnergy);
00684 break;
00685 }
00686 }
00687 }
00688
00689
00690
00691 void G4IonParametrisedLossModel::SampleSecondaries(
00692 std::vector<G4DynamicParticle*>* secondaries,
00693 const G4MaterialCutsCouple*,
00694 const G4DynamicParticle* particle,
00695 G4double cutKinEnergySec,
00696 G4double userMaxKinEnergySec) {
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724 G4double rossiMaxKinEnergySec = MaxSecondaryKinEnergy(particle);
00725 G4double maxKinEnergySec =
00726 std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
00727
00728 if(cutKinEnergySec >= maxKinEnergySec) return;
00729
00730 G4double kineticEnergy = particle -> GetKineticEnergy();
00731 G4ThreeVector direction = particle ->GetMomentumDirection();
00732
00733 G4double energy = kineticEnergy + cacheMass;
00734 G4double betaSquared = kineticEnergy *
00735 (energy + cacheMass) / (energy * energy);
00736
00737 G4double kinEnergySec;
00738 G4double grej;
00739
00740 do {
00741
00742
00743 G4double xi = G4UniformRand();
00744 kinEnergySec = cutKinEnergySec * maxKinEnergySec /
00745 (maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
00746
00747
00748
00749 grej = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
00750
00751 if(grej > 1.0) {
00752 G4cout << "G4IonParametrisedLossModel::SampleSecondary Warning: "
00753 << "Majorant 1.0 < "
00754 << grej << " for e= " << kinEnergySec
00755 << G4endl;
00756 }
00757
00758 } while( G4UniformRand() >= grej );
00759
00760 G4double momentumSec =
00761 std::sqrt(kinEnergySec * (kinEnergySec + 2.0 * electron_mass_c2));
00762
00763 G4double totMomentum = energy*std::sqrt(betaSquared);
00764 G4double cost = kinEnergySec * (energy + electron_mass_c2) /
00765 (momentumSec * totMomentum);
00766 if(cost > 1.0) cost = 1.0;
00767 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
00768
00769 G4double phi = twopi * G4UniformRand() ;
00770
00771 G4ThreeVector directionSec(sint*std::cos(phi),sint*std::sin(phi), cost) ;
00772 directionSec.rotateUz(direction);
00773
00774
00775 G4DynamicParticle* delta = new G4DynamicParticle(G4Electron::Definition(),
00776 directionSec,
00777 kinEnergySec);
00778
00779 secondaries -> push_back(delta);
00780
00781
00782 kineticEnergy -= kinEnergySec;
00783 G4ThreeVector finalP = direction*totMomentum - directionSec*momentumSec;
00784 finalP = finalP.unit();
00785
00786 particleChangeLoss -> SetProposedKineticEnergy(kineticEnergy);
00787 particleChangeLoss -> SetProposedMomentumDirection(finalP);
00788 }
00789
00790
00791
00792 void G4IonParametrisedLossModel::UpdateRangeCache(
00793 const G4ParticleDefinition* particle,
00794 const G4MaterialCutsCouple* matCutsCouple) {
00795
00796
00797
00798
00799
00800 if(particle == rangeCacheParticle &&
00801 matCutsCouple == rangeCacheMatCutsCouple) {
00802 }
00803 else{
00804 rangeCacheParticle = particle;
00805 rangeCacheMatCutsCouple = matCutsCouple;
00806
00807 const G4Material* material = matCutsCouple -> GetMaterial();
00808 LossTableList::iterator iter = IsApplicable(particle, material);
00809
00810
00811 if(iter != lossTableList.end()) {
00812
00813
00814 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
00815 RangeEnergyTable::iterator iterRange = r.find(ionMatCouple);
00816
00817 if(iterRange == r.end()) BuildRangeVector(particle, matCutsCouple);
00818
00819 rangeCacheEnergyRange = E[ionMatCouple];
00820 rangeCacheRangeEnergy = r[ionMatCouple];
00821 }
00822 else {
00823 rangeCacheEnergyRange = 0;
00824 rangeCacheRangeEnergy = 0;
00825 }
00826 }
00827 }
00828
00829
00830
00831 void G4IonParametrisedLossModel::UpdateDEDXCache(
00832 const G4ParticleDefinition* particle,
00833 const G4Material* material,
00834 G4double cutEnergy) {
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844 if(particle == dedxCacheParticle &&
00845 material == dedxCacheMaterial &&
00846 cutEnergy == dedxCacheEnergyCut) {
00847 }
00848 else {
00849
00850 dedxCacheParticle = particle;
00851 dedxCacheMaterial = material;
00852 dedxCacheEnergyCut = cutEnergy;
00853
00854 G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
00855 dedxCacheGenIonMassRatio = massRatio;
00856
00857 LossTableList::iterator iter = IsApplicable(particle, material);
00858 dedxCacheIter = iter;
00859
00860
00861 if(iter != lossTableList.end()) {
00862
00863
00864 G4double transitionEnergy =
00865 (*iter) -> GetUpperEnergyEdge(particle, material);
00866 dedxCacheTransitionEnergy = transitionEnergy;
00867
00868
00869
00870 G4double dEdxParam = (*iter) -> GetDEDX(particle, material,
00871 transitionEnergy);
00872
00873 G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material,
00874 particle,
00875 transitionEnergy,
00876 cutEnergy);
00877 dEdxParam -= dEdxDeltaRays;
00878
00879
00880
00881 G4double transitionChargeSquare =
00882 GetChargeSquareRatio(particle, material, transitionEnergy);
00883
00884 G4double scaledTransitionEnergy = transitionEnergy * massRatio;
00885
00886 G4double dEdxBetheBloch =
00887 betheBlochModel -> ComputeDEDXPerVolume(
00888 material, genericIon,
00889 scaledTransitionEnergy, cutEnergy);
00890 dEdxBetheBloch *= transitionChargeSquare;
00891
00892
00893 dEdxBetheBloch +=
00894 corrections -> ComputeIonCorrections(particle,
00895 material, transitionEnergy);
00896
00897
00898 dedxCacheTransitionFactor =
00899 (dEdxParam - dEdxBetheBloch)/dEdxBetheBloch
00900 * transitionEnergy;
00901 }
00902 else {
00903
00904 dedxCacheParticle = particle;
00905 dedxCacheMaterial = material;
00906 dedxCacheEnergyCut = cutEnergy;
00907
00908 dedxCacheGenIonMassRatio =
00909 genericIonPDGMass / particle -> GetPDGMass();
00910
00911 dedxCacheTransitionEnergy = 0.0;
00912 dedxCacheTransitionFactor = 0.0;
00913 }
00914 }
00915 }
00916
00917
00918
00919 void G4IonParametrisedLossModel::CorrectionsAlongStep(
00920 const G4MaterialCutsCouple* couple,
00921 const G4DynamicParticle* dynamicParticle,
00922 G4double& eloss,
00923 G4double&,
00924 G4double length) {
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938 const G4ParticleDefinition* particle = dynamicParticle -> GetDefinition();
00939 const G4Material* material = couple -> GetMaterial();
00940
00941 G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
00942
00943 if(kineticEnergy == eloss) { return; }
00944
00945 G4double cutEnergy = DBL_MAX;
00946 size_t cutIndex = couple -> GetIndex();
00947 cutEnergy = cutEnergies[cutIndex];
00948
00949 UpdateDEDXCache(particle, material, cutEnergy);
00950
00951 LossTableList::iterator iter = dedxCacheIter;
00952
00953
00954
00955 if(iter != lossTableList.end()) {
00956
00957
00958
00959
00960 eloss =
00961 length * ComputeDEDXPerVolume(material, particle,
00962 kineticEnergy, cutEnergy);
00963
00964 #ifdef PRINT_DEBUG
00965 G4cout.precision(6);
00966 G4cout << "########################################################"
00967 << G4endl
00968 << "# G4IonParametrisedLossModel::CorrectionsAlongStep"
00969 << G4endl
00970 << "# cut(MeV) = " << cutEnergy/MeV
00971 << G4endl;
00972
00973 G4cout << "#"
00974 << std::setw(13) << std::right << "E(MeV)"
00975 << std::setw(14) << "l(um)"
00976 << std::setw(14) << "l*dE/dx(MeV)"
00977 << std::setw(14) << "(l*dE/dx)/E"
00978 << G4endl
00979 << "# ------------------------------------------------------"
00980 << G4endl;
00981
00982 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
00983 << std::setw(14) << length / um
00984 << std::setw(14) << eloss / MeV
00985 << std::setw(14) << eloss / kineticEnergy * 100.0
00986 << G4endl;
00987 #endif
00988
00989
00990
00991
00992
00993 if(eloss > energyLossLimit * kineticEnergy) {
00994
00995 eloss = ComputeLossForStep(couple, particle,
00996 kineticEnergy,length);
00997
00998 #ifdef PRINT_DEBUG
00999 G4cout << "# Correction applied:"
01000 << G4endl;
01001
01002 G4cout << std::setw(14) << std::right << kineticEnergy / MeV
01003 << std::setw(14) << length / um
01004 << std::setw(14) << eloss / MeV
01005 << std::setw(14) << eloss / kineticEnergy * 100.0
01006 << G4endl;
01007 #endif
01008
01009 }
01010 }
01011
01012
01013
01014 G4double energy = kineticEnergy - eloss * 0.5;
01015 if(energy < 0.0) energy = kineticEnergy * 0.5;
01016
01017 G4double chargeSquareRatio = corrections ->
01018 EffectiveChargeSquareRatio(particle,
01019 material,
01020 energy);
01021 GetModelOfFluctuations() -> SetParticleAndCharge(particle,
01022 chargeSquareRatio);
01023
01024
01025
01026
01027
01028
01029 G4double transitionEnergy = dedxCacheTransitionEnergy;
01030
01031 if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
01032 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
01033 material,
01034 energy);
01035
01036 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
01037 eloss *= chargeSquareRatioCorr;
01038 }
01039 else if (iter == lossTableList.end()) {
01040
01041 chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
01042 material,
01043 energy);
01044
01045 G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
01046 eloss *= chargeSquareRatioCorr;
01047 }
01048
01049
01050
01051
01052 if(iter == lossTableList.end()) {
01053
01054 G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
01055 G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
01056
01057
01058 if(scaledKineticEnergy > lowEnergyLimit)
01059 eloss += length *
01060 corrections -> IonHighOrderCorrections(particle, couple, energy);
01061 }
01062 }
01063
01064
01065
01066 void G4IonParametrisedLossModel::BuildRangeVector(
01067 const G4ParticleDefinition* particle,
01068 const G4MaterialCutsCouple* matCutsCouple) {
01069
01070 G4double cutEnergy = DBL_MAX;
01071 size_t cutIndex = matCutsCouple -> GetIndex();
01072 cutEnergy = cutEnergies[cutIndex];
01073
01074 const G4Material* material = matCutsCouple -> GetMaterial();
01075
01076 G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
01077
01078 G4double lowerEnergy = lowerEnergyEdgeIntegr / massRatio;
01079 G4double upperEnergy = upperEnergyEdgeIntegr / massRatio;
01080
01081 G4double logLowerEnergyEdge = std::log(lowerEnergy);
01082 G4double logUpperEnergyEdge = std::log(upperEnergy);
01083
01084 G4double logDeltaEnergy = (logUpperEnergyEdge - logLowerEnergyEdge) /
01085 G4double(nmbBins);
01086
01087 G4double logDeltaIntegr = logDeltaEnergy / G4double(nmbSubBins);
01088
01089 G4LPhysicsFreeVector* energyRangeVector =
01090 new G4LPhysicsFreeVector(nmbBins+1,
01091 lowerEnergy,
01092 upperEnergy);
01093 energyRangeVector -> SetSpline(true);
01094
01095 G4double dedxLow = ComputeDEDXPerVolume(material,
01096 particle,
01097 lowerEnergy,
01098 cutEnergy);
01099
01100 G4double range = 2.0 * lowerEnergy / dedxLow;
01101
01102 energyRangeVector -> PutValues(0, lowerEnergy, range);
01103
01104 G4double logEnergy = std::log(lowerEnergy);
01105 for(size_t i = 1; i < nmbBins+1; i++) {
01106
01107 G4double logEnergyIntegr = logEnergy;
01108
01109 for(size_t j = 0; j < nmbSubBins; j++) {
01110
01111 G4double binLowerBoundary = std::exp(logEnergyIntegr);
01112 logEnergyIntegr += logDeltaIntegr;
01113
01114 G4double binUpperBoundary = std::exp(logEnergyIntegr);
01115 G4double deltaIntegr = binUpperBoundary - binLowerBoundary;
01116
01117 G4double energyIntegr = binLowerBoundary + 0.5 * deltaIntegr;
01118
01119 G4double dedxValue = ComputeDEDXPerVolume(material,
01120 particle,
01121 energyIntegr,
01122 cutEnergy);
01123
01124 if(dedxValue > 0.0) range += deltaIntegr / dedxValue;
01125
01126 #ifdef PRINT_DEBUG_DETAILS
01127 G4cout << " E = "<< energyIntegr/MeV
01128 << " MeV -> dE = " << deltaIntegr/MeV
01129 << " MeV -> dE/dx = " << dedxValue/MeV*mm
01130 << " MeV/mm -> dE/(dE/dx) = " << deltaIntegr /
01131 dedxValue / mm
01132 << " mm -> range = " << range / mm
01133 << " mm " << G4endl;
01134 #endif
01135 }
01136
01137 logEnergy += logDeltaEnergy;
01138
01139 G4double energy = std::exp(logEnergy);
01140
01141 energyRangeVector -> PutValues(i, energy, range);
01142
01143 #ifdef PRINT_DEBUG_DETAILS
01144 G4cout << "G4IonParametrisedLossModel::BuildRangeVector() bin = "
01145 << i <<", E = "
01146 << energy / MeV << " MeV, R = "
01147 << range / mm << " mm"
01148 << G4endl;
01149 #endif
01150
01151 }
01152
01153 G4bool b;
01154
01155 G4double lowerRangeEdge =
01156 energyRangeVector -> GetValue(lowerEnergy, b);
01157 G4double upperRangeEdge =
01158 energyRangeVector -> GetValue(upperEnergy, b);
01159
01160 G4LPhysicsFreeVector* rangeEnergyVector
01161 = new G4LPhysicsFreeVector(nmbBins+1,
01162 lowerRangeEdge,
01163 upperRangeEdge);
01164 rangeEnergyVector -> SetSpline(true);
01165
01166 for(size_t i = 0; i < nmbBins+1; i++) {
01167 G4double energy = energyRangeVector -> GetLowEdgeEnergy(i);
01168 rangeEnergyVector ->
01169 PutValues(i, energyRangeVector -> GetValue(energy, b), energy);
01170 }
01171
01172 #ifdef PRINT_DEBUG_TABLES
01173 G4cout << *energyLossVector
01174 << *energyRangeVector
01175 << *rangeEnergyVector << G4endl;
01176 #endif
01177
01178 IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
01179
01180 E[ionMatCouple] = energyRangeVector;
01181 r[ionMatCouple] = rangeEnergyVector;
01182 }
01183
01184
01185
01186 G4double G4IonParametrisedLossModel::ComputeLossForStep(
01187 const G4MaterialCutsCouple* matCutsCouple,
01188 const G4ParticleDefinition* particle,
01189 G4double kineticEnergy,
01190 G4double stepLength) {
01191
01192 G4double loss = 0.0;
01193
01194 UpdateRangeCache(particle, matCutsCouple);
01195
01196 G4PhysicsVector* energyRange = rangeCacheEnergyRange;
01197 G4PhysicsVector* rangeEnergy = rangeCacheRangeEnergy;
01198
01199 if(energyRange != 0 && rangeEnergy != 0) {
01200 G4bool b;
01201
01202 G4double lowerEnEdge = energyRange -> GetLowEdgeEnergy( 0 );
01203 G4double lowerRangeEdge = rangeEnergy -> GetLowEdgeEnergy( 0 );
01204
01205
01206 G4double range = energyRange -> GetValue(kineticEnergy, b);
01207
01208
01209 if(kineticEnergy < lowerEnEdge) {
01210
01211 range = energyRange -> GetValue(lowerEnEdge, b);
01212 range *= std::sqrt(kineticEnergy / lowerEnEdge);
01213 }
01214
01215 #ifdef PRINT_DEBUG
01216 G4cout << "G4IonParametrisedLossModel::ComputeLossForStep() range = "
01217 << range / mm << " mm, step = " << stepLength / mm << " mm"
01218 << G4endl;
01219 #endif
01220
01221
01222 G4double remRange = range - stepLength;
01223
01224
01225
01226 if(remRange < 0.0) loss = kineticEnergy;
01227 else if(remRange < lowerRangeEdge) {
01228
01229 G4double ratio = remRange / lowerRangeEdge;
01230 loss = kineticEnergy - ratio * ratio * lowerEnEdge;
01231 }
01232 else {
01233
01234 G4double energy = rangeEnergy -> GetValue(range - stepLength, b);
01235 loss = kineticEnergy - energy;
01236 }
01237 }
01238
01239 if(loss < 0.0) loss = 0.0;
01240
01241 return loss;
01242 }
01243
01244
01245
01246 G4bool G4IonParametrisedLossModel::AddDEDXTable(
01247 const G4String& nam,
01248 G4VIonDEDXTable* table,
01249 G4VIonDEDXScalingAlgorithm* algorithm) {
01250
01251 if(table == 0) {
01252 G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
01253 << " add table: Invalid pointer."
01254 << G4endl;
01255
01256 return false;
01257 }
01258
01259
01260 LossTableList::iterator iter = lossTableList.begin();
01261 LossTableList::iterator iter_end = lossTableList.end();
01262
01263 for(;iter != iter_end; iter++) {
01264 G4String tableName = (*iter) -> GetName();
01265
01266 if(tableName == nam) {
01267 G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
01268 << " add table: Name already exists."
01269 << G4endl;
01270
01271 return false;
01272 }
01273 }
01274
01275 G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm;
01276 if(scalingAlgorithm == 0)
01277 scalingAlgorithm = new G4VIonDEDXScalingAlgorithm;
01278
01279 G4IonDEDXHandler* handler =
01280 new G4IonDEDXHandler(table, scalingAlgorithm, nam);
01281
01282 lossTableList.push_front(handler);
01283
01284 return true;
01285 }
01286
01287
01288
01289 G4bool G4IonParametrisedLossModel::RemoveDEDXTable(
01290 const G4String& nam) {
01291
01292 LossTableList::iterator iter = lossTableList.begin();
01293 LossTableList::iterator iter_end = lossTableList.end();
01294
01295 for(;iter != iter_end; iter++) {
01296 G4String tableName = (*iter) -> GetName();
01297
01298 if(tableName == nam) {
01299 delete (*iter);
01300
01301
01302 lossTableList.erase(iter);
01303
01304
01305 RangeEnergyTable::iterator iterRange = r.begin();
01306 RangeEnergyTable::iterator iterRange_end = r.end();
01307
01308 for(;iterRange != iterRange_end; iterRange++)
01309 delete iterRange -> second;
01310 r.clear();
01311
01312 EnergyRangeTable::iterator iterEnergy = E.begin();
01313 EnergyRangeTable::iterator iterEnergy_end = E.end();
01314
01315 for(;iterEnergy != iterEnergy_end; iterEnergy++)
01316 delete iterEnergy -> second;
01317 E.clear();
01318
01319 return true;
01320 }
01321 }
01322
01323 return false;
01324 }
01325
01326
01327
01328 void G4IonParametrisedLossModel::DeactivateICRU73Scaling() {
01329
01330 RemoveDEDXTable("ICRU73");
01331 AddDEDXTable("ICRU73", new G4IonStoppingData("ion_stopping_data/icru73"));
01332 }
01333
01334