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 #include "G4MuBetheBlochModel.hh"
00060 #include "G4PhysicalConstants.hh"
00061 #include "G4SystemOfUnits.hh"
00062 #include "Randomize.hh"
00063 #include "G4Electron.hh"
00064 #include "G4LossTableManager.hh"
00065 #include "G4EmCorrections.hh"
00066 #include "G4ParticleChangeForLoss.hh"
00067
00068 G4double G4MuBetheBlochModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083, 0.5917,
00069 0.7628, 0.8983, 0.9801 };
00070
00071 G4double G4MuBetheBlochModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813, 0.1813,
00072 0.1569, 0.1112, 0.0506 };
00073
00074
00075
00076 using namespace std;
00077
00078 G4MuBetheBlochModel::G4MuBetheBlochModel(const G4ParticleDefinition* p,
00079 const G4String& nam)
00080 : G4VEmModel(nam),
00081 particle(0),
00082 limitKinEnergy(100.*keV),
00083 logLimitKinEnergy(log(limitKinEnergy)),
00084 twoln10(2.0*log(10.0)),
00085 bg2lim(0.0169),
00086 taulim(8.4146e-3),
00087 alphaprime(fine_structure_const/twopi)
00088 {
00089 theElectron = G4Electron::Electron();
00090 corr = G4LossTableManager::Instance()->EmCorrections();
00091 fParticleChange = 0;
00092
00093
00094
00095 mass = massSquare = ratio = 1.0;
00096
00097 if(p) { SetParticle(p); }
00098 }
00099
00100
00101
00102 G4MuBetheBlochModel::~G4MuBetheBlochModel()
00103 {}
00104
00105
00106
00107 G4double G4MuBetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
00108 const G4MaterialCutsCouple* couple)
00109 {
00110 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
00111 }
00112
00113
00114
00115 G4double G4MuBetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
00116 G4double kinEnergy)
00117 {
00118 G4double tau = kinEnergy/mass;
00119 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
00120 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
00121 return tmax;
00122 }
00123
00124
00125
00126 void G4MuBetheBlochModel::Initialise(const G4ParticleDefinition* p,
00127 const G4DataVector&)
00128 {
00129 if(p) { SetParticle(p); }
00130 if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
00131 }
00132
00133
00134
00135 G4double G4MuBetheBlochModel::ComputeCrossSectionPerElectron(
00136 const G4ParticleDefinition* p,
00137 G4double kineticEnergy,
00138 G4double cutEnergy,
00139 G4double maxKinEnergy)
00140 {
00141 G4double cross = 0.0;
00142 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00143 G4double maxEnergy = min(tmax,maxKinEnergy);
00144 if(cutEnergy < maxEnergy) {
00145
00146 G4double totEnergy = kineticEnergy + mass;
00147 G4double energy2 = totEnergy*totEnergy;
00148 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
00149
00150 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax
00151 + 0.5*(maxEnergy - cutEnergy)/energy2;
00152
00153
00154 if (maxEnergy > limitKinEnergy) {
00155
00156 G4double logtmax = log(maxEnergy);
00157 G4double logtmin = log(max(cutEnergy,limitKinEnergy));
00158 G4double logstep = logtmax - logtmin;
00159 G4double dcross = 0.0;
00160
00161 for (G4int ll=0; ll<8; ll++)
00162 {
00163 G4double ep = exp(logtmin + xgi[ll]*logstep);
00164 G4double a1 = log(1.0 + 2.0*ep/electron_mass_c2);
00165 G4double a3 = log(4.0*totEnergy*(totEnergy - ep)/massSquare);
00166 dcross += wgi[ll]*(1.0/ep - beta2/tmax + 0.5*ep/energy2)*a1*(a3 - a1);
00167 }
00168
00169 cross += dcross*logstep*alphaprime;
00170 }
00171
00172 cross *= twopi_mc2_rcl2/beta2;
00173
00174 }
00175
00176
00177
00178
00179 return cross;
00180 }
00181
00182
00183
00184 G4double G4MuBetheBlochModel::ComputeCrossSectionPerAtom(
00185 const G4ParticleDefinition* p,
00186 G4double kineticEnergy,
00187 G4double Z, G4double,
00188 G4double cutEnergy,
00189 G4double maxEnergy)
00190 {
00191 G4double cross = Z*ComputeCrossSectionPerElectron
00192 (p,kineticEnergy,cutEnergy,maxEnergy);
00193 return cross;
00194 }
00195
00196
00197
00198 G4double G4MuBetheBlochModel::CrossSectionPerVolume(
00199 const G4Material* material,
00200 const G4ParticleDefinition* p,
00201 G4double kineticEnergy,
00202 G4double cutEnergy,
00203 G4double maxEnergy)
00204 {
00205 G4double eDensity = material->GetElectronDensity();
00206 G4double cross = eDensity*ComputeCrossSectionPerElectron
00207 (p,kineticEnergy,cutEnergy,maxEnergy);
00208 return cross;
00209 }
00210
00211
00212
00213 G4double G4MuBetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
00214 const G4ParticleDefinition* p,
00215 G4double kineticEnergy,
00216 G4double cut)
00217 {
00218 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00219 G4double tau = kineticEnergy/mass;
00220 G4double cutEnergy = min(cut,tmax);
00221 G4double gam = tau + 1.0;
00222 G4double bg2 = tau * (tau+2.0);
00223 G4double beta2 = bg2/(gam*gam);
00224
00225 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
00226 G4double eexc2 = eexc*eexc;
00227
00228
00229
00230
00231
00232
00233 G4double eDensity = material->GetElectronDensity();
00234
00235 G4double dedx = log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
00236 -(1.0 + cutEnergy/tmax)*beta2;
00237
00238 G4double totEnergy = kineticEnergy + mass;
00239 G4double del = 0.5*cutEnergy/totEnergy;
00240 dedx += del*del;
00241
00242
00243 G4double x = log(bg2)/twoln10;
00244
00245
00246
00247
00248 dedx -= material->GetIonisation()->DensityCorrection(x);
00249
00250
00251 dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
00252
00253
00254
00255 if (dedx < 0.0) dedx = 0.0 ;
00256
00257
00258 if (cutEnergy > limitKinEnergy) {
00259
00260 G4double logtmax = log(cutEnergy);
00261 G4double logstep = logtmax - logLimitKinEnergy;
00262 G4double dloss = 0.0;
00263 G4double ftot2= 0.5/(totEnergy*totEnergy);
00264
00265 for (G4int ll=0; ll<8; ll++)
00266 {
00267 G4double ep = exp(logLimitKinEnergy + xgi[ll]*logstep);
00268 G4double a1 = log(1.0 + 2.0*ep/electron_mass_c2);
00269 G4double a3 = log(4.0*totEnergy*(totEnergy - ep)/massSquare);
00270 dloss += wgi[ll]*(1.0 - beta2*ep/tmax + ep*ep*ftot2)*a1*(a3 - a1);
00271 }
00272 dedx += dloss*logstep*alphaprime;
00273 }
00274
00275 dedx *= twopi_mc2_rcl2*eDensity/beta2;
00276
00277
00278 dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
00279
00280 return dedx;
00281 }
00282
00283
00284
00285 void G4MuBetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
00286 const G4MaterialCutsCouple*,
00287 const G4DynamicParticle* dp,
00288 G4double minKinEnergy,
00289 G4double maxEnergy)
00290 {
00291 G4double tmax = MaxSecondaryKinEnergy(dp);
00292 G4double maxKinEnergy = min(maxEnergy,tmax);
00293 if(minKinEnergy >= maxKinEnergy) { return; }
00294
00295 G4double kineticEnergy = dp->GetKineticEnergy();
00296 G4double totEnergy = kineticEnergy + mass;
00297 G4double etot2 = totEnergy*totEnergy;
00298 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
00299
00300 G4double grej = 1.;
00301 if(tmax > limitKinEnergy) {
00302 G4double a0 = log(2.*totEnergy/mass);
00303 grej += alphaprime*a0*a0;
00304 }
00305
00306 G4double deltaKinEnergy, f;
00307
00308
00309 do {
00310 G4double q = G4UniformRand();
00311 deltaKinEnergy = minKinEnergy*maxKinEnergy
00312 /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
00313
00314
00315 f = 1.0 - beta2*deltaKinEnergy/tmax
00316 + 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
00317
00318 if(deltaKinEnergy > limitKinEnergy) {
00319 G4double a1 = log(1.0 + 2.0*deltaKinEnergy/electron_mass_c2);
00320 G4double a3 = log(4.0*totEnergy*(totEnergy - deltaKinEnergy)/massSquare);
00321 f *= (1. + alphaprime*a1*(a3 - a1));
00322 }
00323
00324 if(f > grej) {
00325 G4cout << "G4MuBetheBlochModel::SampleSecondary Warning! "
00326 << "Majorant " << grej << " < "
00327 << f << " for edelta= " << deltaKinEnergy
00328 << " tmin= " << minKinEnergy << " max= " << maxKinEnergy
00329 << G4endl;
00330 }
00331
00332
00333 } while( grej*G4UniformRand() > f );
00334
00335 G4double deltaMomentum =
00336 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
00337 G4double totalMomentum = totEnergy*sqrt(beta2);
00338 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
00339 (deltaMomentum * totalMomentum);
00340
00341 G4double sint = sqrt(1.0 - cost*cost);
00342
00343 G4double phi = twopi * G4UniformRand() ;
00344
00345 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
00346 G4ThreeVector direction = dp->GetMomentumDirection();
00347 deltaDirection.rotateUz(direction);
00348
00349
00350 kineticEnergy -= deltaKinEnergy;
00351 G4ThreeVector dir = totalMomentum*direction - deltaMomentum*deltaDirection;
00352 direction = dir.unit();
00353 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00354 fParticleChange->SetProposedMomentumDirection(direction);
00355
00356
00357 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
00358 deltaDirection,deltaKinEnergy);
00359 vdp->push_back(delta);
00360 }
00361
00362