G4eBremsstrahlungModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4eBremsstrahlungModel
00034 //
00035 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
00036 //
00037 // Creation date: 03.01.2002
00038 //
00039 // Modifications:
00040 //
00041 // 11-11-02  Fix division by 0 (V.Ivanchenko)
00042 // 04-12-02  Change G4DynamicParticle constructor in PostStep (V.Ivanchenko)
00043 // 23-12-02  Change interface in order to move to cut per region (V.Ivanchenko)
00044 // 24-01-03  Fix for compounds (V.Ivanchenko)
00045 // 27-01-03  Make models region aware (V.Ivanchenko)
00046 // 13-02-03  Add name (V.Ivanchenko)
00047 // 09-05-03  Fix problem of supression function + optimise sampling (V.Ivanchenko)
00048 // 20-05-04  Correction to ensure unit independence (L.Urban)
00049 // 08-04-05  Major optimisation of internal interfaces (V.Ivantchenko)
00050 // 03-08-05  Add extra protection at initialisation (V.Ivantchenko)
00051 // 07-02-06  public function ComputeCrossSectionPerAtom() (mma)
00052 // 21-03-06  Fix problem of initialisation in case when cuts are not defined (VI)
00053 // 27-03-06  Fix calculation of fl parameter at low energy (energy loss) (VI)
00054 // 15-02-07  correct LPMconstant by a factor 2, thanks to G. Depaola (mma)
00055 // 09-09-08  MigdalConstant increased in (2pi)^2 times (A.Schaelicke) 
00056 // 13-10-10  Add angular distributon interface (VI)
00057 //
00058 // Class Description:
00059 //
00060 //
00061 // -------------------------------------------------------------------
00062 //
00063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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   //  loop for elements in the material
00191   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
00192 
00193     G4double Z     = (*theElementVector)[i]->GetZ();
00194     G4double natom = theAtomicNumDensityVector[i];
00195 
00196     // loss for MinKinEnergy<KineticEnergy<=100 GeV
00197     if (kineticEnergy <= thigh) {
00198 
00199       //      x = log(totalEnergy/electron_mass_c2);
00200       loss = ComputeBremLoss(Z, kineticEnergy, cut) ;
00201       if (!isElectron) loss *= PositronCorrFactorLoss(Z, kineticEnergy, cut);
00202 
00203     // extrapolation for KineticEnergy>100 GeV
00204     } else {
00205 
00206       //      G4double xhigh = log(thigh/electron_mass_c2);
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     // now compute the correction due to the supression(s)
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       // correct the loss
00268       loss *= floss;
00269     }
00270     dedx += loss;
00271   }
00272   if(dedx < 0.) { dedx = 0.; }
00273   return dedx;
00274 }
00275 
00276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00277 
00278 G4double G4eBremsstrahlungModel::ComputeBremLoss(G4double Z, G4double T,
00279                                                  G4double Cut)
00280 
00281 // compute loss due to soft brems
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   // Z=2
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   // Z=4
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   // Z=6
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   // Z=14
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   // Z=26
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   // Z=50
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   // Z=82
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   // Z=92
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   //  correction
00371   loss *= (aaa+bbb*T/Tlim)/(1.+ccc*T/Tlim);
00372   loss *= fl;
00373   loss /= Avogadro;
00374 
00375   return loss;
00376 }
00377 
00378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00379 
00380 G4double G4eBremsstrahlungModel::PositronCorrFactorLoss(G4double Z,
00381                                  G4double kineticEnergy, G4double cut)
00382 
00383 //calculates the correction factor for the energy loss due to bremsstrahlung for positrons
00384 //the same correction is in the (discrete) bremsstrahlung
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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   // now compute the correction due to the supression(s)
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   // correct the cross section
00477   cross *= fsig;
00478 
00479   return cross;
00480 }
00481 
00482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00483 
00484 G4double G4eBremsstrahlungModel::ComputeCrossSectionPerAtom(
00485                                               const G4ParticleDefinition*,
00486                                                     G4double kineticEnergy, 
00487                                                     G4double Z,   G4double,
00488                                                     G4double cut, G4double)
00489  
00490 // Calculates the cross section per atom in GEANT4 internal units.
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   // Z=2
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   // Z=4
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   // Z=6
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   // Z=14
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   // Z=26
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   // Z=50
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   // Z=82
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   // Z=92
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00590  
00591 G4double G4eBremsstrahlungModel::PositronCorrFactorSigma( G4double Z,
00592                                  G4double kineticEnergy, G4double cut)
00593 
00594 //Calculates the correction factor for the total cross section of the positron
00595 //  bremsstrahl.
00596 // Eta is the ratio of positron to electron energy loss by bremstrahlung. 
00597 // A parametrized formula from L. Urban is used to estimate eta. It is a fit to
00598 // the results of L. Kim & al: Phys Rev. A33,3002 (1986)
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00613 
00614 G4DataVector* G4eBremsstrahlungModel::ComputePartialSumSigma(
00615               const G4Material* material,
00616                     G4double kineticEnergy,
00617                     G4double cut)
00618 
00619 // Build the table of cross section per element. 
00620 //The table is built for MATERIALS.
00621 // This table is used by DoIt to select randomly an element in the material.
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00642 
00643 void G4eBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp, 
00644                                                const G4MaterialCutsCouple* couple,
00645                                                const G4DynamicParticle* dp,
00646                                                G4double tmin,
00647                                                G4double maxEnergy)
00648 // The emitted gamma energy is sampled using a parametrized formula 
00649 // from L. Urban.
00650 // This parametrization is derived from :
00651 //    cross-section values of Seltzer and Berger for electron energies
00652 //    1 keV - 10 GeV,
00653 //    screened Bethe Heilter differential cross section above 10 GeV,
00654 //    Migdal corrections in both case.
00655 //  Seltzer & Berger: Nim B 12:95 (1985)
00656 //  Nelson, Hirayama & Rogers: Technical report 265 SLAC (1985)
00657 //  Migdal: Phys Rev 103:1811 (1956); Messel & Crawford: Pergamon Press (1970)
00658 //
00659 // A modified version of the random number techniques of Butcher&Messel is used
00660 //    (Nuc Phys 20(1960),15).
00661 {
00662   G4double kineticEnergy = dp->GetKineticEnergy();
00663   G4double tmax = min(maxEnergy, kineticEnergy);
00664   if(tmin >= tmax) { return; }
00665 
00666 //
00667 // GEANT4 internal units.
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   // select randomly one element constituing the material
00696   const G4Element* anElement = SelectRandomAtom(couple);
00697 
00698   // Extract Z factors for this Element
00699   G4double lnZ = 3.*(anElement->GetIonisation()->GetlogZ3());
00700   G4double FZ = lnZ* (4.- 0.55*lnZ);
00701   G4double ZZ = anElement->GetIonisation()->GetZZ3();
00702 
00703   // limits of the energy sampling
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   // Migdal factor
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   // precalculated parameters
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     // limit of the screening variable
00740     screenfac =
00741       136.*electron_mass_c2/((anElement->GetIonisation()->GetZ3())*totalEnergy);
00742     G4double screenmin = screenfac*epsilmin/(1.-epsilmin);
00743 
00744     // Compute the maximum of the rejection function
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     // Compute the maximum of the rejection function
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   //  sample the energy rate of the emitted gamma for e- kin energy > 1 MeV
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         if ( greject > grejmax ) {
00785             G4cout << "### G4eBremsstrahlungModel Warning: Majoranta exceeded! "
00786                    << greject << " > " << grejmax
00787                    << " x= " << x 
00788                    << " e= " << kineticEnergy
00789                    << G4endl;
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         if ( greject > grejmax ) {
00803           G4cout << "### G4eBremsstrahlungModel Warning: Majoranta exceeded! " 
00804                  << greject << " > " << grejmax 
00805                  << " x= " << x 
00806                  << " e= " << kineticEnergy
00807                  << G4endl;
00808         }
00809         */
00810       } while( greject < G4UniformRand()*grejmax );
00811     }
00812     gammaEnergy = x*kineticEnergy; 
00813 
00814     if (LPMFlag()) {
00815      // take into account the supression due to the LPM effect
00816       if (G4UniformRand() <= SupressionFunction(material,kineticEnergy,
00817                                                            gammaEnergy))
00818         LPMOK = true;
00819       }
00820     else LPMOK = true;
00821 
00822   } while (!LPMOK);
00823 
00824   //
00825   // angles of the emitted gamma. ( Z - axis along the parent particle)
00826   // use general interface
00827   //
00828 
00829   G4ThreeVector gammaDirection = 
00830     GetAngularDistribution()->SampleDirection(dp, totalEnergy-gammaEnergy,
00831                                               G4lrint(anElement->GetZ()), 
00832                                               couple->GetMaterial());
00833 
00834   // create G4DynamicParticle object for the Gamma
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   // energy of primary
00845   G4double finalE = kineticEnergy - gammaEnergy;
00846 
00847   // stop tracking and create new secondary instead of primary
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     // continue tracking
00857   } else {
00858     fParticleChange->SetProposedMomentumDirection(direction);
00859     fParticleChange->SetProposedKineticEnergy(finalE);
00860   }
00861 }
00862 
00863 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00864 
00865 const G4Element* G4eBremsstrahlungModel::SelectRandomAtom(
00866            const G4MaterialCutsCouple* couple) 
00867 {
00868   // select randomly 1 element within the material
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00896 
00897 G4double G4eBremsstrahlungModel::SupressionFunction(const G4Material* material,
00898                                  G4double kineticEnergy, G4double gammaEnergy)
00899 {
00900   // supression due to the LPM effect+polarisation of the medium/
00901   // supression due to the polarisation alone
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00942 
00943 

Generated on Mon May 27 17:48:05 2013 for Geant4 by  doxygen 1.4.7