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 #include "G4BraggIonModel.hh"
00062 #include "G4PhysicalConstants.hh"
00063 #include "G4SystemOfUnits.hh"
00064 #include "Randomize.hh"
00065 #include "G4Electron.hh"
00066 #include "G4ParticleChangeForLoss.hh"
00067 #include "G4LossTableManager.hh"
00068 #include "G4EmCorrections.hh"
00069
00070
00071
00072 using namespace std;
00073
00074 G4BraggIonModel::G4BraggIonModel(const G4ParticleDefinition* p,
00075 const G4String& nam)
00076 : G4VEmModel(nam),
00077 corr(0),
00078 particle(0),
00079 fParticleChange(0),
00080 currentMaterial(0),
00081 iMolecula(-1),
00082 iASTAR(-1),
00083 isIon(false),
00084 isInitialised(false)
00085 {
00086 SetHighEnergyLimit(2.0*MeV);
00087
00088 HeMass = 3.727417*GeV;
00089 rateMassHe2p = HeMass/proton_mass_c2;
00090 lowestKinEnergy = 1.0*keV/rateMassHe2p;
00091 massFactor = 1000.*amu_c2/HeMass;
00092 theZieglerFactor = eV*cm2*1.0e-15;
00093 theElectron = G4Electron::Electron();
00094 corrFactor = 1.0;
00095 if(p) { SetParticle(p); }
00096 else { SetParticle(theElectron); }
00097 }
00098
00099
00100
00101 G4BraggIonModel::~G4BraggIonModel()
00102 {}
00103
00104
00105
00106 void G4BraggIonModel::Initialise(const G4ParticleDefinition* p,
00107 const G4DataVector&)
00108 {
00109 if(p != particle) { SetParticle(p); }
00110
00111 corrFactor = chargeSquare;
00112
00113
00114 SetDeexcitationFlag(false);
00115
00116 if(!isInitialised) {
00117 isInitialised = true;
00118
00119 G4String pname = particle->GetParticleName();
00120 if(particle->GetParticleType() == "nucleus" &&
00121 pname != "deuteron" && pname != "triton" &&
00122 pname != "alpha+" && pname != "helium" &&
00123 pname != "hydrogen") { isIon = true; }
00124
00125 corr = G4LossTableManager::Instance()->EmCorrections();
00126
00127 fParticleChange = GetParticleChangeForLoss();
00128 }
00129 }
00130
00131
00132
00133 G4double G4BraggIonModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
00134 const G4Material* mat,
00135 G4double kineticEnergy)
00136 {
00137
00138
00139 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
00140 corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
00141 return corrFactor;
00142 }
00143
00144
00145
00146 G4double G4BraggIonModel::GetParticleCharge(const G4ParticleDefinition* p,
00147 const G4Material* mat,
00148 G4double kineticEnergy)
00149 {
00150
00151
00152 return corr->GetParticleCharge(p,mat,kineticEnergy);
00153 }
00154
00155
00156
00157 G4double G4BraggIonModel::ComputeCrossSectionPerElectron(
00158 const G4ParticleDefinition* p,
00159 G4double kineticEnergy,
00160 G4double cutEnergy,
00161 G4double maxKinEnergy)
00162 {
00163 G4double cross = 0.0;
00164 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00165 G4double maxEnergy = std::min(tmax,maxKinEnergy);
00166 if(cutEnergy < tmax) {
00167
00168 G4double energy = kineticEnergy + mass;
00169 G4double energy2 = energy*energy;
00170 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
00171 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax;
00172
00173 cross *= twopi_mc2_rcl2*chargeSquare/beta2;
00174 }
00175
00176
00177
00178 return cross;
00179 }
00180
00181
00182
00183 G4double G4BraggIonModel::ComputeCrossSectionPerAtom(
00184 const G4ParticleDefinition* p,
00185 G4double kineticEnergy,
00186 G4double Z, G4double,
00187 G4double cutEnergy,
00188 G4double maxEnergy)
00189 {
00190 G4double cross = Z*ComputeCrossSectionPerElectron
00191 (p,kineticEnergy,cutEnergy,maxEnergy);
00192 return cross;
00193 }
00194
00195
00196
00197 G4double G4BraggIonModel::CrossSectionPerVolume(
00198 const G4Material* material,
00199 const G4ParticleDefinition* p,
00200 G4double kineticEnergy,
00201 G4double cutEnergy,
00202 G4double maxEnergy)
00203 {
00204 G4double eDensity = material->GetElectronDensity();
00205 G4double cross = eDensity*ComputeCrossSectionPerElectron
00206 (p,kineticEnergy,cutEnergy,maxEnergy);
00207 return cross;
00208 }
00209
00210
00211
00212 G4double G4BraggIonModel::ComputeDEDXPerVolume(const G4Material* material,
00213 const G4ParticleDefinition* p,
00214 G4double kineticEnergy,
00215 G4double cutEnergy)
00216 {
00217 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00218 G4double tmin = min(cutEnergy, tmax);
00219 G4double tkin = kineticEnergy/massRate;
00220 G4double dedx = 0.0;
00221
00222 if(tkin < lowestKinEnergy) {
00223 dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
00224 } else {
00225 dedx = DEDX(material, tkin);
00226 }
00227
00228 if (cutEnergy < tmax) {
00229
00230 G4double tau = kineticEnergy/mass;
00231 G4double gam = tau + 1.0;
00232 G4double bg2 = tau * (tau+2.0);
00233 G4double beta2 = bg2/(gam*gam);
00234 G4double x = tmin/tmax;
00235
00236 dedx += (log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
00237 * (material->GetElectronDensity())/beta2;
00238 }
00239
00240
00241
00242 if (dedx < 0.0) dedx = 0.0 ;
00243
00244 dedx *= chargeSquare;
00245
00246
00247
00248
00249
00250 return dedx;
00251 }
00252
00253
00254
00255 void G4BraggIonModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
00256 const G4DynamicParticle* dp,
00257 G4double& eloss,
00258 G4double&,
00259 G4double )
00260 {
00261
00262 const G4ParticleDefinition* p = dp->GetDefinition();
00263 const G4Material* mat = couple->GetMaterial();
00264 G4double preKinEnergy = dp->GetKineticEnergy();
00265 G4double e = preKinEnergy - eloss*0.5;
00266 if(e < 0.0) { e = preKinEnergy*0.5; }
00267
00268 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
00269 GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
00270 G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
00271 eloss *= qfactor;
00272
00273
00274
00275 }
00276
00277
00278
00279 void G4BraggIonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
00280 const G4MaterialCutsCouple*,
00281 const G4DynamicParticle* dp,
00282 G4double xmin,
00283 G4double maxEnergy)
00284 {
00285 G4double tmax = MaxSecondaryKinEnergy(dp);
00286 G4double xmax = std::min(tmax, maxEnergy);
00287 if(xmin >= xmax) { return; }
00288
00289 G4double kineticEnergy = dp->GetKineticEnergy();
00290 G4double energy = kineticEnergy + mass;
00291 G4double energy2 = energy*energy;
00292 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
00293 G4double grej = 1.0;
00294 G4double deltaKinEnergy, f;
00295
00296 G4ThreeVector direction = dp->GetMomentumDirection();
00297
00298
00299 do {
00300 G4double q = G4UniformRand();
00301 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
00302
00303 f = 1.0 - beta2*deltaKinEnergy/tmax;
00304
00305 if(f > grej) {
00306 G4cout << "G4BraggIonModel::SampleSecondary Warning! "
00307 << "Majorant " << grej << " < "
00308 << f << " for e= " << deltaKinEnergy
00309 << G4endl;
00310 }
00311
00312 } while( grej*G4UniformRand() >= f );
00313
00314 G4double deltaMomentum =
00315 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
00316 G4double totMomentum = energy*sqrt(beta2);
00317 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
00318 (deltaMomentum * totMomentum);
00319 if(cost > 1.0) { cost = 1.0; }
00320 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
00321
00322 G4double phi = twopi * G4UniformRand() ;
00323
00324 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
00325 deltaDirection.rotateUz(direction);
00326
00327
00328 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,
00329 deltaKinEnergy);
00330
00331 vdp->push_back(delta);
00332
00333
00334 kineticEnergy -= deltaKinEnergy;
00335 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
00336 finalP = finalP.unit();
00337
00338 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00339 fParticleChange->SetProposedMomentumDirection(finalP);
00340 }
00341
00342
00343
00344 G4double G4BraggIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
00345 G4double kinEnergy)
00346 {
00347 if(pd != particle) { SetParticle(pd); }
00348 G4double tau = kinEnergy/mass;
00349 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
00350 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
00351 return tmax;
00352 }
00353
00354
00355
00356 G4bool G4BraggIonModel::HasMaterial(const G4Material* material)
00357 {
00358 G4String chFormula = material->GetChemicalFormula();
00359 if("" == chFormula) { return false; }
00360
00361
00362 const size_t numberOfMolecula = 11;
00363 const G4String molName[numberOfMolecula] = {
00364 "CaF_2", "Cellulose_Nitrate", "LiF", "Policarbonate",
00365 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polymethly_Methacralate",
00366 "Polysterene", "SiO_2", "NaI", "H_2O",
00367 "Graphite" };
00368 const G4int idxASTAR[numberOfMolecula] = {
00369 17, 19, 33, 51,
00370 52, 54,
00371 56, 62, 43, 71,
00372 13};
00373
00374
00375 for (size_t i=0; i<numberOfMolecula; ++i) {
00376 if (chFormula == molName[i]) {
00377 iMolecula = -1;
00378 iASTAR = idxASTAR[i];
00379 return true;
00380 }
00381 }
00382 return false ;
00383 }
00384
00385
00386
00387 G4double G4BraggIonModel::StoppingPower(const G4Material* material,
00388 G4double kineticEnergy)
00389 {
00390 G4double ionloss = 0.0 ;
00391
00392 if (iMolecula >= 0) {
00393
00394
00395
00396
00397
00398 G4double T = kineticEnergy*rateMassHe2p/MeV ;
00399
00400 static G4double a[11][5] = {
00401 {9.43672, 0.54398, 84.341, 1.3705, 57.422},
00402 {67.1503, 0.41409, 404.512, 148.97, 20.99},
00403 {5.11203, 0.453, 36.718, 50.6, 28.058},
00404 {61.793, 0.48445, 361.537, 57.889, 50.674},
00405 {7.83464, 0.49804, 160.452, 3.192, 0.71922},
00406 {19.729, 0.52153, 162.341, 58.35, 25.668},
00407 {26.4648, 0.50112, 188.913, 30.079, 16.509},
00408 {7.8655, 0.5205, 63.96, 51.32, 67.775},
00409 {8.8965, 0.5148, 339.36, 1.7205, 0.70423},
00410 {2.959, 0.53255, 34.247, 60.655, 15.153},
00411 {3.80133, 0.41590, 12.9966, 117.83, 242.28} };
00412
00413 static G4double atomicWeight[11] = {
00414 101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
00415 104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};
00416
00417 G4int i = iMolecula;
00418
00419
00420 if ( T < 0.001 ) {
00421 G4double slow = a[i][0] ;
00422 G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
00423 * a[i][2]*1000.0 ;
00424 ionloss = slow*shigh / (slow + shigh) ;
00425 ionloss *= sqrt(T*1000.0) ;
00426
00427
00428 } else {
00429 G4double slow = a[i][0] * pow((T*1000.0), a[i][1]) ;
00430 G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
00431 ionloss = slow*shigh / (slow + shigh) ;
00432
00433
00434
00435
00436
00437
00438
00439 }
00440 if ( ionloss < 0.0) ionloss = 0.0 ;
00441
00442
00443 G4double aa = atomicWeight[iMolecula];
00444 ionloss /= (HeEffChargeSquare(0.5*aa, T)*aa);
00445
00446
00447 } else if(1 == (material->GetNumberOfElements())) {
00448 G4double z = material->GetZ() ;
00449 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
00450 }
00451
00452 return ionloss;
00453 }
00454
00455
00456
00457 G4double G4BraggIonModel::ElectronicStoppingPower(G4double z,
00458 G4double kineticEnergy) const
00459 {
00460 G4double ionloss ;
00461 G4int i = G4int(z)-1 ;
00462 if(i < 0) i = 0 ;
00463 if(i > 91) i = 91 ;
00464
00465
00466
00467
00468
00469
00470 G4double T = kineticEnergy*rateMassHe2p/MeV ;
00471
00472 static G4double a[92][5] = {
00473 {0.35485, 0.6456, 6.01525, 20.8933, 4.3515
00474 },{ 0.58, 0.59, 6.3, 130.0, 44.07
00475 },{ 1.42, 0.49, 12.25, 32.0, 9.161
00476 },{ 2.206, 0.51, 15.32, 0.25, 8.995
00477
00478 },{ 3.691, 0.4128, 18.48, 50.72, 9.0
00479 },{ 3.83523, 0.42993,12.6125, 227.41, 188.97
00480 },{ 1.9259, 0.5550, 27.15125, 26.0665, 6.2768
00481 },{ 2.81015, 0.4759, 50.0253, 10.556, 1.0382
00482 },{ 1.533, 0.531, 40.44, 18.41, 2.718
00483 },{ 2.303, 0.4861, 37.01, 37.96, 5.092
00484
00485 },{ 9.894, 0.3081, 23.65, 0.384, 92.93
00486 },{ 4.3, 0.47, 34.3, 3.3, 12.74
00487 },{ 2.5, 0.625, 45.7, 0.1, 4.359
00488 },{ 2.1, 0.65, 49.34, 1.788, 4.133
00489 },{ 1.729, 0.6562, 53.41, 2.405, 3.845
00490 },{ 1.402, 0.6791, 58.98, 3.528, 3.211
00491 },{ 1.117, 0.7044, 69.69, 3.705, 2.156
00492 },{ 2.291, 0.6284, 73.88, 4.478, 2.066
00493 },{ 8.554, 0.3817, 83.61, 11.84, 1.875
00494 },{ 6.297, 0.4622, 65.39, 10.14, 5.036
00495
00496 },{ 5.307, 0.4918, 61.74, 12.4, 6.665
00497 },{ 4.71, 0.5087, 65.28, 8.806, 5.948
00498 },{ 6.151, 0.4524, 83.0, 18.31, 2.71
00499 },{ 6.57, 0.4322, 84.76, 15.53, 2.779
00500 },{ 5.738, 0.4492, 84.6, 14.18, 3.101
00501 },{ 5.013, 0.4707, 85.8, 16.55, 3.211
00502 },{ 4.32, 0.4947, 76.14, 10.85, 5.441
00503 },{ 4.652, 0.4571, 80.73, 22.0, 4.952
00504 },{ 3.114, 0.5236, 76.67, 7.62, 6.385
00505 },{ 3.114, 0.5236, 76.67, 7.62, 7.502
00506
00507 },{ 3.114, 0.5236, 76.67, 7.62, 8.514
00508 },{ 5.746, 0.4662, 79.24, 1.185, 7.993
00509 },{ 2.792, 0.6346, 106.1, 0.2986, 2.331
00510 },{ 4.667, 0.5095, 124.3, 2.102, 1.667
00511 },{ 2.44, 0.6346, 105.0, 0.83, 2.851
00512 },{ 1.413, 0.7377, 147.9, 1.466, 1.016
00513 },{ 11.72, 0.3826, 102.8, 9.231, 4.371
00514 },{ 7.126, 0.4804, 119.3, 5.784, 2.454
00515 },{ 11.61, 0.3955, 146.7, 7.031, 1.423
00516 },{ 10.99, 0.41, 163.9, 7.1, 1.052
00517
00518 },{ 9.241, 0.4275, 163.1, 7.954, 1.102
00519 },{ 9.276, 0.418, 157.1, 8.038, 1.29
00520 },{ 3.999, 0.6152, 97.6, 1.297, 5.792
00521 },{ 4.306, 0.5658, 97.99, 5.514, 5.754
00522 },{ 3.615, 0.6197, 86.26, 0.333, 8.689
00523 },{ 5.8, 0.49, 147.2, 6.903, 1.289
00524 },{ 5.6, 0.49, 130.0, 10.0, 2.844
00525 },{ 3.55, 0.6068, 124.7, 1.112, 3.119
00526 },{ 3.6, 0.62, 105.8, 0.1692, 6.026
00527 },{ 5.4, 0.53, 103.1, 3.931, 7.767
00528
00529 },{ 3.97, 0.6459, 131.8, 0.2233, 2.723
00530 },{ 3.65, 0.64, 126.8, 0.6834, 3.411
00531 },{ 3.118, 0.6519, 164.9, 1.208, 1.51
00532 },{ 3.949, 0.6209, 200.5, 1.878, 0.9126
00533 },{ 14.4, 0.3923, 152.5, 8.354, 2.597
00534 },{ 10.99, 0.4599, 138.4, 4.811, 3.726
00535 },{ 16.6, 0.3773, 224.1, 6.28, 0.9121
00536 },{ 10.54, 0.4533, 159.3, 4.832, 2.529
00537 },{ 10.33, 0.4502, 162.0, 5.132, 2.444
00538 },{ 10.15, 0.4471, 165.6, 5.378, 2.328
00539
00540 },{ 9.976, 0.4439, 168.0, 5.721, 2.258
00541 },{ 9.804, 0.4408, 176.2, 5.675, 1.997
00542 },{ 14.22, 0.363, 228.4, 7.024, 1.016
00543 },{ 9.952, 0.4318, 233.5, 5.065, 0.9244
00544 },{ 9.272, 0.4345, 210.0, 4.911, 1.258
00545 },{ 10.13, 0.4146, 225.7, 5.525, 1.055
00546 },{ 8.949, 0.4304, 213.3, 5.071, 1.221
00547 },{ 11.94, 0.3783, 247.2, 6.655, 0.849
00548 },{ 8.472, 0.4405, 195.5, 4.051, 1.604
00549 },{ 8.301, 0.4399, 203.7, 3.667, 1.459
00550
00551 },{ 6.567, 0.4858, 193.0, 2.65, 1.66
00552 },{ 5.951, 0.5016, 196.1, 2.662, 1.589
00553 },{ 7.495, 0.4523, 251.4, 3.433, 0.8619
00554 },{ 6.335, 0.4825, 255.1, 2.834, 0.8228
00555 },{ 4.314, 0.5558, 214.8, 2.354, 1.263
00556 },{ 4.02, 0.5681, 219.9, 2.402, 1.191
00557 },{ 3.836, 0.5765, 210.2, 2.742, 1.305
00558 },{ 4.68, 0.5247, 244.7, 2.749, 0.8962
00559 },{ 2.892, 0.6204, 208.6, 2.415, 1.416
00560
00561 },{ 2.892, 0.6204, 208.6, 2.415, 1.416
00562
00563 },{ 4.728, 0.5522, 217.0, 3.091, 1.386
00564 },{ 6.18, 0.52, 170.0, 4.0, 3.224
00565 },{ 9.0, 0.47, 198.0, 3.8, 2.032
00566 },{ 2.324, 0.6997, 216.0, 1.599, 1.399
00567 },{ 1.961, 0.7286, 223.0, 1.621, 1.296
00568 },{ 1.75, 0.7427, 350.1, 0.9789, 0.5507
00569 },{ 10.31, 0.4613, 261.2, 4.738, 0.9899
00570 },{ 7.962, 0.519, 235.7, 4.347, 1.313
00571 },{ 6.227, 0.5645, 231.9, 3.961, 1.379
00572 },{ 5.246, 0.5947, 228.6, 4.027, 1.432
00573
00574 },{ 5.408, 0.5811, 235.7, 3.961, 1.358
00575 },{ 5.218, 0.5828, 245.0, 3.838, 1.25}
00576 };
00577
00578
00579 if ( T < 0.001 ) {
00580 G4double slow = a[i][0] ;
00581 G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
00582 * a[i][2]*1000.0 ;
00583 ionloss = slow*shigh / (slow + shigh) ;
00584 ionloss *= sqrt(T*1000.0) ;
00585
00586
00587 } else {
00588 G4double slow = a[i][0] * pow((T*1000.0), a[i][1]) ;
00589 G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
00590 ionloss = slow*shigh / (slow + shigh) ;
00591
00592
00593
00594
00595
00596
00597
00598 }
00599 if ( ionloss < 0.0) { ionloss = 0.0; }
00600
00601
00602 ionloss /= HeEffChargeSquare(z, T);
00603
00604 return ionloss;
00605 }
00606
00607
00608
00609 G4double G4BraggIonModel::DEDX(const G4Material* material,
00610 G4double kineticEnergy)
00611 {
00612 G4double eloss = 0.0;
00613
00614 if(material != currentMaterial) {
00615 currentMaterial = material;
00616 iASTAR = -1;
00617 iMolecula = -1;
00618 if( !HasMaterial(material) ) { iASTAR = astar.GetIndex(material); }
00619 }
00620
00621 const G4int numberOfElements = material->GetNumberOfElements();
00622 const G4double* theAtomicNumDensityVector =
00623 material->GetAtomicNumDensityVector();
00624
00625 if( iASTAR >= 0 ) {
00626 G4double T = kineticEnergy*rateMassHe2p;
00627 return astar.GetElectronicDEDX(iASTAR, T)*material->GetDensity()/
00628 HeEffChargeSquare(astar.GetEffectiveZ(iASTAR), T/MeV);
00629
00630 } else if(iMolecula >= 0) {
00631
00632 eloss = StoppingPower(material, kineticEnergy)*
00633 material->GetDensity()/amu;
00634
00635
00636 } else if(1 == numberOfElements) {
00637
00638 G4double z = material->GetZ();
00639 eloss = ElectronicStoppingPower(z, kineticEnergy)
00640 * (material->GetTotNbOfAtomsPerVolume());
00641
00642
00643 } else {
00644 const G4ElementVector* theElementVector =
00645 material->GetElementVector() ;
00646
00647
00648 for (G4int i=0; i<numberOfElements; i++)
00649 {
00650 const G4Element* element = (*theElementVector)[i] ;
00651 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
00652 * theAtomicNumDensityVector[i];
00653 }
00654 }
00655 return eloss*theZieglerFactor;
00656 }
00657
00658
00659
00660 G4double G4BraggIonModel::HeEffChargeSquare(G4double z,
00661 G4double kinEnergyHeInMeV) const
00662 {
00663
00664
00665
00666
00667
00668 static G4double c[6] = {0.2865, 0.1266, -0.001429,
00669 0.02402,-0.01135, 0.001475};
00670
00671 G4double e = std::max(0.0,std::log(kinEnergyHeInMeV*massFactor));
00672 G4double x = c[0] ;
00673 G4double y = 1.0 ;
00674 for (G4int i=1; i<6; i++) {
00675 y *= e ;
00676 x += y * c[i] ;
00677 }
00678
00679 G4double w = 7.6 - e ;
00680 w = 1.0 + (0.007 + 0.00005*z) * exp( -w*w ) ;
00681 w = 4.0 * (1.0 - exp(-x)) * w * w ;
00682
00683 return w;
00684 }
00685
00686
00687