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 #include "G4BraggModel.hh"
00067 #include "G4PhysicalConstants.hh"
00068 #include "G4SystemOfUnits.hh"
00069 #include "Randomize.hh"
00070 #include "G4Electron.hh"
00071 #include "G4ParticleChangeForLoss.hh"
00072 #include "G4LossTableManager.hh"
00073 #include "G4EmCorrections.hh"
00074
00075
00076
00077 using namespace std;
00078
00079 G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam)
00080 : G4VEmModel(nam),
00081 particle(0),
00082 currentMaterial(0),
00083 protonMassAMU(1.007276),
00084 iMolecula(-1),
00085 iPSTAR(-1),
00086 isIon(false),
00087 isInitialised(false)
00088 {
00089 fParticleChange = 0;
00090 SetHighEnergyLimit(2.0*MeV);
00091
00092 lowestKinEnergy = 1.0*keV;
00093 theZieglerFactor = eV*cm2*1.0e-15;
00094 theElectron = G4Electron::Electron();
00095 expStopPower125 = 0.0;
00096
00097 corr = G4LossTableManager::Instance()->EmCorrections();
00098 if(p) { SetParticle(p); }
00099 else { SetParticle(theElectron); }
00100 }
00101
00102
00103
00104 G4BraggModel::~G4BraggModel()
00105 {}
00106
00107
00108
00109 void G4BraggModel::Initialise(const G4ParticleDefinition* p,
00110 const G4DataVector&)
00111 {
00112 if(p != particle) { SetParticle(p); }
00113
00114
00115 SetDeexcitationFlag(false);
00116
00117 if(!isInitialised) {
00118 isInitialised = true;
00119
00120 G4String pname = particle->GetParticleName();
00121 if(particle->GetParticleType() == "nucleus" &&
00122 pname != "deuteron" && pname != "triton" &&
00123 pname != "alpha+" && pname != "helium" &&
00124 pname != "hydrogen") { isIon = true; }
00125
00126 fParticleChange = GetParticleChangeForLoss();
00127 }
00128 }
00129
00130
00131
00132 G4double G4BraggModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
00133 const G4Material* mat,
00134 G4double kineticEnergy)
00135 {
00136
00137 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
00138 GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
00139 return q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
00140 }
00141
00142
00143
00144 G4double G4BraggModel::GetParticleCharge(const G4ParticleDefinition* p,
00145 const G4Material* mat,
00146 G4double kineticEnergy)
00147 {
00148
00149 return corr->GetParticleCharge(p,mat,kineticEnergy);
00150 }
00151
00152
00153
00154 G4double G4BraggModel::ComputeCrossSectionPerElectron(
00155 const G4ParticleDefinition* p,
00156 G4double kineticEnergy,
00157 G4double cutEnergy,
00158 G4double maxKinEnergy)
00159 {
00160 G4double cross = 0.0;
00161 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00162 G4double maxEnergy = std::min(tmax,maxKinEnergy);
00163 if(cutEnergy < maxEnergy) {
00164
00165 G4double energy = kineticEnergy + mass;
00166 G4double energy2 = energy*energy;
00167 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
00168 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax;
00169
00170 cross *= twopi_mc2_rcl2*chargeSquare/beta2;
00171 }
00172
00173
00174
00175 return cross;
00176 }
00177
00178
00179
00180 G4double G4BraggModel::ComputeCrossSectionPerAtom(
00181 const G4ParticleDefinition* p,
00182 G4double kineticEnergy,
00183 G4double Z, G4double,
00184 G4double cutEnergy,
00185 G4double maxEnergy)
00186 {
00187 G4double cross = Z*ComputeCrossSectionPerElectron
00188 (p,kineticEnergy,cutEnergy,maxEnergy);
00189 return cross;
00190 }
00191
00192
00193
00194 G4double G4BraggModel::CrossSectionPerVolume(
00195 const G4Material* material,
00196 const G4ParticleDefinition* p,
00197 G4double kineticEnergy,
00198 G4double cutEnergy,
00199 G4double maxEnergy)
00200 {
00201 G4double eDensity = material->GetElectronDensity();
00202 G4double cross = eDensity*ComputeCrossSectionPerElectron
00203 (p,kineticEnergy,cutEnergy,maxEnergy);
00204 return cross;
00205 }
00206
00207
00208
00209 G4double G4BraggModel::ComputeDEDXPerVolume(const G4Material* material,
00210 const G4ParticleDefinition* p,
00211 G4double kineticEnergy,
00212 G4double cutEnergy)
00213 {
00214 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00215 G4double tkin = kineticEnergy/massRate;
00216 G4double dedx = 0.0;
00217
00218 if(tkin < lowestKinEnergy) {
00219 dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
00220 } else {
00221 dedx = DEDX(material, tkin);
00222 }
00223
00224 if (cutEnergy < tmax) {
00225
00226 G4double tau = kineticEnergy/mass;
00227 G4double gam = tau + 1.0;
00228 G4double bg2 = tau * (tau+2.0);
00229 G4double beta2 = bg2/(gam*gam);
00230 G4double x = cutEnergy/tmax;
00231
00232 dedx += (log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
00233 * (material->GetElectronDensity())/beta2;
00234 }
00235
00236
00237
00238 if (dedx < 0.0) { dedx = 0.0; }
00239
00240 dedx *= chargeSquare;
00241
00242
00243
00244
00245 return dedx;
00246 }
00247
00248
00249
00250 void G4BraggModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
00251 const G4MaterialCutsCouple*,
00252 const G4DynamicParticle* dp,
00253 G4double xmin,
00254 G4double maxEnergy)
00255 {
00256 G4double tmax = MaxSecondaryKinEnergy(dp);
00257 G4double xmax = std::min(tmax, maxEnergy);
00258 if(xmin >= xmax) { return; }
00259
00260 G4double kineticEnergy = dp->GetKineticEnergy();
00261 G4double energy = kineticEnergy + mass;
00262 G4double energy2 = energy*energy;
00263 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
00264 G4double grej = 1.0;
00265 G4double deltaKinEnergy, f;
00266
00267 G4ThreeVector direction = dp->GetMomentumDirection();
00268
00269
00270 do {
00271 G4double q = G4UniformRand();
00272 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
00273
00274 f = 1.0 - beta2*deltaKinEnergy/tmax;
00275
00276 if(f > grej) {
00277 G4cout << "G4BraggModel::SampleSecondary Warning! "
00278 << "Majorant " << grej << " < "
00279 << f << " for e= " << deltaKinEnergy
00280 << G4endl;
00281 }
00282
00283 } while( grej*G4UniformRand() >= f );
00284
00285 G4double deltaMomentum =
00286 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
00287 G4double totMomentum = energy*sqrt(beta2);
00288 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
00289 (deltaMomentum * totMomentum);
00290 if(cost > 1.0) cost = 1.0;
00291 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
00292
00293 G4double phi = twopi * G4UniformRand() ;
00294
00295 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
00296 deltaDirection.rotateUz(direction);
00297
00298
00299 kineticEnergy -= deltaKinEnergy;
00300 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
00301 finalP = finalP.unit();
00302
00303 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00304 fParticleChange->SetProposedMomentumDirection(finalP);
00305
00306
00307 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,
00308 deltaKinEnergy);
00309
00310 vdp->push_back(delta);
00311 }
00312
00313
00314
00315 G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
00316 G4double kinEnergy)
00317 {
00318 if(pd != particle) { SetParticle(pd); }
00319 G4double tau = kinEnergy/mass;
00320 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
00321 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
00322 return tmax;
00323 }
00324
00325
00326
00327 G4bool G4BraggModel::HasMaterial(const G4Material* material)
00328 {
00329 G4String chFormula = material->GetChemicalFormula();
00330 if("" == chFormula) { return false; }
00331
00332
00333 const size_t numberOfMolecula = 11;
00334 const G4String molName[numberOfMolecula] = {
00335 "Al_2O_3", "CO_2", "CH_4",
00336 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene", "(C_8H_8)_N",
00337 "C_3H_8", "SiO_2", "H_2O",
00338 "H_2O-Gas", "Graphite" } ;
00339 const G4int idxPSTAR[numberOfMolecula] = {
00340 6, 16, 36,
00341 52, 55, 56,
00342 59, 62, 71,
00343 72, 13};
00344
00345
00346 const G4State theState = material->GetState() ;
00347
00348 if( theState == kStateGas && "H_2O" == chFormula) {
00349 chFormula = G4String("H_2O-Gas");
00350 }
00351
00352
00353 for (size_t i=0; i<numberOfMolecula; ++i) {
00354 if (chFormula == molName[i]) {
00355 iMolecula = -1;
00356 iPSTAR = idxPSTAR[i];
00357 return true;
00358 }
00359 }
00360 return false;
00361 }
00362
00363
00364
00365 G4double G4BraggModel::StoppingPower(const G4Material* material,
00366 G4double kineticEnergy)
00367 {
00368 G4double ionloss = 0.0 ;
00369
00370 if (iMolecula >= 0) {
00371
00372
00373
00374
00375
00376 G4double T = kineticEnergy/(keV*protonMassAMU) ;
00377
00378 static G4double a[11][5] = {
00379 {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2},
00380 {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3},
00381 {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3},
00382 {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3},
00383 {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2},
00384 {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1},
00385 {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2},
00386 {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2},
00387 {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3},
00388 {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3},
00389 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2} };
00390
00391 static G4double atomicWeight[11] = {
00392 101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
00393 104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};
00394
00395 if ( T < 10.0 ) {
00396 ionloss = a[iMolecula][0] * sqrt(T) ;
00397
00398 } else if ( T < 10000.0 ) {
00399 G4double slow = a[iMolecula][1] * pow(T, 0.45) ;
00400 G4double shigh = log( 1.0 + a[iMolecula][3]/T
00401 + a[iMolecula][4]*T ) * a[iMolecula][2]/T ;
00402 ionloss = slow*shigh / (slow + shigh) ;
00403 }
00404
00405 if ( ionloss < 0.0) ionloss = 0.0 ;
00406 if ( 10 == iMolecula ) {
00407 if (T < 100.0) {
00408 ionloss *= (1.0+0.023+0.0066*log10(T));
00409 }
00410 else if (T < 700.0) {
00411 ionloss *=(1.0+0.089-0.0248*log10(T-99.));
00412 }
00413 else if (T < 10000.0) {
00414 ionloss *=(1.0+0.089-0.0248*log10(700.-99.));
00415 }
00416 }
00417 ionloss /= atomicWeight[iMolecula];
00418
00419
00420 } else if(1 == (material->GetNumberOfElements())) {
00421 G4double z = material->GetZ() ;
00422 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
00423 }
00424
00425 return ionloss;
00426 }
00427
00428
00429
00430 G4double G4BraggModel::ElectronicStoppingPower(G4double z,
00431 G4double kineticEnergy) const
00432 {
00433 G4double ionloss ;
00434 G4int i = G4int(z)-1 ;
00435 if(i < 0) i = 0 ;
00436 if(i > 91) i = 91 ;
00437
00438
00439
00440
00441
00442 G4double T = kineticEnergy/(keV*protonMassAMU) ;
00443
00444 static G4double a[92][5] = {
00445 {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1},
00446 {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2},
00447 {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2},
00448 {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2},
00449 {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2},
00450 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2},
00451 {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2},
00452 {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2},
00453 {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2},
00454 {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2},
00455
00456 {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2},
00457 {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2},
00458 {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2},
00459 {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2},
00460 {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2},
00461 {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2},
00462 {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2},
00463 {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2},
00464 {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2},
00465 {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2},
00466
00467 {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3},
00468 {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3},
00469 {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3},
00470 {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3},
00471 {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3},
00472 {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3},
00473 {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3},
00474 {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3},
00475 {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2},
00476 {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3},
00477
00478 {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3},
00479 {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3},
00480 {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3},
00481 {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3},
00482 {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3},
00483 {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3},
00484 {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3},
00485 {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3},
00486 {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3},
00487 {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3},
00488
00489 {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3},
00490 {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3},
00491 {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3},
00492 {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3},
00493 {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3},
00494 {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3},
00495
00496 {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2},
00497 {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3},
00498 {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3},
00499 {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3},
00500
00501 {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3},
00502 {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3},
00503 {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3},
00504 {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3},
00505 {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3},
00506 {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3},
00507 {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3},
00508 {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3},
00509 {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3},
00510 {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3},
00511
00512 {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3},
00513 {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3},
00514 {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3},
00515 {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3},
00516 {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3},
00517 {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3},
00518 {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3},
00519 {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3},
00520 {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3},
00521 {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3},
00522
00523 {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3},
00524 {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3},
00525 {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3},
00526 {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3},
00527 {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3},
00528 {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3},
00529 {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3},
00530 {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3},
00531
00532 {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2},
00533 {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3},
00534
00535 {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3},
00536 {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3},
00537 {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3},
00538 {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3},
00539 {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3},
00540 {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3},
00541 {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3},
00542 {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3},
00543 {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3},
00544 {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3},
00545
00546 {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3},
00547 {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3}
00548 };
00549
00550 G4double fac = 1.0 ;
00551
00552
00553 if ( T < 40.0 && 5 == i) {
00554 fac = sqrt(T/40.0) ;
00555 T = 40.0 ;
00556
00557
00558 } else if ( T < 10.0 ) {
00559 fac = sqrt(T*0.1) ;
00560 T =10.0 ;
00561 }
00562
00563
00564 G4double slow = a[i][1] * pow(T, 0.45) ;
00565 G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
00566 ionloss = slow*shigh*fac / (slow + shigh) ;
00567
00568 if ( ionloss < 0.0) { ionloss = 0.0; }
00569
00570 return ionloss;
00571 }
00572
00573
00574
00575 G4double G4BraggModel::DEDX(const G4Material* material, G4double kineticEnergy)
00576 {
00577 G4double eloss = 0.0;
00578
00579
00580 if(material != currentMaterial) {
00581 currentMaterial = material;
00582 iPSTAR = -1;
00583 iMolecula = -1;
00584 if( !HasMaterial(material) ) { iPSTAR = pstar.GetIndex(material); }
00585
00586
00587
00588
00589 }
00590
00591 const G4int numberOfElements = material->GetNumberOfElements();
00592 const G4double* theAtomicNumDensityVector =
00593 material->GetAtomicNumDensityVector();
00594
00595 if( iPSTAR >= 0 ) {
00596 return pstar.GetElectronicDEDX(iPSTAR, kineticEnergy)*material->GetDensity();
00597
00598 } else if(iMolecula >= 0) {
00599
00600 eloss = StoppingPower(material, kineticEnergy)*
00601 material->GetDensity()/amu;
00602
00603
00604 } else if(1 == numberOfElements) {
00605
00606 G4double z = material->GetZ();
00607 eloss = ElectronicStoppingPower(z, kineticEnergy)
00608 * (material->GetTotNbOfAtomsPerVolume());
00609
00610
00611
00612 } else if( MolecIsInZiegler1988(material) ) {
00613
00614
00615 G4double eloss125 = 0.0 ;
00616 const G4ElementVector* theElementVector =
00617 material->GetElementVector();
00618
00619
00620 for (G4int i=0; i<numberOfElements; i++) {
00621 const G4Element* element = (*theElementVector)[i] ;
00622 G4double z = element->GetZ() ;
00623 eloss += ElectronicStoppingPower(z,kineticEnergy)
00624 * theAtomicNumDensityVector[i] ;
00625 eloss125 += ElectronicStoppingPower(z,125.0*keV)
00626 * theAtomicNumDensityVector[i] ;
00627 }
00628
00629
00630 eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
00631
00632
00633 } else {
00634 const G4ElementVector* theElementVector =
00635 material->GetElementVector() ;
00636
00637
00638 for (G4int i=0; i<numberOfElements; i++)
00639 {
00640 const G4Element* element = (*theElementVector)[i] ;
00641 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
00642 * theAtomicNumDensityVector[i];
00643 }
00644 }
00645 return eloss*theZieglerFactor;
00646 }
00647
00648
00649
00650 G4bool G4BraggModel::MolecIsInZiegler1988(const G4Material* material)
00651 {
00652
00653
00654
00655
00656 G4String myFormula = G4String(" ") ;
00657 const G4String chFormula = material->GetChemicalFormula() ;
00658 if (myFormula == chFormula ) { return false; }
00659
00660
00661
00662
00663
00664
00665
00666 myFormula = G4String("H_2O") ;
00667 const G4State theState = material->GetState() ;
00668 if( theState == kStateGas && myFormula == chFormula) return false ;
00669
00670 const size_t numberOfMolecula = 53 ;
00671
00672
00673 const G4double HeEff = 2.8735 ;
00674
00675 static G4String nameOfMol[numberOfMolecula] = {
00676 "H_2O", "C_2H_4O", "C_3H_6O", "C_2H_2", "C_H_3OH",
00677 "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH_3", "C_14H_10",
00678 "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
00679 "CF_4", "C_6H_8", "C_6H_12", "C_6H_10O", "C_6H_10",
00680 "C_8H_16", "C_5H_10", "C_5H_8", "C_3H_6-Cyclopropane","C_2H_4F_2",
00681 "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_2F_6", "C_2H_6O",
00682 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
00683 "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
00684 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_8H_6", "(CH_2)_N",
00685 "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_3H_6-Propylene", "C_3H_6O",
00686 "C_3H_6S", "C_4H_4S", "C_7H_8"
00687 } ;
00688
00689 static G4double expStopping[numberOfMolecula] = {
00690 66.1, 190.4, 258.7, 42.2, 141.5,
00691 210.9, 279.6, 198.8, 31.0, 267.5,
00692 122.8, 311.4, 260.3, 328.9, 391.3,
00693 206.6, 374.0, 422.0, 432.0, 398.0,
00694 554.0, 353.0, 326.0, 74.6, 220.5,
00695 197.4, 362.0, 170.0, 330.5, 211.3,
00696 262.3, 349.6, 51.3, 187.0, 236.9,
00697 121.9, 35.8, 247.0, 292.6, 268.0,
00698 262.3, 49.0, 398.9, 444.0, 22.91,
00699 68.0, 155.0, 84.0, 74.2, 254.7,
00700 306.8, 324.4, 420.0
00701 } ;
00702
00703 static G4double expCharge[numberOfMolecula] = {
00704 HeEff, HeEff, HeEff, 1.0, HeEff,
00705 HeEff, HeEff, HeEff, 1.0, 1.0,
00706 1.0, HeEff, HeEff, HeEff, HeEff,
00707 HeEff, HeEff, HeEff, HeEff, HeEff,
00708 HeEff, HeEff, HeEff, 1.0, HeEff,
00709 HeEff, HeEff, HeEff, HeEff, HeEff,
00710 HeEff, HeEff, 1.0, HeEff, HeEff,
00711 HeEff, 1.0, HeEff, HeEff, HeEff,
00712 HeEff, 1.0, HeEff, HeEff, 1.0,
00713 1.0, 1.0, 1.0, 1.0, HeEff,
00714 HeEff, HeEff, HeEff
00715 } ;
00716
00717 static G4double numberOfAtomsPerMolecula[numberOfMolecula] = {
00718 3.0, 7.0, 10.0, 4.0, 6.0,
00719 9.0, 12.0, 7.0, 4.0, 24.0,
00720 12.0, 14.0, 10.0, 13.0, 5.0,
00721 5.0, 14.0, 18.0, 17.0, 17.0,
00722 24.0, 15.0, 13.0, 9.0, 8.0,
00723 6.0, 14.0, 8.0, 8.0, 9.0,
00724 10.0, 15.0, 6.0, 7.0, 7.0,
00725 3.0, 5.0, 5.0, 5.0, 5.0,
00726 9.0, 3.0, 16.0, 14.0, 3.0,
00727 9.0, 16.0, 11.0, 9.0, 10.0,
00728 10.0, 9.0, 15.0
00729 } ;
00730
00731
00732 for (size_t i=0; i<numberOfMolecula; i++)
00733 {
00734 if(chFormula == nameOfMol[i]) {
00735 G4double exp125 = expStopping[i] *
00736 (material->GetTotNbOfAtomsPerVolume()) /
00737 (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
00738 SetExpStopPower125(exp125);
00739 return true;
00740 }
00741 }
00742
00743 return false;
00744 }
00745
00746
00747
00748 G4double G4BraggModel::ChemicalFactor(G4double kineticEnergy,
00749 G4double eloss125) const
00750 {
00751
00752
00753
00754
00755 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2 ;
00756 G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2 ;
00757 G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ;
00758 G4double beta = sqrt(1.0 - 1.0/(gamma*gamma)) ;
00759 G4double beta25 = sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
00760 G4double beta125 = sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
00761
00762 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
00763 (1.0 + exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
00764 (1.0 + exp( 1.48 * ( beta/beta25 - 7.0 ) ) ) ;
00765
00766 return factor ;
00767 }
00768
00769
00770