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 #include "G4eBremsstrahlungRelModel.hh"
00057 #include "G4PhysicalConstants.hh"
00058 #include "G4SystemOfUnits.hh"
00059 #include "G4Electron.hh"
00060 #include "G4Positron.hh"
00061 #include "G4Gamma.hh"
00062 #include "Randomize.hh"
00063 #include "G4Material.hh"
00064 #include "G4Element.hh"
00065 #include "G4ElementVector.hh"
00066 #include "G4ProductionCutsTable.hh"
00067 #include "G4ParticleChangeForLoss.hh"
00068 #include "G4LossTableManager.hh"
00069 #include "G4ModifiedTsai.hh"
00070 #include "G4DipBustGenerator.hh"
00071
00072
00073
00074 const G4double G4eBremsstrahlungRelModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
00075 0.5917, 0.7628, 0.8983, 0.9801 };
00076 const G4double G4eBremsstrahlungRelModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
00077 0.1813, 0.1569, 0.1112, 0.0506 };
00078 const G4double G4eBremsstrahlungRelModel::Fel_light[] = {0., 5.31 , 4.79 , 4.74 , 4.71} ;
00079 const G4double G4eBremsstrahlungRelModel::Finel_light[] = {0., 6.144 , 5.621 , 5.805 , 5.924} ;
00080
00081 using namespace std;
00082
00083 G4eBremsstrahlungRelModel::G4eBremsstrahlungRelModel(const G4ParticleDefinition* p,
00084 const G4String& nam)
00085 : G4VEmModel(nam),
00086 particle(0),
00087 bremFactor(fine_structure_const*classic_electr_radius*classic_electr_radius*16./3.),
00088 isElectron(true),
00089 fMigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi),
00090 fLPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)*0.5),
00091 fXiLPM(0), fPhiLPM(0), fGLPM(0),
00092 use_completescreening(false),isInitialised(false)
00093 {
00094 fParticleChange = 0;
00095 theGamma = G4Gamma::Gamma();
00096
00097 lowKinEnergy = 0.1*GeV;
00098 SetLowEnergyLimit(lowKinEnergy);
00099
00100 nist = G4NistManager::Instance();
00101
00102 SetLPMFlag(true);
00103
00104 SetAngularDistribution(new G4DipBustGenerator());
00105
00106 particleMass = kinEnergy = totalEnergy = currentZ = z13 = z23 = lnZ = Fel
00107 = Finel = fCoulomb = fMax = densityFactor = densityCorr = lpmEnergy
00108 = xiLPM = phiLPM = gLPM = klpm = kp = 0.0;
00109
00110 energyThresholdLPM = 1.e39;
00111
00112 InitialiseConstants();
00113 if(p) { SetParticle(p); }
00114 }
00115
00116
00117
00118 void G4eBremsstrahlungRelModel::InitialiseConstants()
00119 {
00120 facFel = log(184.15);
00121 facFinel = log(1194.);
00122
00123 preS1 = 1./(184.15*184.15);
00124 logTwo = log(2.);
00125 }
00126
00127
00128
00129 G4eBremsstrahlungRelModel::~G4eBremsstrahlungRelModel()
00130 {
00131 }
00132
00133
00134
00135 void G4eBremsstrahlungRelModel::SetParticle(const G4ParticleDefinition* p)
00136 {
00137 particle = p;
00138 particleMass = p->GetPDGMass();
00139 if(p == G4Electron::Electron()) { isElectron = true; }
00140 else { isElectron = false;}
00141 }
00142
00143
00144
00145 void G4eBremsstrahlungRelModel::SetupForMaterial(const G4ParticleDefinition*,
00146 const G4Material* mat,
00147 G4double kineticEnergy)
00148 {
00149 densityFactor = mat->GetElectronDensity()*fMigdalConstant;
00150 lpmEnergy = mat->GetRadlen()*fLPMconstant;
00151
00152
00153 if (LPMFlag()) {
00154 energyThresholdLPM=sqrt(densityFactor)*lpmEnergy;
00155 } else {
00156 energyThresholdLPM=1.e39;
00157 }
00158
00159 kinEnergy = kineticEnergy;
00160 totalEnergy = kineticEnergy + particleMass;
00161 densityCorr = densityFactor*totalEnergy*totalEnergy;
00162
00163
00164 klpm=totalEnergy*totalEnergy/lpmEnergy;
00165 kp=sqrt(densityCorr);
00166
00167 }
00168
00169
00170
00171
00172 void G4eBremsstrahlungRelModel::Initialise(const G4ParticleDefinition* p,
00173 const G4DataVector& cuts)
00174 {
00175 if(p) { SetParticle(p); }
00176
00177 lowKinEnergy = LowEnergyLimit();
00178
00179 currentZ = 0.;
00180
00181 InitialiseElementSelectors(p, cuts);
00182
00183 if(isInitialised) { return; }
00184 fParticleChange = GetParticleChangeForLoss();
00185 isInitialised = true;
00186 }
00187
00188
00189
00190 G4double G4eBremsstrahlungRelModel::ComputeDEDXPerVolume(
00191 const G4Material* material,
00192 const G4ParticleDefinition* p,
00193 G4double kineticEnergy,
00194 G4double cutEnergy)
00195 {
00196 if(!particle) { SetParticle(p); }
00197 if(kineticEnergy < lowKinEnergy) { return 0.0; }
00198 G4double cut = std::min(cutEnergy, kineticEnergy);
00199 if(cut == 0.0) { return 0.0; }
00200
00201 SetupForMaterial(particle, material,kineticEnergy);
00202
00203 const G4ElementVector* theElementVector = material->GetElementVector();
00204 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
00205
00206 G4double dedx = 0.0;
00207
00208
00209 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
00210
00211 G4VEmModel::SetCurrentElement((*theElementVector)[i]);
00212 SetCurrentElement((*theElementVector)[i]->GetZ());
00213
00214 dedx += theAtomicNumDensityVector[i]*currentZ*currentZ*ComputeBremLoss(cut);
00215 }
00216 dedx *= bremFactor;
00217
00218
00219 return dedx;
00220 }
00221
00222
00223
00224 G4double G4eBremsstrahlungRelModel::ComputeBremLoss(G4double cut)
00225 {
00226 G4double loss = 0.0;
00227
00228
00229 G4double vcut = cut/totalEnergy;
00230 G4int n = (G4int)(20*vcut) + 3;
00231 G4double delta = vcut/G4double(n);
00232
00233 G4double e0 = 0.0;
00234 G4double xs;
00235
00236
00237 for(G4int l=0; l<n; l++) {
00238
00239 for(G4int i=0; i<8; i++) {
00240
00241 G4double eg = (e0 + xgi[i]*delta)*totalEnergy;
00242
00243 if(totalEnergy > energyThresholdLPM) {
00244 xs = ComputeRelDXSectionPerAtom(eg);
00245 } else {
00246 xs = ComputeDXSectionPerAtom(eg);
00247 }
00248 loss += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
00249 }
00250 e0 += delta;
00251 }
00252
00253 loss *= delta*totalEnergy;
00254
00255 return loss;
00256 }
00257
00258
00259
00260 G4double G4eBremsstrahlungRelModel::ComputeCrossSectionPerAtom(
00261 const G4ParticleDefinition* p,
00262 G4double kineticEnergy,
00263 G4double Z, G4double,
00264 G4double cutEnergy,
00265 G4double maxEnergy)
00266 {
00267 if(!particle) { SetParticle(p); }
00268 if(kineticEnergy < lowKinEnergy) { return 0.0; }
00269 G4double cut = std::min(cutEnergy, kineticEnergy);
00270 G4double tmax = std::min(maxEnergy, kineticEnergy);
00271
00272 if(cut >= tmax) { return 0.0; }
00273
00274 SetCurrentElement(Z);
00275
00276 G4double cross = ComputeXSectionPerAtom(cut);
00277
00278
00279 if(tmax < kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
00280
00281 cross *= Z*Z*bremFactor;
00282
00283 return cross;
00284 }
00285
00286
00287
00288
00289 G4double G4eBremsstrahlungRelModel::ComputeXSectionPerAtom(G4double cut)
00290 {
00291 G4double cross = 0.0;
00292
00293
00294 G4double vcut = log(cut/totalEnergy);
00295 G4double vmax = log(kinEnergy/totalEnergy);
00296 G4int n = (G4int)(0.45*(vmax - vcut)) + 4;
00297
00298 G4double delta = (vmax - vcut)/G4double(n);
00299
00300 G4double e0 = vcut;
00301 G4double xs;
00302
00303
00304 for(G4int l=0; l<n; l++) {
00305
00306 for(G4int i=0; i<8; i++) {
00307
00308 G4double eg = exp(e0 + xgi[i]*delta)*totalEnergy;
00309
00310 if(totalEnergy > energyThresholdLPM) {
00311 xs = ComputeRelDXSectionPerAtom(eg);
00312 } else {
00313 xs = ComputeDXSectionPerAtom(eg);
00314 }
00315 cross += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
00316 }
00317 e0 += delta;
00318 }
00319
00320 cross *= delta;
00321
00322 return cross;
00323 }
00324
00325
00326 void G4eBremsstrahlungRelModel::CalcLPMFunctions(G4double k)
00327 {
00328
00329
00330 G4double sprime = sqrt(0.125*k*lpmEnergy/(totalEnergy*(totalEnergy-k)));
00331
00332 G4double s1 = preS1*z23;
00333 G4double logS1 = 2./3.*lnZ-2.*facFel;
00334 G4double logTS1 = logTwo+logS1;
00335
00336 xiLPM = 2.;
00337
00338 if (sprime>1)
00339 xiLPM = 1.;
00340 else if (sprime>sqrt(2.)*s1) {
00341 G4double h = log(sprime)/logTS1;
00342 xiLPM = 1+h-0.08*(1-h)*(1-sqr(1-h))/logTS1;
00343 }
00344
00345 G4double s0 = sprime/sqrt(xiLPM);
00346
00347
00348
00349 G4double k2 = k*k;
00350 s0 *= (1 + (densityCorr/k2) );
00351
00352
00353
00354 xiLPM = 1.;
00355 if (s0<=s1) xiLPM = 2.;
00356 else if ( (s1<s0) && (s0<=1) ) xiLPM = 1. + log(s0)/logS1;
00357
00358
00359
00360
00361 G4double s2=s0*s0;
00362 G4double s3=s0*s2;
00363 G4double s4=s2*s2;
00364
00365 if (s0<0.1) {
00366
00367 phiLPM = 6.*s0 - 18.84955592153876*s2 + 39.47841760435743*s3
00368 - 57.69873135166053*s4;
00369 gLPM = 37.69911184307752*s2 - 236.8705056261446*s3 + 807.7822389*s4;
00370 }
00371 else if (s0<1.9516) {
00372
00373
00374 phiLPM = 1.-exp(-6.*s0*(1.+(3.-pi)*s0)
00375 +s3/(0.623+0.795*s0+0.658*s2));
00376 if (s0<0.415827397755) {
00377
00378 G4double psiLPM = 1-exp(-4*s0-8*s2/(1+3.936*s0+4.97*s2-0.05*s3+7.50*s4));
00379 gLPM = 3*psiLPM-2*phiLPM;
00380 }
00381 else {
00382
00383 G4double pre = -0.16072300849123999 + s0*3.7550300067531581 + s2*-1.7981383069010097
00384 + s3*0.67282686077812381 + s4*-0.1207722909879257;
00385 gLPM = tanh(pre);
00386 }
00387 }
00388 else {
00389
00390 phiLPM = 1. - 0.0119048/s4;
00391 gLPM = 1. - 0.0230655/s4;
00392 }
00393
00394
00395
00396 if (xiLPM*phiLPM>1. || s0>0.57) { xiLPM=1./phiLPM; }
00397 }
00398
00399
00400
00401
00402 G4double G4eBremsstrahlungRelModel::ComputeRelDXSectionPerAtom(G4double gammaEnergy)
00403
00404
00405
00406 {
00407 if(gammaEnergy < 0.0) { return 0.0; }
00408
00409 G4double y = gammaEnergy/totalEnergy;
00410 G4double y2 = y*y*.25;
00411 G4double yone2 = (1.-y+2.*y2);
00412
00413
00414
00415
00416
00417 CalcLPMFunctions(gammaEnergy);
00418
00419 G4double mainLPM = xiLPM*(y2 * gLPM + yone2*phiLPM) * ( (Fel-fCoulomb) + Finel/currentZ );
00420 G4double secondTerm = (1.-y)/12.*(1.+1./currentZ);
00421
00422 G4double cross = mainLPM+secondTerm;
00423 return cross;
00424 }
00425
00426
00427
00428 G4double G4eBremsstrahlungRelModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
00429
00430
00431
00432
00433 {
00434
00435 if(gammaEnergy < 0.0) { return 0.0; }
00436
00437 G4double y = gammaEnergy/totalEnergy;
00438
00439 G4double main=0.,secondTerm=0.;
00440
00441 if (use_completescreening|| currentZ<5) {
00442
00443 main = (3./4.*y*y - y + 1.) * ( (Fel-fCoulomb) + Finel/currentZ );
00444 secondTerm = (1.-y)/12.*(1.+1./currentZ);
00445 }
00446 else {
00447
00448 G4double dd=100.*electron_mass_c2*y/(totalEnergy-gammaEnergy);
00449 G4double gg=dd/z13;
00450 G4double eps=dd/z23;
00451 G4double phi1=Phi1(gg,currentZ), phi1m2=Phi1M2(gg,currentZ);
00452 G4double psi1=Psi1(eps,currentZ), psi1m2=Psi1M2(eps,currentZ);
00453
00454 main = (3./4.*y*y - y + 1.) * ( (0.25*phi1-1./3.*lnZ-fCoulomb) + (0.25*psi1-2./3.*lnZ)/currentZ );
00455 secondTerm = (1.-y)/8.*(phi1m2+psi1m2/currentZ);
00456 }
00457 G4double cross = main+secondTerm;
00458 return cross;
00459 }
00460
00461
00462
00463 void G4eBremsstrahlungRelModel::SampleSecondaries(
00464 std::vector<G4DynamicParticle*>* vdp,
00465 const G4MaterialCutsCouple* couple,
00466 const G4DynamicParticle* dp,
00467 G4double cutEnergy,
00468 G4double maxEnergy)
00469 {
00470 G4double kineticEnergy = dp->GetKineticEnergy();
00471 if(kineticEnergy < lowKinEnergy) { return; }
00472 G4double cut = std::min(cutEnergy, kineticEnergy);
00473 G4double emax = std::min(maxEnergy, kineticEnergy);
00474 if(cut >= emax) { return; }
00475
00476 SetupForMaterial(particle, couple->GetMaterial(), kineticEnergy);
00477
00478 const G4Element* elm =
00479 SelectRandomAtom(couple,particle,kineticEnergy,cut,emax);
00480 SetCurrentElement(elm->GetZ());
00481
00482 kinEnergy = kineticEnergy;
00483 totalEnergy = kineticEnergy + particleMass;
00484 densityCorr = densityFactor*totalEnergy*totalEnergy;
00485
00486
00487 G4bool highe = true;
00488 if(totalEnergy < energyThresholdLPM) { highe = false; }
00489
00490 G4double xmin = log(cut*cut + densityCorr);
00491 G4double xmax = log(emax*emax + densityCorr);
00492 G4double gammaEnergy, f, x;
00493
00494 do {
00495 x = exp(xmin + G4UniformRand()*(xmax - xmin)) - densityCorr;
00496 if(x < 0.0) { x = 0.0; }
00497 gammaEnergy = sqrt(x);
00498 if(highe) { f = ComputeRelDXSectionPerAtom(gammaEnergy); }
00499 else { f = ComputeDXSectionPerAtom(gammaEnergy); }
00500
00501 if ( f > fMax ) {
00502 G4cout << "### G4eBremsstrahlungRelModel Warning: Majoranta exceeded! "
00503 << f << " > " << fMax
00504 << " Egamma(MeV)= " << gammaEnergy
00505 << " Ee(MeV)= " << kineticEnergy
00506 << " " << GetName()
00507 << G4endl;
00508 }
00509
00510 } while (f < fMax*G4UniformRand());
00511
00512
00513
00514
00515
00516
00517 G4ThreeVector gammaDirection =
00518 GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
00519 G4lrint(currentZ),
00520 couple->GetMaterial());
00521
00522
00523 G4DynamicParticle* gamma = new G4DynamicParticle(theGamma,gammaDirection,
00524 gammaEnergy);
00525 vdp->push_back(gamma);
00526
00527 G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
00528 G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
00529 - gammaEnergy*gammaDirection).unit();
00530
00531
00532 G4double finalE = kineticEnergy - gammaEnergy;
00533
00534
00535 if(gammaEnergy > SecondaryThreshold()) {
00536 fParticleChange->ProposeTrackStatus(fStopAndKill);
00537 fParticleChange->SetProposedKineticEnergy(0.0);
00538 G4DynamicParticle* el =
00539 new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
00540 direction, finalE);
00541 vdp->push_back(el);
00542
00543
00544 } else {
00545 fParticleChange->SetProposedMomentumDirection(direction);
00546 fParticleChange->SetProposedKineticEnergy(finalE);
00547 }
00548 }
00549
00550
00551
00552