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 #include "G4WentzelOKandVIxSection.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "Randomize.hh"
00053 #include "G4Electron.hh"
00054 #include "G4Positron.hh"
00055 #include "G4Proton.hh"
00056 #include "G4LossTableManager.hh"
00057
00058
00059
00060 G4double G4WentzelOKandVIxSection::ScreenRSquareElec[] = {0.0};
00061 G4double G4WentzelOKandVIxSection::ScreenRSquare[] = {0.0};
00062 G4double G4WentzelOKandVIxSection::FormFactor[] = {0.0};
00063
00064 using namespace std;
00065
00066 G4WentzelOKandVIxSection::G4WentzelOKandVIxSection() :
00067 numlimit(0.1),
00068 nwarnings(0),
00069 nwarnlimit(50),
00070 alpha2(fine_structure_const*fine_structure_const)
00071 {
00072 fNistManager = G4NistManager::Instance();
00073 fG4pow = G4Pow::GetInstance();
00074 theElectron = G4Electron::Electron();
00075 thePositron = G4Positron::Positron();
00076 theProton = G4Proton::Proton();
00077 lowEnergyLimit = 1.0*eV;
00078 G4double p0 = electron_mass_c2*classic_electr_radius;
00079 coeff = twopi*p0*p0;
00080 particle = 0;
00081
00082
00083
00084
00085 if(0.0 == ScreenRSquare[0]) {
00086 G4double a0 = electron_mass_c2/0.88534;
00087 G4double constn = 6.937e-6/(MeV*MeV);
00088
00089 ScreenRSquare[0] = alpha2*a0*a0;
00090 ScreenRSquareElec[0] = ScreenRSquare[0];
00091 for(G4int j=1; j<100; ++j) {
00092 G4double x = a0*fG4pow->Z13(j);
00093 if(1 == j) { ScreenRSquare[j] = 0.5*alpha2*a0*a0; }
00094 else {
00095 ScreenRSquare[j] = 0.5*(1 + exp(-j*j*0.001))*alpha2*x*x;
00096 ScreenRSquareElec[j] = 0.5*alpha2*x*x;
00097 }
00098 x = fNistManager->GetA27(j);
00099 FormFactor[j] = constn*x*x;
00100 }
00101 }
00102 currentMaterial = 0;
00103 elecXSRatio = factB = factD = formfactA = screenZ = 0.0;
00104 cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;
00105
00106 factB1= 0.5*CLHEP::pi*fine_structure_const;
00107
00108 tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
00109 ecut = etag = DBL_MAX;
00110 targetZ = 0;
00111 cosThetaMax = 1.0;
00112 targetMass = proton_mass_c2;
00113 particle = 0;
00114 }
00115
00116
00117
00118 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection()
00119 {}
00120
00121
00122
00123 void G4WentzelOKandVIxSection::Initialise(const G4ParticleDefinition* p,
00124 G4double CosThetaLim)
00125 {
00126 SetupParticle(p);
00127 tkin = mom2 = momCM2 = 0.0;
00128 ecut = etag = DBL_MAX;
00129 targetZ = 0;
00130 cosThetaMax = CosThetaLim;
00131 G4double a = G4LossTableManager::Instance()->FactorForAngleLimit()
00132 *CLHEP::hbarc/CLHEP::fermi;
00133 factorA2 = 0.5*a*a;
00134 currentMaterial = 0;
00135
00136
00137
00138
00139 }
00140
00141
00142
00143 void G4WentzelOKandVIxSection::SetupParticle(const G4ParticleDefinition* p)
00144 {
00145 particle = p;
00146 mass = particle->GetPDGMass();
00147 spin = particle->GetPDGSpin();
00148 if(0.0 != spin) { spin = 0.5; }
00149 G4double q = std::fabs(particle->GetPDGCharge()/eplus);
00150 chargeSquare = q*q;
00151 charge3 = chargeSquare*q;
00152 tkin = 0.0;
00153 currentMaterial = 0;
00154 targetZ = 0;
00155 }
00156
00157
00158
00159 G4double
00160 G4WentzelOKandVIxSection::SetupTarget(G4int Z, G4double cut)
00161 {
00162 G4double cosTetMaxNuc2 = cosTetMaxNuc;
00163 if(Z != targetZ || tkin != etag) {
00164 etag = tkin;
00165 targetZ = Z;
00166 if(targetZ > 99) { targetZ = 99; }
00167 SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2);
00168
00169
00170
00171
00172
00173
00174
00175 kinFactor = coeff*Z*chargeSquare*invbeta2/mom2;
00176
00177 if(1 == Z) {
00178 screenZ = ScreenRSquare[targetZ]/mom2;
00179 } else if(mass > MeV) {
00180 screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
00181 ScreenRSquare[targetZ]/mom2;
00182 } else {
00183 G4double tau = tkin/mass;
00184 screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z
00185 *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))*
00186 ScreenRSquareElec[targetZ]/mom2;
00187 }
00188 if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
00189 cosTetMaxNuc2 = 0.0;
00190 }
00191 formfactA = FormFactor[targetZ]*mom2;
00192
00193 cosTetMaxElec = 1.0;
00194 ComputeMaxElectronScattering(cut);
00195 }
00196 return cosTetMaxNuc2;
00197 }
00198
00199
00200
00201 G4double
00202 G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax)
00203 {
00204 G4double xsec = 0.0;
00205 if(cosTMax >= 1.0) { return xsec; }
00206
00207 G4double xSection = 0.0;
00208 G4double x = 0;
00209 G4double y = 0;
00210 G4double x1= 0;
00211 G4double x2= 0;
00212 G4double xlog = 0.0;
00213
00214 G4double costm = std::max(cosTMax,cosTetMaxElec);
00215 G4double fb = screenZ*factB;
00216
00217
00218 if(costm < 1.0) {
00219 x = (1.0 - costm)/screenZ;
00220 if(x < numlimit) {
00221 x2 = 0.5*x*x;
00222 y = x2*(1.0 - 1.3333333*x + 3*x2);
00223 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
00224 } else {
00225 x1= x/(1 + x);
00226 xlog = log(1.0 + x);
00227 y = xlog - x1;
00228 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
00229 }
00230
00231 if(y < 0.0) {
00232 ++nwarnings;
00233 if(nwarnings < nwarnlimit) {
00234 G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
00235 << G4endl;
00236 G4cout << "y= " << y
00237 << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
00238 << " Z= " << targetZ << " "
00239 << particle->GetParticleName() << G4endl;
00240 G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
00241 << " x= " << x << G4endl;
00242 }
00243 y = 0.0;
00244 }
00245 xSection = y;
00246 }
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256 if(cosTMax < 1.0) {
00257 x = (1.0 - cosTMax)/screenZ;
00258 if(x < numlimit) {
00259 x2 = 0.5*x*x;
00260 y = x2*(1.0 - 1.3333333*x + 3*x2);
00261 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
00262 } else {
00263 x1= x/(1 + x);
00264 xlog = log(1.0 + x);
00265 y = xlog - x1;
00266 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
00267 }
00268
00269 if(y < 0.0) {
00270 ++nwarnings;
00271 if(nwarnings < nwarnlimit) {
00272 G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
00273 << G4endl;
00274 G4cout << "y= " << y
00275 << " e(MeV)= " << tkin << " Z= " << targetZ << " "
00276 << particle->GetParticleName() << G4endl;
00277 G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
00278 << " x= " << " x1= " << x1 <<G4endl;
00279 }
00280 y = 0.0;
00281 }
00282 xSection += y*targetZ;
00283 }
00284 xSection *= kinFactor;
00285
00286
00287
00288
00289
00290
00291
00292
00293 return xSection;
00294 }
00295
00296
00297
00298 G4ThreeVector
00299 G4WentzelOKandVIxSection::SampleSingleScattering(G4double cosTMin,
00300 G4double cosTMax,
00301 G4double elecRatio)
00302 {
00303 G4ThreeVector v(0.0,0.0,1.0);
00304
00305 G4double formf = formfactA;
00306 G4double cost1 = cosTMin;
00307 G4double cost2 = cosTMax;
00308 if(elecRatio > 0.0) {
00309 if(G4UniformRand() <= elecRatio) {
00310 formf = 0.0;
00311 cost1 = std::max(cost1,cosTetMaxElec);
00312 cost2 = std::max(cost2,cosTetMaxElec);
00313 }
00314 }
00315 if(cost1 < cost2) { return v; }
00316
00317 G4double w1 = 1. - cost1 + screenZ;
00318 G4double w2 = 1. - cost2 + screenZ;
00319 G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;
00320
00321
00322 G4double fm = 1.0 + formf*z1;
00323
00324 G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm );
00325
00326 if( G4UniformRand() > grej ) { return v; }
00327
00328 G4double cost = 1.0 - z1;
00329
00330 if(cost > 1.0) { cost = 1.0; }
00331 else if(cost < -1.0) { cost =-1.0; }
00332 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
00333
00334 G4double phi = twopi*G4UniformRand();
00335 G4double vx1 = sint*cos(phi);
00336 G4double vy1 = sint*sin(phi);
00337
00338
00339 v.set(vx1,vy1,cost);
00340 return v;
00341 }
00342
00343
00344
00345 void
00346 G4WentzelOKandVIxSection::ComputeMaxElectronScattering(G4double cutEnergy)
00347 {
00348 if(mass > MeV) {
00349 G4double ratio = electron_mass_c2/mass;
00350 G4double tau = tkin/mass;
00351 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
00352 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
00353
00354 cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
00355 } else {
00356
00357 G4double tmax = tkin;
00358 if(particle == theElectron) { tmax *= 0.5; }
00359
00360 G4double t = std::min(cutEnergy, tmax);
00361 G4double mom21 = t*(t + 2.0*electron_mass_c2);
00362 G4double t1 = tkin - t;
00363
00364
00365 if(t1 > 0.0) {
00366 G4double mom22 = t1*(t1 + 2.0*mass);
00367 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
00368 if(ctm < 1.0) { cosTetMaxElec = ctm; }
00369 if(particle == theElectron && cosTetMaxElec < 0.0) { cosTetMaxElec = 0.0; }
00370 }
00371 }
00372 }
00373
00374