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 "G4eBremsstrahlungModel.hh"
00067 #include "G4PhysicalConstants.hh"
00068 #include "G4SystemOfUnits.hh"
00069 #include "G4Electron.hh"
00070 #include "G4Positron.hh"
00071 #include "G4Gamma.hh"
00072 #include "Randomize.hh"
00073 #include "G4Material.hh"
00074 #include "G4Element.hh"
00075 #include "G4ElementVector.hh"
00076 #include "G4ProductionCutsTable.hh"
00077 #include "G4DataVector.hh"
00078 #include "G4ParticleChangeForLoss.hh"
00079 #include "G4ModifiedTsai.hh"
00080
00081
00082
00083 using namespace std;
00084
00085 G4eBremsstrahlungModel::G4eBremsstrahlungModel(const G4ParticleDefinition* p,
00086 const G4String& nam)
00087 : G4VEmModel(nam),
00088 particle(0),
00089 isElectron(true),
00090 probsup(1.0),
00091 MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi),
00092 LPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)),
00093 isInitialised(false)
00094 {
00095 if(p) { SetParticle(p); }
00096 theGamma = G4Gamma::Gamma();
00097 SetAngularDistribution(new G4ModifiedTsai());
00098 highKinEnergy = HighEnergyLimit();
00099 lowKinEnergy = LowEnergyLimit();
00100 fParticleChange = 0;
00101 }
00102
00103
00104
00105 G4eBremsstrahlungModel::~G4eBremsstrahlungModel()
00106 {
00107 size_t n = partialSumSigma.size();
00108 if(n > 0) {
00109 for(size_t i=0; i<n; i++) {
00110 delete partialSumSigma[i];
00111 }
00112 }
00113 }
00114
00115
00116
00117 void G4eBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p)
00118 {
00119 particle = p;
00120 if(p == G4Electron::Electron()) { isElectron = true; }
00121 else { isElectron = false;}
00122 }
00123
00124
00125
00126 void G4eBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
00127 const G4DataVector& cuts)
00128 {
00129 if(p) { SetParticle(p); }
00130 highKinEnergy = HighEnergyLimit();
00131 lowKinEnergy = LowEnergyLimit();
00132 const G4ProductionCutsTable* theCoupleTable=
00133 G4ProductionCutsTable::GetProductionCutsTable();
00134
00135 if(theCoupleTable) {
00136 G4int numOfCouples = theCoupleTable->GetTableSize();
00137
00138 G4int nn = partialSumSigma.size();
00139 G4int nc = cuts.size();
00140 if(nn > 0) {
00141 for (G4int ii=0; ii<nn; ii++){
00142 G4DataVector* a=partialSumSigma[ii];
00143 if ( a ) delete a;
00144 }
00145 partialSumSigma.clear();
00146 }
00147 if(numOfCouples>0) {
00148 for (G4int i=0; i<numOfCouples; i++) {
00149 G4double cute = DBL_MAX;
00150 if(i < nc) cute = cuts[i];
00151 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
00152 const G4Material* material = couple->GetMaterial();
00153 G4DataVector* dv = ComputePartialSumSigma(material, 0.5*highKinEnergy,
00154 std::min(cute, 0.25*highKinEnergy));
00155 partialSumSigma.push_back(dv);
00156 }
00157 }
00158 }
00159 if(isInitialised) return;
00160 fParticleChange = GetParticleChangeForLoss();
00161 isInitialised = true;
00162 }
00163
00164
00165
00166 G4double G4eBremsstrahlungModel::ComputeDEDXPerVolume(
00167 const G4Material* material,
00168 const G4ParticleDefinition* p,
00169 G4double kineticEnergy,
00170 G4double cutEnergy)
00171 {
00172 if(!particle) { SetParticle(p); }
00173 if(kineticEnergy < lowKinEnergy) { return 0.0; }
00174
00175 const G4double thigh = 100.*GeV;
00176
00177 G4double cut = std::min(cutEnergy, kineticEnergy);
00178
00179 G4double rate, loss;
00180 const G4double factorHigh = 36./(1450.*GeV);
00181 const G4double coef1 = -0.5;
00182 const G4double coef2 = 2./9.;
00183
00184 const G4ElementVector* theElementVector = material->GetElementVector();
00185 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
00186
00187 G4double totalEnergy = kineticEnergy + electron_mass_c2;
00188 G4double dedx = 0.0;
00189
00190
00191 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
00192
00193 G4double Z = (*theElementVector)[i]->GetZ();
00194 G4double natom = theAtomicNumDensityVector[i];
00195
00196
00197 if (kineticEnergy <= thigh) {
00198
00199
00200 loss = ComputeBremLoss(Z, kineticEnergy, cut) ;
00201 if (!isElectron) loss *= PositronCorrFactorLoss(Z, kineticEnergy, cut);
00202
00203
00204 } else {
00205
00206
00207 G4double cuthigh = thigh*0.5;
00208
00209 if (cut < thigh) {
00210
00211 loss = ComputeBremLoss(Z, thigh, cut) ;
00212 if (!isElectron) loss *= PositronCorrFactorLoss(Z, thigh, cut) ;
00213 rate = cut/totalEnergy;
00214 loss *= (1. + coef1*rate + coef2*rate*rate);
00215 rate = cut/thigh;
00216 loss /= (1.+coef1*rate+coef2*rate*rate);
00217
00218 } else {
00219
00220 loss = ComputeBremLoss(Z, thigh, cuthigh) ;
00221 if (!isElectron) loss *= PositronCorrFactorLoss(Z, thigh, cuthigh) ;
00222 rate = cut/totalEnergy;
00223 loss *= (1. + coef1*rate + coef2*rate*rate);
00224 loss *= cut*factorHigh;
00225 }
00226 }
00227 loss *= natom;
00228
00229 G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
00230 * (material->GetElectronDensity()) ;
00231
00232
00233 G4double kmin = 1.*eV;
00234 G4double kmax = cut;
00235
00236 if (kmax > kmin) {
00237
00238 G4double floss = 0.;
00239 G4int nmax = 100;
00240
00241 G4double vmin=log(kmin);
00242 G4double vmax=log(kmax) ;
00243 G4int nn = (G4int)(nmax*(vmax-vmin)/(log(highKinEnergy)-vmin)) ;
00244 G4double u,fac,c,v,dv ;
00245 if(nn > 0) {
00246
00247 dv = (vmax-vmin)/nn ;
00248 v = vmin-dv ;
00249
00250 for(G4int n=0; n<=nn; n++) {
00251
00252 v += dv;
00253 u = exp(v);
00254 fac = u*SupressionFunction(material,kineticEnergy,u);
00255 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
00256 if ((n==0)||(n==nn)) c=0.5;
00257 else c=1. ;
00258 fac *= c ;
00259 floss += fac ;
00260 }
00261 floss *=dv/(kmax-kmin);
00262
00263 } else {
00264 floss = 1.;
00265 }
00266 if(floss > 1.) floss = 1.;
00267
00268 loss *= floss;
00269 }
00270 dedx += loss;
00271 }
00272 if(dedx < 0.) { dedx = 0.; }
00273 return dedx;
00274 }
00275
00276
00277
00278 G4double G4eBremsstrahlungModel::ComputeBremLoss(G4double Z, G4double T,
00279 G4double Cut)
00280
00281
00282 {
00283 static const G4double beta=1.0, ksi=2.0;
00284 static const G4double clossh = 0.254 , closslow = 1./3. , alosslow = 1. ;
00285 static const G4double Tlim= 10.*MeV ;
00286
00287 static const G4double xlim = 1.2 ;
00288 static const G4int NZ = 8 ;
00289 static const G4int Nloss = 11 ;
00290 static const G4double ZZ[NZ] =
00291 {2.,4.,6.,14.,26.,50.,82.,92.};
00292 static const G4double coefloss[NZ][Nloss] = {
00293
00294 { 0.98916, 0.47564, -0.2505, -0.45186, 0.14462,
00295 0.21307, -0.013738, -0.045689, -0.0042914, 0.0034429,
00296 0.00064189},
00297
00298
00299 { 1.0626, 0.37662, -0.23646, -0.45188, 0.14295,
00300 0.22906, -0.011041, -0.051398, -0.0055123, 0.0039919,
00301 0.00078003},
00302
00303 { 1.0954, 0.315, -0.24011, -0.43849, 0.15017,
00304 0.23001, -0.012846, -0.052555, -0.0055114, 0.0041283,
00305 0.00080318},
00306
00307
00308 { 1.1649, 0.18976, -0.24972, -0.30124, 0.1555,
00309 0.13565, -0.024765, -0.027047, -0.00059821, 0.0019373,
00310 0.00027647},
00311
00312
00313 { 1.2261, 0.14272, -0.25672, -0.28407, 0.13874,
00314 0.13586, -0.020562, -0.026722, -0.00089557, 0.0018665,
00315 0.00026981},
00316
00317
00318 { 1.3147, 0.020049, -0.35543, -0.13927, 0.17666,
00319 0.073746, -0.036076, -0.013407, 0.0025727, 0.00084005,
00320 -1.4082e-05},
00321
00322
00323 { 1.3986, -0.10586, -0.49187, -0.0048846, 0.23621,
00324 0.031652, -0.052938, -0.0076639, 0.0048181, 0.00056486,
00325 -0.00011995},
00326
00327
00328 { 1.4217, -0.116, -0.55497, -0.044075, 0.27506,
00329 0.081364, -0.058143, -0.023402, 0.0031322, 0.0020201,
00330 0.00017519}
00331
00332 } ;
00333 static G4double aaa = 0.414;
00334 static G4double bbb = 0.345;
00335 static G4double ccc = 0.460;
00336
00337 G4int iz = 0;
00338 G4double delz = 1.e6;
00339 for (G4int ii=0; ii<NZ; ii++)
00340 {
00341 G4double dz = std::abs(Z-ZZ[ii]);
00342 if(dz < delz) {
00343 iz = ii;
00344 delz = dz;
00345 }
00346 }
00347
00348 G4double xx = log10(T/MeV);
00349 G4double fl = 1.;
00350
00351 if (xx <= xlim)
00352 {
00353 xx /= xlim;
00354 G4double yy = 1.0;
00355 fl = 0.0;
00356 for (G4int j=0; j<Nloss; j++) {
00357 fl += yy+coefloss[iz][j];
00358 yy *= xx;
00359 }
00360 if (fl < 0.00001) fl = 0.00001;
00361 else if (fl > 1.0) fl = 1.0;
00362 }
00363
00364 G4double loss;
00365 G4double E = T+electron_mass_c2 ;
00366
00367 loss = Z*(Z+ksi)*E*E/(T+E)*exp(beta*log(Cut/T))*(2.-clossh*exp(log(Z)/4.));
00368 if (T <= Tlim) loss /= exp(closslow*log(Tlim/T));
00369 if( T <= Cut) loss *= exp(alosslow*log(T/Cut));
00370
00371 loss *= (aaa+bbb*T/Tlim)/(1.+ccc*T/Tlim);
00372 loss *= fl;
00373 loss /= Avogadro;
00374
00375 return loss;
00376 }
00377
00378
00379
00380 G4double G4eBremsstrahlungModel::PositronCorrFactorLoss(G4double Z,
00381 G4double kineticEnergy, G4double cut)
00382
00383
00384
00385
00386 {
00387 static const G4double K = 132.9416*eV ;
00388 static const G4double a1=4.15e-1, a3=2.10e-3, a5=54.0e-5 ;
00389
00390 G4double x = log(kineticEnergy/(K*Z*Z)), x2 = x*x, x3 = x2*x;
00391 G4double eta = 0.5+atan(a1*x+a3*x3+a5*x3*x2)/pi;
00392 G4double e0 = cut/kineticEnergy;
00393
00394 G4double factor = 0.0;
00395 if (e0 < 1.0) {
00396 factor=log(1.-e0)/eta;
00397 factor=exp(factor);
00398 }
00399 factor = eta*(1.-factor)/e0;
00400
00401 return factor;
00402 }
00403
00404
00405
00406 G4double G4eBremsstrahlungModel::CrossSectionPerVolume(
00407 const G4Material* material,
00408 const G4ParticleDefinition* p,
00409 G4double kineticEnergy,
00410 G4double cutEnergy,
00411 G4double maxEnergy)
00412 {
00413 if(!particle) { SetParticle(p); }
00414 G4double cross = 0.0;
00415 G4double tmax = min(maxEnergy, kineticEnergy);
00416 G4double cut = min(cutEnergy, kineticEnergy);
00417 if(cut >= tmax) { return cross; }
00418
00419 const G4ElementVector* theElementVector = material->GetElementVector();
00420 const G4double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
00421 G4double dum=0.;
00422
00423 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
00424
00425 cross += theAtomNumDensityVector[i] * ComputeCrossSectionPerAtom(p,
00426 kineticEnergy, (*theElementVector)[i]->GetZ(), dum, cut);
00427 if (tmax < kineticEnergy) {
00428 cross -= theAtomNumDensityVector[i] * ComputeCrossSectionPerAtom(p,
00429 kineticEnergy, (*theElementVector)[i]->GetZ(), dum, tmax);
00430 }
00431 }
00432
00433
00434
00435 G4double kmax = tmax;
00436 G4double kmin = cut;
00437
00438 G4double totalEnergy = kineticEnergy+electron_mass_c2 ;
00439 G4double kp2 = MigdalConstant*totalEnergy*totalEnergy
00440 *(material->GetElectronDensity());
00441
00442 G4double fsig = 0.;
00443 G4int nmax = 100;
00444 G4double vmin=log(kmin);
00445 G4double vmax=log(kmax) ;
00446 G4int nn = (G4int)(nmax*(vmax-vmin)/(log(highKinEnergy)-vmin));
00447 G4double u,fac,c,v,dv,y ;
00448 if(nn > 0) {
00449
00450 dv = (vmax-vmin)/nn ;
00451 v = vmin-dv ;
00452 for(G4int n=0; n<=nn; n++) {
00453
00454 v += dv;
00455 u = exp(v);
00456 fac = SupressionFunction(material, kineticEnergy, u);
00457 y = u/kmax;
00458 fac *= (4.-4.*y+3.*y*y)/3.;
00459 fac *= probsup*(u*u/(u*u+kp2))+1.-probsup;
00460
00461 if ((n==0)||(n==nn)) c=0.5;
00462 else c=1. ;
00463
00464 fac *= c;
00465 fsig += fac;
00466 }
00467 y = kmin/kmax ;
00468 fsig *=dv/(-4.*log(y)/3.-4.*(1.-y)/3.+0.5*(1.-y*y));
00469
00470 } else {
00471
00472 fsig = 1.;
00473 }
00474 if (fsig > 1.) fsig = 1.;
00475
00476
00477 cross *= fsig;
00478
00479 return cross;
00480 }
00481
00482
00483
00484 G4double G4eBremsstrahlungModel::ComputeCrossSectionPerAtom(
00485 const G4ParticleDefinition*,
00486 G4double kineticEnergy,
00487 G4double Z, G4double,
00488 G4double cut, G4double)
00489
00490
00491
00492
00493 {
00494 G4double cross = 0.0 ;
00495 if ( kineticEnergy < keV || kineticEnergy < cut) { return cross; }
00496
00497 static const G4double ksi=2.0, alfa=1.00;
00498 static const G4double csigh = 0.127, csiglow = 0.25, asiglow = 0.020*MeV ;
00499 static const G4double Tlim = 10.*MeV ;
00500
00501 static const G4double xlim = 1.2 ;
00502 static const G4int NZ = 8 ;
00503 static const G4int Nsig = 11 ;
00504 static const G4double ZZ[NZ] =
00505 {2.,4.,6.,14.,26.,50.,82.,92.} ;
00506 static const G4double coefsig[NZ][Nsig] = {
00507
00508 { 0.4638, 0.37748, 0.32249, -0.060362, -0.065004,
00509 -0.033457, -0.004583, 0.011954, 0.0030404, -0.0010077,
00510 -0.00028131},
00511
00512
00513 { 0.50008, 0.33483, 0.34364, -0.086262, -0.055361,
00514 -0.028168, -0.0056172, 0.011129, 0.0027528, -0.00092265,
00515 -0.00024348},
00516
00517
00518 { 0.51587, 0.31095, 0.34996, -0.11623, -0.056167,
00519 -0.0087154, 0.00053943, 0.0054092, 0.00077685, -0.00039635,
00520 -6.7818e-05},
00521
00522
00523 { 0.55058, 0.25629, 0.35854, -0.080656, -0.054308,
00524 -0.049933, -0.00064246, 0.016597, 0.0021789, -0.001327,
00525 -0.00025983},
00526
00527
00528 { 0.5791, 0.26152, 0.38953, -0.17104, -0.099172,
00529 0.024596, 0.023718, -0.0039205, -0.0036658, 0.00041749,
00530 0.00023408},
00531
00532
00533 { 0.62085, 0.27045, 0.39073, -0.37916, -0.18878,
00534 0.23905, 0.095028, -0.068744, -0.023809, 0.0062408,
00535 0.0020407},
00536
00537
00538 { 0.66053, 0.24513, 0.35404, -0.47275, -0.22837,
00539 0.35647, 0.13203, -0.1049, -0.034851, 0.0095046,
00540 0.0030535},
00541
00542
00543 { 0.67143, 0.23079, 0.32256, -0.46248, -0.20013,
00544 0.3506, 0.11779, -0.1024, -0.032013, 0.0092279,
00545 0.0028592}
00546
00547 } ;
00548
00549 G4int iz = 0 ;
00550 G4double delz = 1.e6 ;
00551 for (G4int ii=0; ii<NZ; ii++)
00552 {
00553 G4double absdelz = std::abs(Z-ZZ[ii]);
00554 if(absdelz < delz)
00555 {
00556 iz = ii ;
00557 delz = absdelz;
00558 }
00559 }
00560
00561 G4double xx = log10(kineticEnergy/MeV);
00562 G4double fs = 1. ;
00563
00564 if (xx <= xlim) {
00565
00566 fs = coefsig[iz][Nsig-1] ;
00567 for (G4int j=Nsig-2; j>=0; j--) {
00568
00569 fs = fs*xx+coefsig[iz][j] ;
00570 }
00571 if(fs < 0.) { fs = 0.; }
00572 }
00573
00574 cross = Z*(Z+ksi)*(1.-csigh*exp(log(Z)/4.))*pow(log(kineticEnergy/cut),alfa);
00575
00576 if (kineticEnergy <= Tlim)
00577 cross *= exp(csiglow*log(Tlim/kineticEnergy))
00578 *(1.+asiglow/(sqrt(Z)*kineticEnergy));
00579
00580 if (!isElectron)
00581 cross *= PositronCorrFactorSigma(Z, kineticEnergy, cut);
00582
00583 cross *= fs/Avogadro ;
00584
00585 if (cross < 0.) { cross = 0.; }
00586 return cross;
00587 }
00588
00589
00590
00591 G4double G4eBremsstrahlungModel::PositronCorrFactorSigma( G4double Z,
00592 G4double kineticEnergy, G4double cut)
00593
00594
00595
00596
00597
00598
00599
00600 {
00601 static const G4double K = 132.9416*eV;
00602 static const G4double a1 = 4.15e-1, a3 = 2.10e-3, a5 = 54.0e-5;
00603
00604 G4double x = log(kineticEnergy/(K*Z*Z));
00605 G4double x2 = x*x;
00606 G4double x3 = x2*x;
00607 G4double eta = 0.5 + atan(a1*x + a3*x3 + a5*x3*x2)/pi ;
00608 G4double alfa = (1. - eta)/eta;
00609 return eta*pow((1. - cut/kineticEnergy), alfa);
00610 }
00611
00612
00613
00614 G4DataVector* G4eBremsstrahlungModel::ComputePartialSumSigma(
00615 const G4Material* material,
00616 G4double kineticEnergy,
00617 G4double cut)
00618
00619
00620
00621
00622 {
00623 G4int nElements = material->GetNumberOfElements();
00624 const G4ElementVector* theElementVector = material->GetElementVector();
00625 const G4double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
00626 G4double dum = 0.;
00627
00628 G4DataVector* dv = new G4DataVector();
00629
00630 G4double cross = 0.0;
00631
00632 for (G4int i=0; i<nElements; i++ ) {
00633
00634 cross += theAtomNumDensityVector[i] * ComputeCrossSectionPerAtom( particle,
00635 kineticEnergy, (*theElementVector)[i]->GetZ(), dum, cut);
00636 dv->push_back(cross);
00637 }
00638 return dv;
00639 }
00640
00641
00642
00643 void G4eBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
00644 const G4MaterialCutsCouple* couple,
00645 const G4DynamicParticle* dp,
00646 G4double tmin,
00647 G4double maxEnergy)
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661 {
00662 G4double kineticEnergy = dp->GetKineticEnergy();
00663 G4double tmax = min(maxEnergy, kineticEnergy);
00664 if(tmin >= tmax) { return; }
00665
00666
00667
00668
00669 static const G4double
00670 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
00671 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
00672 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
00673
00674 static const G4double
00675 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
00676 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
00677 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
00678
00679 static const G4double
00680 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
00681 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
00682 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
00683
00684 static const G4double
00685 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
00686 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
00687 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
00688
00689 static const G4double tlow = 1.*MeV;
00690
00691 G4double gammaEnergy;
00692 G4bool LPMOK = false;
00693 const G4Material* material = couple->GetMaterial();
00694
00695
00696 const G4Element* anElement = SelectRandomAtom(couple);
00697
00698
00699 G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3());
00700 G4double FZ = lnZ* (4.- 0.55*lnZ);
00701 G4double ZZ = anElement->GetIonisation()->GetZZ3();
00702
00703
00704 G4double totalEnergy = kineticEnergy + electron_mass_c2;
00705
00706 G4double xmin = tmin/kineticEnergy;
00707 G4double xmax = tmax/kineticEnergy;
00708 G4double kappa = 0.0;
00709 if(xmax >= 1.) { xmax = 1.; }
00710 else { kappa = log(xmax)/log(xmin); }
00711 G4double epsilmin = tmin/totalEnergy;
00712 G4double epsilmax = tmax/totalEnergy;
00713
00714
00715 G4double MigdalFactor = (material->GetElectronDensity())*MigdalConstant
00716 / (epsilmax*epsilmax);
00717
00718 G4double x, epsil, greject, migdal, grejmax, q;
00719 G4double U = log(kineticEnergy/electron_mass_c2);
00720 G4double U2 = U*U;
00721
00722
00723 G4double ah, bh;
00724 G4double screenfac = 0.0;
00725
00726 if (kineticEnergy > tlow) {
00727
00728 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12);
00729 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22);
00730 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
00731
00732 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12);
00733 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22);
00734 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
00735
00736 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
00737 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
00738
00739
00740 screenfac =
00741 136.*electron_mass_c2/((anElement->GetIonisation()->GetZ3())*totalEnergy);
00742 G4double screenmin = screenfac*epsilmin/(1.-epsilmin);
00743
00744
00745 G4double F1 = max(ScreenFunction1(screenmin) - FZ ,0.);
00746 G4double F2 = max(ScreenFunction2(screenmin) - FZ ,0.);
00747 grejmax = (F1 - epsilmin* (F1*ah - bh*epsilmin*F2))/(42.392 - FZ);
00748
00749 } else {
00750
00751 G4double al0 = al00 + ZZ* (al01 + ZZ* al02);
00752 G4double al1 = al10 + ZZ* (al11 + ZZ* al12);
00753 G4double al2 = al20 + ZZ* (al21 + ZZ* al22);
00754
00755 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02);
00756 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12);
00757 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
00758
00759 ah = al0 + al1*U + al2*U2;
00760 bh = bl0 + bl1*U + bl2*U2;
00761
00762
00763 grejmax = max(1. + xmin* (ah + bh*xmin), 1.+ah+bh);
00764 G4double xm = -ah/(2.*bh);
00765 if ( xmin < xm && xm < xmax) grejmax = max(grejmax, 1.+ xm* (ah + bh*xm));
00766 }
00767
00768
00769
00770
00771
00772 do {
00773 if (kineticEnergy > tlow) {
00774 do {
00775 q = G4UniformRand();
00776 x = pow(xmin, q + kappa*(1.0 - q));
00777 epsil = x*kineticEnergy/totalEnergy;
00778 G4double screenvar = screenfac*epsil/(1.0-epsil);
00779 G4double F1 = max(ScreenFunction1(screenvar) - FZ ,0.);
00780 G4double F2 = max(ScreenFunction2(screenvar) - FZ ,0.);
00781 migdal = (1. + MigdalFactor)/(1. + MigdalFactor/(x*x));
00782 greject = migdal*(F1 - epsil* (ah*F1 - bh*epsil*F2))/(42.392 - FZ);
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792 } while( greject < G4UniformRand()*grejmax );
00793
00794 } else {
00795
00796 do {
00797 q = G4UniformRand();
00798 x = pow(xmin, q + kappa*(1.0 - q));
00799 migdal = (1. + MigdalFactor)/(1. + MigdalFactor/(x*x));
00800 greject = migdal*(1. + x* (ah + bh*x));
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810 } while( greject < G4UniformRand()*grejmax );
00811 }
00812 gammaEnergy = x*kineticEnergy;
00813
00814 if (LPMFlag()) {
00815
00816 if (G4UniformRand() <= SupressionFunction(material,kineticEnergy,
00817 gammaEnergy))
00818 LPMOK = true;
00819 }
00820 else LPMOK = true;
00821
00822 } while (!LPMOK);
00823
00824
00825
00826
00827
00828
00829 G4ThreeVector gammaDirection =
00830 GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
00831 G4lrint(anElement->GetZ()),
00832 couple->GetMaterial());
00833
00834
00835 G4DynamicParticle* gamma = new G4DynamicParticle(theGamma,gammaDirection,
00836 gammaEnergy);
00837 vdp->push_back(gamma);
00838
00839 G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
00840
00841 G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
00842 - gammaEnergy*gammaDirection).unit();
00843
00844
00845 G4double finalE = kineticEnergy - gammaEnergy;
00846
00847
00848 if(gammaEnergy > SecondaryThreshold()) {
00849 fParticleChange->ProposeTrackStatus(fStopAndKill);
00850 fParticleChange->SetProposedKineticEnergy(0.0);
00851 G4DynamicParticle* el =
00852 new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
00853 direction, finalE);
00854 vdp->push_back(el);
00855
00856
00857 } else {
00858 fParticleChange->SetProposedMomentumDirection(direction);
00859 fParticleChange->SetProposedKineticEnergy(finalE);
00860 }
00861 }
00862
00863
00864
00865 const G4Element* G4eBremsstrahlungModel::SelectRandomAtom(
00866 const G4MaterialCutsCouple* couple)
00867 {
00868
00869
00870 const G4Material* material = couple->GetMaterial();
00871 G4int nElements = material->GetNumberOfElements();
00872 const G4ElementVector* theElementVector = material->GetElementVector();
00873
00874 const G4Element* elm = 0;
00875
00876 if(1 < nElements) {
00877
00878 --nElements;
00879 G4DataVector* dv = partialSumSigma[couple->GetIndex()];
00880 G4double rval = G4UniformRand()*((*dv)[nElements]);
00881
00882 elm = (*theElementVector)[nElements];
00883 for (G4int i=0; i<nElements; ++i) {
00884 if (rval <= (*dv)[i]) {
00885 elm = (*theElementVector)[i];
00886 break;
00887 }
00888 }
00889 } else { elm = (*theElementVector)[0]; }
00890
00891 SetCurrentElement(elm);
00892 return elm;
00893 }
00894
00895
00896
00897 G4double G4eBremsstrahlungModel::SupressionFunction(const G4Material* material,
00898 G4double kineticEnergy, G4double gammaEnergy)
00899 {
00900
00901
00902
00903
00904 G4double totEnergy = kineticEnergy+electron_mass_c2 ;
00905 G4double totEnergySquare = totEnergy*totEnergy ;
00906
00907 G4double LPMEnergy = LPMconstant*(material->GetRadlen()) ;
00908
00909 G4double gammaEnergySquare = gammaEnergy*gammaEnergy ;
00910
00911 G4double electronDensity = material->GetElectronDensity();
00912
00913 G4double sp = gammaEnergySquare/
00914 (gammaEnergySquare+MigdalConstant*totEnergySquare*electronDensity);
00915
00916 G4double supr = 1.0;
00917
00918 if (LPMFlag()) {
00919
00920 G4double s2lpm = LPMEnergy*gammaEnergy/totEnergySquare;
00921
00922 if (s2lpm < 1.) {
00923
00924 G4double LPMgEnergyLimit = totEnergySquare/LPMEnergy ;
00925 G4double LPMgEnergyLimit2 = LPMgEnergyLimit*LPMgEnergyLimit;
00926 G4double splim = LPMgEnergyLimit2/
00927 (LPMgEnergyLimit2+MigdalConstant*totEnergySquare*electronDensity);
00928 G4double w = 1.+1./splim ;
00929
00930 if ((1.-sp) < 1.e-6) w = s2lpm*(3.-sp);
00931 else w = s2lpm*(1.+1./sp);
00932
00933 supr = (sqrt(w*w+4.*s2lpm)-w)/(sqrt(w*w+4.)-w) ;
00934 supr /= sp;
00935 }
00936
00937 }
00938 return supr;
00939 }
00940
00941
00942
00943