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 #include "G4BetheBlochModel.hh"
00061 #include "Randomize.hh"
00062 #include "G4PhysicalConstants.hh"
00063 #include "G4SystemOfUnits.hh"
00064 #include "G4Electron.hh"
00065 #include "G4LossTableManager.hh"
00066 #include "G4EmCorrections.hh"
00067 #include "G4ParticleChangeForLoss.hh"
00068
00069
00070
00071 using namespace std;
00072
00073 G4BetheBlochModel::G4BetheBlochModel(const G4ParticleDefinition* p,
00074 const G4String& nam)
00075 : G4VEmModel(nam),
00076 particle(0),
00077 tlimit(DBL_MAX),
00078 twoln10(2.0*log(10.0)),
00079 bg2lim(0.0169),
00080 taulim(8.4146e-3),
00081 isIon(false),
00082 isInitialised(false)
00083 {
00084 fParticleChange = 0;
00085 theElectron = G4Electron::Electron();
00086 if(p) {
00087 SetGenericIon(p);
00088 SetParticle(p);
00089 } else {
00090 SetParticle(theElectron);
00091 }
00092 corr = G4LossTableManager::Instance()->EmCorrections();
00093 nist = G4NistManager::Instance();
00094 SetLowEnergyLimit(2.0*MeV);
00095 }
00096
00097
00098
00099 G4BetheBlochModel::~G4BetheBlochModel()
00100 {}
00101
00102
00103
00104 void G4BetheBlochModel::Initialise(const G4ParticleDefinition* p,
00105 const G4DataVector&)
00106 {
00107 SetGenericIon(p);
00108 SetParticle(p);
00109
00110
00111
00112
00113
00114
00115 SetDeexcitationFlag(false);
00116
00117 if(!isInitialised) {
00118 isInitialised = true;
00119 fParticleChange = GetParticleChangeForLoss();
00120 }
00121 }
00122
00123
00124
00125 G4double G4BetheBlochModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
00126 const G4Material* mat,
00127 G4double kineticEnergy)
00128 {
00129
00130 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
00131 corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
00132 return corrFactor;
00133 }
00134
00135
00136
00137 G4double G4BetheBlochModel::GetParticleCharge(const G4ParticleDefinition* p,
00138 const G4Material* mat,
00139 G4double kineticEnergy)
00140 {
00141
00142 return corr->GetParticleCharge(p,mat,kineticEnergy);
00143 }
00144
00145
00146
00147 void G4BetheBlochModel::SetupParameters()
00148 {
00149 mass = particle->GetPDGMass();
00150 spin = particle->GetPDGSpin();
00151 G4double q = particle->GetPDGCharge()/eplus;
00152 chargeSquare = q*q;
00153 corrFactor = chargeSquare;
00154 ratio = electron_mass_c2/mass;
00155 G4double magmom =
00156 particle->GetPDGMagneticMoment()*mass/(0.5*eplus*hbar_Planck*c_squared);
00157 magMoment2 = magmom*magmom - 1.0;
00158 formfact = 0.0;
00159 if(particle->GetLeptonNumber() == 0) {
00160 G4double x = 0.8426*GeV;
00161 if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
00162 else if(mass > GeV) {
00163 x /= nist->GetZ13(mass/proton_mass_c2);
00164
00165 }
00166 formfact = 2.0*electron_mass_c2/(x*x);
00167 tlimit = 2.0/formfact;
00168 }
00169 }
00170
00171
00172
00173 G4double
00174 G4BetheBlochModel::ComputeCrossSectionPerElectron(const G4ParticleDefinition* p,
00175 G4double kineticEnergy,
00176 G4double cutEnergy,
00177 G4double maxKinEnergy)
00178 {
00179 G4double cross = 0.0;
00180 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00181 G4double maxEnergy = min(tmax,maxKinEnergy);
00182 if(cutEnergy < maxEnergy) {
00183
00184 G4double totEnergy = kineticEnergy + mass;
00185 G4double energy2 = totEnergy*totEnergy;
00186 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
00187
00188 cross = 1.0/cutEnergy - 1.0/maxEnergy
00189 - beta2*log(maxEnergy/cutEnergy)/tmax;
00190
00191
00192 if( 0.5 == spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
00193
00194
00195
00196
00197
00198
00199
00200 cross *= twopi_mc2_rcl2*chargeSquare/beta2;
00201 }
00202
00203
00204
00205
00206 return cross;
00207 }
00208
00209
00210
00211 G4double G4BetheBlochModel::ComputeCrossSectionPerAtom(
00212 const G4ParticleDefinition* p,
00213 G4double kineticEnergy,
00214 G4double Z, G4double,
00215 G4double cutEnergy,
00216 G4double maxEnergy)
00217 {
00218 G4double cross = Z*ComputeCrossSectionPerElectron
00219 (p,kineticEnergy,cutEnergy,maxEnergy);
00220 return cross;
00221 }
00222
00223
00224
00225 G4double G4BetheBlochModel::CrossSectionPerVolume(
00226 const G4Material* material,
00227 const G4ParticleDefinition* p,
00228 G4double kineticEnergy,
00229 G4double cutEnergy,
00230 G4double maxEnergy)
00231 {
00232 G4double eDensity = material->GetElectronDensity();
00233 G4double cross = eDensity*ComputeCrossSectionPerElectron
00234 (p,kineticEnergy,cutEnergy,maxEnergy);
00235 return cross;
00236 }
00237
00238
00239
00240 G4double G4BetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
00241 const G4ParticleDefinition* p,
00242 G4double kineticEnergy,
00243 G4double cut)
00244 {
00245 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00246 G4double cutEnergy = std::min(cut,tmax);
00247
00248 G4double tau = kineticEnergy/mass;
00249 G4double gam = tau + 1.0;
00250 G4double bg2 = tau * (tau+2.0);
00251 G4double beta2 = bg2/(gam*gam);
00252
00253 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
00254 G4double eexc2 = eexc*eexc;
00255
00256 G4double eDensity = material->GetElectronDensity();
00257
00258 G4double dedx = log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
00259 - (1.0 + cutEnergy/tmax)*beta2;
00260
00261 if(0.5 == spin) {
00262 G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
00263 dedx += del*del;
00264 }
00265
00266
00267 G4double x = log(bg2)/twoln10;
00268 dedx -= material->GetIonisation()->DensityCorrection(x);
00269
00270
00271 dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
00272
00273
00274 dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
00275
00276
00277 if(isIon) {
00278 dedx += corr->IonBarkasCorrection(p,material,kineticEnergy);
00279 } else {
00280 dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
00281 }
00282
00283 if (dedx < 0.0) { dedx = 0.0; }
00284
00285
00286
00287
00288 return dedx;
00289 }
00290
00291
00292
00293 void G4BetheBlochModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
00294 const G4DynamicParticle* dp,
00295 G4double& eloss,
00296 G4double&,
00297 G4double length)
00298 {
00299 if(isIon) {
00300 const G4ParticleDefinition* p = dp->GetDefinition();
00301 const G4Material* mat = couple->GetMaterial();
00302 G4double preKinEnergy = dp->GetKineticEnergy();
00303 G4double e = preKinEnergy - eloss*0.5;
00304 if(e < preKinEnergy*0.75) { e = preKinEnergy*0.75; }
00305
00306 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
00307 GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
00308 G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
00309 G4double highOrder = length*corr->IonHighOrderCorrections(p,couple,e);
00310 G4double elossnew = eloss*qfactor + highOrder;
00311 if(elossnew > preKinEnergy) { elossnew = preKinEnergy; }
00312 else if(elossnew < eloss*0.5) { elossnew = eloss*0.5; }
00313 eloss = elossnew;
00314
00315
00316
00317 }
00318 }
00319
00320
00321
00322 void G4BetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
00323 const G4MaterialCutsCouple*,
00324 const G4DynamicParticle* dp,
00325 G4double minKinEnergy,
00326 G4double maxEnergy)
00327 {
00328 G4double kineticEnergy = dp->GetKineticEnergy();
00329 G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
00330
00331 G4double maxKinEnergy = std::min(maxEnergy,tmax);
00332 if(minKinEnergy >= maxKinEnergy) { return; }
00333
00334 G4double totEnergy = kineticEnergy + mass;
00335 G4double etot2 = totEnergy*totEnergy;
00336 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
00337
00338 G4double deltaKinEnergy, f;
00339 G4double f1 = 0.0;
00340 G4double fmax = 1.0;
00341 if( 0.5 == spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
00342
00343
00344 do {
00345 G4double q = G4UniformRand();
00346 deltaKinEnergy = minKinEnergy*maxKinEnergy
00347 /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
00348
00349 f = 1.0 - beta2*deltaKinEnergy/tmax;
00350 if( 0.5 == spin ) {
00351 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
00352 f += f1;
00353 }
00354
00355 } while( fmax*G4UniformRand() > f);
00356
00357
00358
00359
00360 G4double x = formfact*deltaKinEnergy;
00361 if(x > 1.e-6) {
00362
00363 G4double x1 = 1.0 + x;
00364 G4double grej = 1.0/(x1*x1);
00365 if( 0.5 == spin ) {
00366 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
00367 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
00368 }
00369 if(grej > 1.1) {
00370 G4cout << "### G4BetheBlochModel WARNING: grej= " << grej
00371 << " " << dp->GetDefinition()->GetParticleName()
00372 << " Ekin(MeV)= " << kineticEnergy
00373 << " delEkin(MeV)= " << deltaKinEnergy
00374 << G4endl;
00375 }
00376 if(G4UniformRand() > grej) return;
00377 }
00378
00379
00380 G4double totMomentum = totEnergy*sqrt(beta2);
00381 G4double deltaMomentum =
00382 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
00383 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
00384 (deltaMomentum * totMomentum);
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
00402
00403 G4double phi = twopi * G4UniformRand() ;
00404
00405
00406 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost);
00407 G4ThreeVector direction = dp->GetMomentumDirection();
00408 deltaDirection.rotateUz(direction);
00409
00410
00411 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
00412 deltaDirection,deltaKinEnergy);
00413
00414 vdp->push_back(delta);
00415
00416
00417 kineticEnergy -= deltaKinEnergy;
00418 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
00419 finalP = finalP.unit();
00420
00421 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00422 fParticleChange->SetProposedMomentumDirection(finalP);
00423 }
00424
00425
00426
00427 G4double G4BetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
00428 G4double kinEnergy)
00429 {
00430
00431 SetParticle(pd);
00432 G4double tau = kinEnergy/mass;
00433 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
00434 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
00435 return std::min(tmax,tlimit);
00436 }
00437
00438