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 "G4eBremParametrizedModel.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4Electron.hh"
00053 #include "G4Positron.hh"
00054 #include "G4Gamma.hh"
00055 #include "Randomize.hh"
00056 #include "G4Material.hh"
00057 #include "G4Element.hh"
00058 #include "G4ElementVector.hh"
00059 #include "G4ProductionCutsTable.hh"
00060 #include "G4ParticleChangeForLoss.hh"
00061 #include "G4LossTableManager.hh"
00062 #include "G4ModifiedTsai.hh"
00063
00064
00065
00066 const G4double G4eBremParametrizedModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
00067 0.5917, 0.7628, 0.8983, 0.9801 };
00068 const G4double G4eBremParametrizedModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
00069 0.1813, 0.1569, 0.1112, 0.0506 };
00070
00071
00072 using namespace std;
00073
00074 G4eBremParametrizedModel::G4eBremParametrizedModel(const G4ParticleDefinition* p,
00075 const G4String& nam)
00076 : G4VEmModel(nam),
00077 particle(0),
00078 isElectron(true),
00079 fMigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi),
00080 bremFactor(fine_structure_const*classic_electr_radius*classic_electr_radius*16./3.),
00081 isInitialised(false)
00082 {
00083 theGamma = G4Gamma::Gamma();
00084
00085 minThreshold = 0.1*keV;
00086 lowKinEnergy = 10.*MeV;
00087 SetLowEnergyLimit(lowKinEnergy);
00088
00089 nist = G4NistManager::Instance();
00090
00091 SetAngularDistribution(new G4ModifiedTsai());
00092
00093 particleMass = kinEnergy = totalEnergy = currentZ = z13 = z23 = lnZ = Fel = Finel
00094 = densityFactor = densityCorr = fMax = fCoulomb = 0.;
00095
00096 InitialiseConstants();
00097 if(p) { SetParticle(p); }
00098 }
00099
00100
00101
00102 void G4eBremParametrizedModel::InitialiseConstants()
00103 {
00104 facFel = log(184.15);
00105 facFinel = log(1194.);
00106 }
00107
00108
00109
00110 G4eBremParametrizedModel::~G4eBremParametrizedModel()
00111 {
00112 }
00113
00114
00115
00116 void G4eBremParametrizedModel::SetParticle(const G4ParticleDefinition* p)
00117 {
00118 particle = p;
00119 particleMass = p->GetPDGMass();
00120 if(p == G4Electron::Electron()) { isElectron = true; }
00121 else { isElectron = false;}
00122 }
00123
00124
00125
00126 G4double G4eBremParametrizedModel::MinEnergyCut(const G4ParticleDefinition*,
00127 const G4MaterialCutsCouple*)
00128 {
00129 return minThreshold;
00130 }
00131
00132
00133
00134 void G4eBremParametrizedModel::SetupForMaterial(const G4ParticleDefinition*,
00135 const G4Material* mat,
00136 G4double kineticEnergy)
00137 {
00138 densityFactor = mat->GetElectronDensity()*fMigdalConstant;
00139
00140
00141 kinEnergy = kineticEnergy;
00142 totalEnergy = kineticEnergy + particleMass;
00143 densityCorr = densityFactor*totalEnergy*totalEnergy;
00144 }
00145
00146
00147
00148
00149 void G4eBremParametrizedModel::Initialise(const G4ParticleDefinition* p,
00150 const G4DataVector& cuts)
00151 {
00152 if(p) { SetParticle(p); }
00153
00154 lowKinEnergy = LowEnergyLimit();
00155
00156 currentZ = 0.;
00157
00158 InitialiseElementSelectors(p, cuts);
00159
00160 if(isInitialised) { return; }
00161 fParticleChange = GetParticleChangeForLoss();
00162 isInitialised = true;
00163 }
00164
00165
00166
00167 G4double G4eBremParametrizedModel::ComputeDEDXPerVolume(
00168 const G4Material* material,
00169 const G4ParticleDefinition* p,
00170 G4double kineticEnergy,
00171 G4double cutEnergy)
00172 {
00173 if(!particle) { SetParticle(p); }
00174 if(kineticEnergy < lowKinEnergy) { return 0.0; }
00175 G4double cut = std::min(cutEnergy, kineticEnergy);
00176 if(cut == 0.0) { return 0.0; }
00177
00178 SetupForMaterial(particle, material,kineticEnergy);
00179
00180 const G4ElementVector* theElementVector = material->GetElementVector();
00181 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
00182
00183 G4double dedx = 0.0;
00184
00185
00186 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
00187
00188 G4VEmModel::SetCurrentElement((*theElementVector)[i]);
00189 SetCurrentElement((*theElementVector)[i]->GetZ());
00190
00191 dedx += theAtomicNumDensityVector[i]*currentZ*currentZ*ComputeBremLoss(cut);
00192 }
00193 dedx *= bremFactor;
00194
00195
00196 return dedx;
00197 }
00198
00199
00200
00201 G4double G4eBremParametrizedModel::ComputeBremLoss(G4double cut)
00202 {
00203 G4double loss = 0.0;
00204
00205
00206 G4double vcut = cut/totalEnergy;
00207 G4int n = (G4int)(20*vcut) + 3;
00208 G4double delta = vcut/G4double(n);
00209
00210 G4double e0 = 0.0;
00211 G4double xs;
00212
00213
00214 for(G4int l=0; l<n; l++) {
00215
00216 for(G4int i=0; i<8; i++) {
00217
00218 G4double eg = (e0 + xgi[i]*delta)*totalEnergy;
00219
00220 xs = ComputeDXSectionPerAtom(eg);
00221
00222 loss += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
00223 }
00224 e0 += delta;
00225 }
00226
00227 loss *= delta*totalEnergy;
00228
00229 return loss;
00230 }
00231
00232
00233
00234 G4double G4eBremParametrizedModel::ComputeCrossSectionPerAtom(
00235 const G4ParticleDefinition* p,
00236 G4double kineticEnergy,
00237 G4double Z, G4double,
00238 G4double cutEnergy,
00239 G4double maxEnergy)
00240 {
00241 if(!particle) { SetParticle(p); }
00242 if(kineticEnergy < lowKinEnergy) { return 0.0; }
00243 G4double cut = std::min(cutEnergy, kineticEnergy);
00244 G4double tmax = std::min(maxEnergy, kineticEnergy);
00245
00246 if(cut >= tmax) { return 0.0; }
00247
00248 SetCurrentElement(Z);
00249
00250 G4double cross = ComputeXSectionPerAtom(cut);
00251
00252
00253 if(tmax < kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
00254
00255 cross *= Z*Z*bremFactor;
00256
00257 return cross;
00258 }
00259
00260
00261
00262
00263 G4double G4eBremParametrizedModel::ComputeXSectionPerAtom(G4double cut)
00264 {
00265 G4double cross = 0.0;
00266
00267
00268 G4double vcut = log(cut/totalEnergy);
00269 G4double vmax = log(kinEnergy/totalEnergy);
00270 G4int n = (G4int)(0.45*(vmax - vcut)) + 4;
00271
00272 G4double delta = (vmax - vcut)/G4double(n);
00273
00274 G4double e0 = vcut;
00275 G4double xs;
00276
00277
00278 for(G4int l=0; l<n; l++) {
00279
00280 for(G4int i=0; i<8; i++) {
00281
00282 G4double eg = exp(e0 + xgi[i]*delta)*totalEnergy;
00283
00284 xs = ComputeDXSectionPerAtom(eg);
00285
00286 cross += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
00287 }
00288 e0 += delta;
00289 }
00290
00291 cross *= delta;
00292
00293 return cross;
00294 }
00295
00296
00297
00298
00299
00300
00301 static const G4double
00302 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
00303 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
00304 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
00305
00306 static const G4double
00307 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
00308 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
00309 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
00310
00311 static const G4double
00312 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
00313 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
00314 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
00315
00316 static const G4double
00317 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
00318 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
00319 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
00320
00321 static const G4double tlow = 1.*MeV;
00322
00323 G4double ScreenFunction1(G4double ScreenVariable)
00324
00325
00326
00327 {
00328 G4double screenVal;
00329
00330 if (ScreenVariable > 1.)
00331 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
00332 else
00333 screenVal = 42.392 - ScreenVariable* (7.796 - 1.961*ScreenVariable);
00334
00335 return screenVal;
00336 }
00337
00338
00339
00340 G4double ScreenFunction2(G4double ScreenVariable)
00341
00342
00343
00344 {
00345 G4double screenVal;
00346
00347 if (ScreenVariable > 1.)
00348 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
00349 else
00350 screenVal = 41.734 - ScreenVariable* (6.484 - 1.250*ScreenVariable);
00351
00352 return screenVal;
00353 }
00354
00355
00356
00357 G4double ComputeParametrizedDXSectionPerAtom(G4double kineticEnergy, G4double gammaEnergy, G4double Z) {
00358 G4double lnZ = std::log(Z);
00359 G4double FZ = lnZ* (4.- 0.55*lnZ);
00360 G4double ZZ = std::pow (Z*(Z+1.),1./3.);
00361 G4double Z3 = std::pow (Z,1./3.);
00362
00363 G4double totalEnergy = kineticEnergy + electron_mass_c2;
00364
00365
00366 G4double epsil, greject;
00367 G4double U = log(kineticEnergy/electron_mass_c2);
00368 G4double U2 = U*U;
00369
00370
00371 G4double ah, bh;
00372
00373 if (kineticEnergy > tlow) {
00374
00375 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12);
00376 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22);
00377 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
00378
00379 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12);
00380 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22);
00381 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
00382
00383 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
00384 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
00385
00386
00387 G4double screenfac =
00388 136.*electron_mass_c2/(Z3*totalEnergy);
00389
00390 epsil = gammaEnergy/totalEnergy;
00391 G4double screenvar = screenfac*epsil/(1.0-epsil);
00392 G4double F1 = max(ScreenFunction1(screenvar) - FZ ,0.);
00393 G4double F2 = max(ScreenFunction2(screenvar) - FZ ,0.);
00394
00395
00396 greject = (F1 - epsil* (ah*F1 - bh*epsil*F2))/8.;
00397
00398 std::cout << " yy = "<<epsil<<std::endl;
00399 std::cout << " F1/(...) "<<F1/(42.392 - FZ)<<std::endl;
00400 std::cout << " F2/(...) "<<F2/(42.392 - FZ)<<std::endl;
00401 std::cout << " (42.392 - FZ) " << (42.392 - FZ) <<std::endl;
00402
00403 } else {
00404
00405 G4double al0 = al00 + ZZ* (al01 + ZZ* al02);
00406 G4double al1 = al10 + ZZ* (al11 + ZZ* al12);
00407 G4double al2 = al20 + ZZ* (al21 + ZZ* al22);
00408
00409 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02);
00410 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12);
00411 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
00412
00413 ah = al0 + al1*U + al2*U2;
00414 bh = bl0 + bl1*U + bl2*U2;
00415
00416 G4double x=gammaEnergy/kineticEnergy;
00417 greject=(1. + x* (ah + bh*x));
00418
00419
00420
00421
00422
00423
00424
00425 }
00426
00427 return greject;
00428 }
00429
00430
00431
00432 G4double G4eBremParametrizedModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
00433 {
00434
00435 if(gammaEnergy < 0.0) { return 0.0; }
00436
00437 G4double y = gammaEnergy/totalEnergy;
00438
00439 G4double main=0.;
00440
00441
00442
00443
00444 main = (3./4.*y*y - y + 1.) * ( (Fel-fCoulomb) + Finel/currentZ );
00445
00446
00447 std::cout<<" F1(0) "<<ScreenFunction1(0.) <<std::endl;
00448 std::cout<<" F1(0) "<<ScreenFunction2(0.) <<std::endl;
00449 std::cout<<"Ekin = "<<kinEnergy<<std::endl;
00450 std::cout<<"Z = "<<currentZ<<std::endl;
00451 std::cout<<"main = "<<main<<std::endl;
00452 std::cout<<" y = "<<y<<std::endl;
00453 std::cout<<" Fel-fCoulomb "<< (Fel-fCoulomb) <<std::endl;
00454
00455 G4double main2 = ComputeParametrizedDXSectionPerAtom(kinEnergy,gammaEnergy,currentZ);
00456 std::cout<<"main2 = "<<main2<<std::endl;
00457 std::cout<<"main2tot = "<<main2 * ( (Fel-fCoulomb) + Finel/currentZ ) / (Fel-fCoulomb);
00458
00459
00460 G4double cross = main2;
00461 return cross;
00462 }
00463
00464
00465
00466 void G4eBremParametrizedModel::SampleSecondaries(
00467 std::vector<G4DynamicParticle*>* vdp,
00468 const G4MaterialCutsCouple* couple,
00469 const G4DynamicParticle* dp,
00470 G4double cutEnergy,
00471 G4double maxEnergy)
00472 {
00473 G4double kineticEnergy = dp->GetKineticEnergy();
00474 if(kineticEnergy < lowKinEnergy) { return; }
00475 G4double cut = std::min(cutEnergy, kineticEnergy);
00476 G4double emax = std::min(maxEnergy, kineticEnergy);
00477 if(cut >= emax) { return; }
00478
00479 SetupForMaterial(particle, couple->GetMaterial(),kineticEnergy);
00480
00481 const G4Element* elm =
00482 SelectRandomAtom(couple,particle,kineticEnergy,cut,emax);
00483 SetCurrentElement(elm->GetZ());
00484
00485 kinEnergy = kineticEnergy;
00486 totalEnergy = kineticEnergy + particleMass;
00487 densityCorr = densityFactor*totalEnergy*totalEnergy;
00488
00489 G4double xmin = log(cut*cut + densityCorr);
00490 G4double xmax = log(emax*emax + densityCorr);
00491 G4double gammaEnergy, f, x;
00492
00493 do {
00494 x = exp(xmin + G4UniformRand()*(xmax - xmin)) - densityCorr;
00495 if(x < 0.0) x = 0.0;
00496 gammaEnergy = sqrt(x);
00497 f = ComputeDXSectionPerAtom(gammaEnergy);
00498
00499 if ( f > fMax ) {
00500 G4cout << "### G4eBremParametrizedModel Warning: Majoranta exceeded! "
00501 << f << " > " << fMax
00502 << " Egamma(MeV)= " << gammaEnergy
00503 << " E(mEV)= " << kineticEnergy
00504 << G4endl;
00505 }
00506
00507 } while (f < fMax*G4UniformRand());
00508
00509
00510
00511
00512
00513 G4ThreeVector gammaDirection =
00514 GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
00515 G4lrint(currentZ),
00516 couple->GetMaterial());
00517
00518
00519 G4DynamicParticle* gamma = new G4DynamicParticle(theGamma,gammaDirection,
00520 gammaEnergy);
00521 vdp->push_back(gamma);
00522
00523 G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
00524 G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
00525 - gammaEnergy*gammaDirection).unit();
00526
00527
00528 G4double finalE = kineticEnergy - gammaEnergy;
00529
00530
00531 if(gammaEnergy > SecondaryThreshold()) {
00532 fParticleChange->ProposeTrackStatus(fStopAndKill);
00533 fParticleChange->SetProposedKineticEnergy(0.0);
00534 G4DynamicParticle* el =
00535 new G4DynamicParticle(const_cast<G4ParticleDefinition*>(particle),
00536 direction, finalE);
00537 vdp->push_back(el);
00538
00539
00540 } else {
00541 fParticleChange->SetProposedMomentumDirection(direction);
00542 fParticleChange->SetProposedKineticEnergy(finalE);
00543 }
00544 }
00545
00546
00547
00548