G4eBremsstrahlungModel Class Reference

#include <G4eBremsstrahlungModel.hh>

Inheritance diagram for G4eBremsstrahlungModel:

G4VEmModel G4ePolarizedBremsstrahlungModel

Public Member Functions

 G4eBremsstrahlungModel (const G4ParticleDefinition *p=0, const G4String &nam="eBrem")
virtual ~G4eBremsstrahlungModel ()
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double tkin, G4double Z, G4double, G4double cut, G4double maxE=DBL_MAX)
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)

Protected Member Functions

const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *couple)

Protected Attributes

const G4ParticleDefinitionparticle
G4ParticleDefinitiontheGamma
G4ParticleChangeForLossfParticleChange
G4bool isElectron

Detailed Description

Definition at line 64 of file G4eBremsstrahlungModel.hh.


Constructor & Destructor Documentation

G4eBremsstrahlungModel::G4eBremsstrahlungModel ( const G4ParticleDefinition p = 0,
const G4String nam = "eBrem" 
)

Definition at line 85 of file G4eBremsstrahlungModel.cc.

References fParticleChange, G4Gamma::Gamma(), G4VEmModel::HighEnergyLimit(), G4VEmModel::LowEnergyLimit(), G4VEmModel::SetAngularDistribution(), and theGamma.

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 }

G4eBremsstrahlungModel::~G4eBremsstrahlungModel (  )  [virtual]

Definition at line 105 of file G4eBremsstrahlungModel.cc.

References CLHEP::detail::n.

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 }


Member Function Documentation

G4double G4eBremsstrahlungModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition ,
G4double  tkin,
G4double  Z,
G4double  ,
G4double  cut,
G4double  maxE = DBL_MAX 
) [virtual]

Reimplemented from G4VEmModel.

Definition at line 484 of file G4eBremsstrahlungModel.cc.

References isElectron.

Referenced by CrossSectionPerVolume().

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 }

G4double G4eBremsstrahlungModel::ComputeDEDXPerVolume ( const G4Material ,
const G4ParticleDefinition ,
G4double  kineticEnergy,
G4double  cutEnergy 
) [virtual]

Reimplemented from G4VEmModel.

Definition at line 166 of file G4eBremsstrahlungModel.cc.

References G4Material::GetAtomicNumDensityVector(), G4Material::GetElectronDensity(), G4Material::GetElementVector(), G4Material::GetNumberOfElements(), isElectron, CLHEP::detail::n, G4InuclParticleNames::nn, and particle.

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 }

G4double G4eBremsstrahlungModel::CrossSectionPerVolume ( const G4Material ,
const G4ParticleDefinition ,
G4double  kineticEnergy,
G4double  cutEnergy,
G4double  maxEnergy 
) [virtual]

Reimplemented from G4VEmModel.

Definition at line 406 of file G4eBremsstrahlungModel.cc.

References ComputeCrossSectionPerAtom(), G4Material::GetAtomicNumDensityVector(), G4Material::GetElectronDensity(), G4Material::GetElementVector(), G4Material::GetNumberOfElements(), CLHEP::detail::n, G4InuclParticleNames::nn, and particle.

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 }

void G4eBremsstrahlungModel::Initialise ( const G4ParticleDefinition ,
const G4DataVector  
) [virtual]

Implements G4VEmModel.

Reimplemented in G4ePolarizedBremsstrahlungModel.

Definition at line 126 of file G4eBremsstrahlungModel.cc.

References DBL_MAX, fParticleChange, G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4VEmModel::GetParticleChangeForLoss(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4VEmModel::HighEnergyLimit(), G4VEmModel::LowEnergyLimit(), and G4InuclParticleNames::nn.

Referenced by G4ePolarizedBremsstrahlungModel::Initialise().

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 }

void G4eBremsstrahlungModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  ,
const G4MaterialCutsCouple ,
const G4DynamicParticle ,
G4double  tmin,
G4double  maxEnergy 
) [virtual]

Implements G4VEmModel.

Reimplemented in G4ePolarizedBremsstrahlungModel.

Definition at line 643 of file G4eBremsstrahlungModel.cc.

References fParticleChange, fStopAndKill, G4lrint(), G4UniformRand, G4VEmModel::GetAngularDistribution(), G4Material::GetElectronDensity(), G4Element::GetIonisation(), G4DynamicParticle::GetKineticEnergy(), G4IonisParamElm::GetlogZ3(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4Element::GetZ(), G4IonisParamElm::GetZ3(), G4IonisParamElm::GetZZ3(), G4VEmModel::LPMFlag(), G4VParticleChange::ProposeTrackStatus(), G4VEmAngularDistribution::SampleDirection(), G4VEmModel::SecondaryThreshold(), SelectRandomAtom(), G4ParticleChangeForLoss::SetProposedKineticEnergy(), G4ParticleChangeForLoss::SetProposedMomentumDirection(), and theGamma.

Referenced by G4ePolarizedBremsstrahlungModel::SampleSecondaries().

00650                                         :
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 }

const G4Element * G4eBremsstrahlungModel::SelectRandomAtom ( const G4MaterialCutsCouple couple  )  [protected]

Definition at line 865 of file G4eBremsstrahlungModel.cc.

References G4UniformRand, G4Material::GetElementVector(), G4MaterialCutsCouple::GetIndex(), G4MaterialCutsCouple::GetMaterial(), G4Material::GetNumberOfElements(), and G4VEmModel::SetCurrentElement().

Referenced by SampleSecondaries().

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 }


Field Documentation

G4ParticleChangeForLoss* G4eBremsstrahlungModel::fParticleChange [protected]

Definition at line 131 of file G4eBremsstrahlungModel.hh.

Referenced by G4eBremsstrahlungModel(), Initialise(), G4ePolarizedBremsstrahlungModel::SampleSecondaries(), and SampleSecondaries().

G4bool G4eBremsstrahlungModel::isElectron [protected]

Definition at line 133 of file G4eBremsstrahlungModel.hh.

Referenced by ComputeCrossSectionPerAtom(), and ComputeDEDXPerVolume().

const G4ParticleDefinition* G4eBremsstrahlungModel::particle [protected]

Definition at line 129 of file G4eBremsstrahlungModel.hh.

Referenced by ComputeDEDXPerVolume(), and CrossSectionPerVolume().

G4ParticleDefinition* G4eBremsstrahlungModel::theGamma [protected]

Definition at line 130 of file G4eBremsstrahlungModel.hh.

Referenced by G4eBremsstrahlungModel(), and SampleSecondaries().


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:49 2013 for Geant4 by  doxygen 1.4.7