G4UrbanMscModel92 Class Reference

#include <G4UrbanMscModel92.hh>

Inheritance diagram for G4UrbanMscModel92:

G4VMscModel G4VEmModel

Public Member Functions

 G4UrbanMscModel92 (const G4String &nam="UrbanMsc92")
virtual ~G4UrbanMscModel92 ()
void Initialise (const G4ParticleDefinition *, const G4DataVector &)
void StartTracking (G4Track *)
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=0., G4double emax=DBL_MAX)
G4ThreeVector & SampleScattering (const G4ThreeVector &, G4double safety)
G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep)
G4double ComputeGeomPathLength (G4double truePathLength)
G4double ComputeTrueStepLength (G4double geomStepLength)
G4double ComputeTheta0 (G4double truePathLength, G4double KineticEnergy)

Detailed Description

Definition at line 70 of file G4UrbanMscModel92.hh.

Constructor & Destructor Documentation

G4UrbanMscModel92::G4UrbanMscModel92 ( const G4String &  nam = "UrbanMsc92"  ) 

Definition at line 119 of file G4UrbanMscModel92.cc.

References G4cout, G4endl, G4LossTableManager::Instance(), G4INCL::Math::pi, and G4VMscModel::skin.

00120   : G4VMscModel(nam)
00121 {
00122   masslimite    = 0.6*MeV;
00123   lambdalimit   = 1.*mm;
00124   fr            = 0.02;
00125   //  facsafety     = 0.3;
00126   taubig        = 8.0;
00127   tausmall      = 1.e-16;
00128   taulim        = 1.e-6;
00129   currentTau    = taulim;
00130   tlimitminfix  = 1.e-6*mm;            
00131   stepmin       = tlimitminfix;
00132   smallstep     = 1.e10;
00133   currentRange  = 0. ;
00134   rangeinit     = 0.;
00135   tlimit        = 1.e10*mm;
00136   tlimitmin     = 10.*tlimitminfix;            
00137   tgeom         = 1.e50*mm;
00138   geombig       = 1.e50*mm;
00139   geommin       = 1.e-3*mm;
00140   geomlimit     = geombig;
00141   presafety     = 0.*mm;
00142                           
00143   y             = 0.;
00144 
00145   Zold          = 0.;
00146   Zeff          = 1.;
00147   Z2            = 1.;                
00148   Z23           = 1.;                    
00149   lnZ           = 0.;
00150   coeffth1      = 0.;
00151   coeffth2      = 0.;
00152   coeffc1       = 0.;
00153   coeffc2       = 0.;
00154   scr1ini       = fine_structure_const*fine_structure_const*
00155                   electron_mass_c2*electron_mass_c2/(0.885*0.885*4.*pi);
00156   scr2ini       = 3.76*fine_structure_const*fine_structure_const;
00157   scr1          = 0.;
00158   scr2          = 0.;
00159 
00160   theta0max     = pi/6.;
00161   rellossmax    = 0.50;
00162   third         = 1./3.;
00163   particle      = 0;
00164   theManager    = G4LossTableManager::Instance(); 
00165   firstStep     = true; 
00166   inside        = false;  
00167   insideskin    = false;
00168 
00169   skindepth = skin*stepmin;
00170 
00171   mass = proton_mass_c2;
00172   charge = ChargeSquare = 1.0;
00173   currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength 
00174     = zPathLength = par1 = par2 = par3 = 0;
00175 
00176   currentMaterialIndex = -1;
00177   fParticleChange = 0;
00178   couple = 0;
00179   G4cout << "### G4UrbanMscModel92 is obsolete and will be removed for "
00180          << "the next Geant4 version" << G4endl;
00181 }

G4UrbanMscModel92::~G4UrbanMscModel92 (  )  [virtual]

Definition at line 185 of file G4UrbanMscModel92.cc.

00186 {}

Member Function Documentation

G4double G4UrbanMscModel92::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  particle,
G4double  KineticEnergy,
G4double  AtomicNumber,
G4double  AtomicWeight = 0.,
G4double  cut = 0.,
G4double  emax = DBL_MAX 
) [virtual]

Reimplemented from G4VEmModel.

Definition at line 209 of file G4UrbanMscModel92.cc.

00214 {
00215   const G4double sigmafactor = twopi*classic_electr_radius*classic_electr_radius;
00216   const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2*
00217                             Bohr_radius*Bohr_radius/(hbarc*hbarc);
00218   const G4double epsmin = 1.e-4 , epsmax = 1.e10;
00219 
00220   const G4double Zdat[15] = { 4.,  6., 13., 20., 26., 29., 32., 38., 47.,
00221                              50., 56., 64., 74., 79., 82. };
00222 
00223   const G4double Tdat[22] = { 100*eV,  200*eV,  400*eV,  700*eV,
00224                                1*keV,   2*keV,   4*keV,   7*keV,
00225                               10*keV,  20*keV,  40*keV,  70*keV,
00226                              100*keV, 200*keV, 400*keV, 700*keV,
00227                                1*MeV,   2*MeV,   4*MeV,   7*MeV,
00228                               10*MeV,  20*MeV};
00229 
00230   // corr. factors for e-/e+ lambda for T <= Tlim
00231           G4double celectron[15][22] =
00232           {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
00233             1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
00234             1.112,1.108,1.100,1.093,1.089,1.087            },
00235            {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
00236             1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
00237             1.109,1.105,1.097,1.090,1.086,1.082            },
00238            {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
00239             1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
00240             1.131,1.124,1.113,1.104,1.099,1.098            },
00241            {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
00242             1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
00243             1.112,1.105,1.096,1.089,1.085,1.098            },
00244            {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
00245             1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
00246             1.073,1.070,1.064,1.059,1.056,1.056            },
00247            {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
00248             1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
00249             1.074,1.070,1.063,1.059,1.056,1.052            },
00250            {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
00251             1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
00252             1.068,1.064,1.059,1.054,1.051,1.050            },
00253            {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
00254             1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
00255             1.039,1.037,1.034,1.031,1.030,1.036            },
00256            {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
00257             1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
00258             1.031,1.028,1.024,1.022,1.021,1.024            },
00259            {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
00260             1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
00261             1.020,1.017,1.015,1.013,1.013,1.020            },
00262            {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
00263             1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
00264             0.995,0.993,0.993,0.993,0.993,1.011            },
00265            {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
00266             1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
00267             0.974,0.972,0.973,0.974,0.975,0.987            },
00268            {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
00269             1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
00270             0.950,0.947,0.949,0.952,0.954,0.963            },
00271            {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
00272             1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
00273             0.941,0.938,0.940,0.944,0.946,0.954            },
00274            {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
00275             1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
00276             0.933,0.930,0.933,0.936,0.939,0.949            }};
00277             
00278            G4double cpositron[15][22] = {
00279            {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
00280             1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
00281             1.131,1.126,1.117,1.108,1.103,1.100            },
00282            {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
00283             1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
00284             1.138,1.132,1.122,1.113,1.108,1.102            },
00285            {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
00286             1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
00287             1.203,1.190,1.173,1.159,1.151,1.145            },
00288            {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
00289             1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
00290             1.225,1.210,1.191,1.175,1.166,1.174            },
00291            {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
00292             1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
00293             1.217,1.203,1.184,1.169,1.160,1.151            },
00294            {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
00295             1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
00296             1.237,1.222,1.201,1.184,1.174,1.159            },
00297            {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
00298             1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
00299             1.252,1.234,1.212,1.194,1.183,1.170            },
00300            {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
00301             2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
00302             1.254,1.237,1.214,1.195,1.185,1.179            },
00303            {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
00304             2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
00305             1.312,1.288,1.258,1.235,1.221,1.205            },
00306            {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
00307             2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
00308             1.320,1.294,1.264,1.240,1.226,1.214            },
00309            {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
00310             2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
00311             1.328,1.302,1.270,1.245,1.231,1.233            },
00312            {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
00313             2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
00314             1.361,1.330,1.294,1.267,1.251,1.239            },
00315            {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
00316             3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
00317             1.409,1.372,1.330,1.298,1.280,1.258            },
00318            {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
00319             3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
00320             1.442,1.400,1.354,1.319,1.299,1.272            },
00321            {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
00322             3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
00323             1.456,1.412,1.364,1.328,1.307,1.282            }};
00324 
00325   //data/corrections for T > Tlim  
00326   G4double Tlim = 10.*MeV;
00327   G4double beta2lim = Tlim*(Tlim+2.*electron_mass_c2)/
00328                       ((Tlim+electron_mass_c2)*(Tlim+electron_mass_c2));
00329   G4double bg2lim   = Tlim*(Tlim+2.*electron_mass_c2)/
00330                       (electron_mass_c2*electron_mass_c2);
00331 
00332   G4double sig0[15] = {0.2672*barn,  0.5922*barn, 2.653*barn,  6.235*barn,
00333                       11.69*barn  , 13.24*barn  , 16.12*barn, 23.00*barn ,
00334                       35.13*barn  , 39.95*barn  , 50.85*barn, 67.19*barn ,
00335                       91.15*barn  , 104.4*barn  , 113.1*barn};
00336                                        
00337   G4double hecorr[15] = {120.70, 117.50, 105.00, 92.92, 79.23,  74.510,  68.29,
00338                           57.39,  41.97,  36.14, 24.53, 10.21,  -7.855, -16.84,
00339                          -22.30};
00340 
00341   G4double sigma;
00342   SetParticle(part);
00343 
00344   Z23 = pow(AtomicNumber,2./3.); 
00345 
00346   // correction if particle .ne. e-/e+
00347   // compute equivalent kinetic energy
00348   // lambda depends on p*beta ....
00349 
00350   G4double eKineticEnergy = KineticEnergy;
00351 
00352   if(mass > electron_mass_c2)
00353   {
00354      G4double TAU = KineticEnergy/mass ;
00355      G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
00356      G4double w = c-2. ;
00357      G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
00358      eKineticEnergy = electron_mass_c2*tau ;
00359   }
00360 
00361   G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
00362   G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
00363                                  /(eTotalEnergy*eTotalEnergy);
00364   G4double bg2   = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
00365                                  /(electron_mass_c2*electron_mass_c2);
00366 
00367   G4double eps = epsfactor*bg2/Z23;
00368 
00369   if     (eps<epsmin)  sigma = 2.*eps*eps;
00370   else if(eps<epsmax)  sigma = log(1.+2.*eps)-2.*eps/(1.+2.*eps);
00371   else                 sigma = log(2.*eps)-1.+1./eps;
00372 
00373   sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
00374 
00375   // interpolate in AtomicNumber and beta2 
00376   G4double c1,c2,cc1,cc2,corr;
00377 
00378   // get bin number in Z
00379   G4int iZ = 14;
00380   while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
00381   if (iZ==14)                               iZ = 13;
00382   if (iZ==-1)                               iZ = 0 ;
00383 
00384   G4double ZZ1 = Zdat[iZ];
00385   G4double ZZ2 = Zdat[iZ+1];
00386   G4double ratZ = (AtomicNumber-ZZ1)*(AtomicNumber+ZZ1)/
00387                   ((ZZ2-ZZ1)*(ZZ2+ZZ1));
00388 
00389   if(eKineticEnergy <= Tlim) 
00390   {
00391     // get bin number in T (beta2)
00392     G4int iT = 21;
00393     while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
00394     if(iT==21)                                  iT = 20;
00395     if(iT==-1)                                  iT = 0 ;
00396 
00397     //  calculate betasquare values
00398     G4double T = Tdat[iT],   E = T + electron_mass_c2;
00399     G4double b2small = T*(E+electron_mass_c2)/(E*E);
00400 
00401     T = Tdat[iT+1]; E = T + electron_mass_c2;
00402     G4double b2big = T*(E+electron_mass_c2)/(E*E);
00403     G4double ratb2 = (beta2-b2small)/(b2big-b2small);
00404 
00405     if (charge < 0.)
00406     {
00407        c1 = celectron[iZ][iT];
00408        c2 = celectron[iZ+1][iT];
00409        cc1 = c1+ratZ*(c2-c1);
00410 
00411        c1 = celectron[iZ][iT+1];
00412        c2 = celectron[iZ+1][iT+1];
00413        cc2 = c1+ratZ*(c2-c1);
00414 
00415        corr = cc1+ratb2*(cc2-cc1);
00416 
00417        sigma *= sigmafactor/corr;
00418     }
00419     else              
00420     {
00421        c1 = cpositron[iZ][iT];
00422        c2 = cpositron[iZ+1][iT];
00423        cc1 = c1+ratZ*(c2-c1);
00424 
00425        c1 = cpositron[iZ][iT+1];
00426        c2 = cpositron[iZ+1][iT+1];
00427        cc2 = c1+ratZ*(c2-c1);
00428 
00429        corr = cc1+ratb2*(cc2-cc1);
00430 
00431        sigma *= sigmafactor/corr;
00432     }
00433   }
00434   else
00435   {
00436     c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
00437     c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
00438     if((AtomicNumber >= ZZ1) && (AtomicNumber <= ZZ2))
00439       sigma = c1+ratZ*(c2-c1) ;
00440     else if(AtomicNumber < ZZ1)
00441       sigma = AtomicNumber*AtomicNumber*c1/(ZZ1*ZZ1);
00442     else if(AtomicNumber > ZZ2)
00443       sigma = AtomicNumber*AtomicNumber*c2/(ZZ2*ZZ2);
00444   }
00445   return sigma;
00446 
00447 }

G4double G4UrbanMscModel92::ComputeGeomPathLength ( G4double  truePathLength  )  [virtual]

Reimplemented from G4VMscModel.

Definition at line 654 of file G4UrbanMscModel92.cc.

References G4VMscModel::dtrl, G4UniformRand, G4VMscModel::GetEnergy(), G4VMscModel::GetTransportMeanFreePath(), and G4VMscModel::samplez.

00655 {
00656   firstStep = false; 
00657   lambdaeff = lambda0;
00658   par1 = -1. ;  
00659   par2 = par3 = 0. ;  
00660 
00661   //  do the true -> geom transformation
00662   zPathLength = tPathLength;
00663 
00664   // z = t for very small tPathLength
00665   if(tPathLength < tlimitminfix) return zPathLength;
00666 
00667   // this correction needed to run MSC with eIoni and eBrem inactivated
00668   // and makes no harm for a normal run
00669   if(tPathLength > currentRange)
00670     tPathLength = currentRange ;
00671 
00672   G4double tau   = tPathLength/lambda0 ;
00673 
00674   if ((tau <= tausmall) || insideskin) {
00675     zPathLength  = tPathLength;
00676     if(zPathLength > lambda0) zPathLength = lambda0;
00677     return zPathLength;
00678   }
00679 
00680   G4double zmean = tPathLength;
00681   if (tPathLength < currentRange*dtrl) {
00682     if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
00683     else             zmean = lambda0*(1.-exp(-tau));
00684   } else if(currentKinEnergy < mass || tPathLength == currentRange)  {
00685     par1 = 1./currentRange ;
00686     par2 = 1./(par1*lambda0) ;
00687     par3 = 1.+par2 ;
00688     if(tPathLength < currentRange)
00689       zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
00690     else
00691       zmean = 1./(par1*par3) ;
00692   } else {
00693     G4double T1 = GetEnergy(particle,currentRange-tPathLength,couple);
00694     G4double lambda1 = GetTransportMeanFreePath(particle,T1);
00695 
00696     par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
00697     par2 = 1./(par1*lambda0) ;
00698     par3 = 1.+par2 ;
00699     zmean = (1.-exp(par3*log(lambda1/lambda0)))/(par1*par3) ;
00700   }
00701 
00702   zPathLength = zmean ;
00703 
00704   //  sample z
00705   if(samplez)
00706   {
00707     const G4double  ztmax = 0.99 ;
00708     G4double zt = zmean/tPathLength ;
00709 
00710     if (tPathLength > stepmin && zt < ztmax)              
00711     {
00712       G4double u,cz1;
00713       if(zt >= third)
00714       {
00715         G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
00716         cz1 = 1.+cz ;
00717         G4double u0 = cz/cz1 ;
00718         G4double grej ;
00719         do {
00720             u = exp(log(G4UniformRand())/cz1) ;
00721             grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
00722            } while (grej < G4UniformRand()) ;
00723       }
00724       else
00725       {
00726         cz1 = 1./zt-1.;
00727         u = 1.-exp(log(G4UniformRand())/cz1) ;
00728       }
00729       zPathLength = tPathLength*u ;
00730     }
00731   }
00732 
00733   if(zPathLength > lambda0) zPathLength = lambda0;
00734 
00735   return zPathLength;
00736 }

G4double G4UrbanMscModel92::ComputeTheta0 ( G4double  truePathLength,
G4double  KineticEnergy 
)

Definition at line 771 of file G4UrbanMscModel92.cc.

00773 {
00774   // for all particles take the width of the central part
00775   //  from a  parametrization similar to the Highland formula
00776   // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
00777   const G4double c_highland = 13.6*MeV ;
00778   G4double betacp = sqrt(currentKinEnergy*(currentKinEnergy+2.*mass)*
00779                          KineticEnergy*(KineticEnergy+2.*mass)/
00780                       ((currentKinEnergy+mass)*(KineticEnergy+mass)));
00781   y = trueStepLength/currentRadLength;
00782   G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp;
00783   y = log(y);
00784   // correction factor from e-/proton scattering data
00785   G4double corr = coeffth1+coeffth2*y;                
00786   if(y < -6.5) corr -= 0.011*(6.5+y);
00787   theta0 *= corr ;                                               
00788 
00789   return theta0;
00790 }

G4double G4UrbanMscModel92::ComputeTruePathLengthLimit ( const G4Track &  track,
G4double &  currentMinimalStep 
) [virtual]

Reimplemented from G4VMscModel.

Definition at line 462 of file G4UrbanMscModel92.cc.

References G4VMscModel::ComputeGeomLimit(), G4VMscModel::ComputeSafety(), G4VMscModel::ConvertTrueToGeom(), G4VMscModel::facgeom, G4VMscModel::facrange, G4VMscModel::facsafety, fGeomBoundary, fUseDistanceToBoundary, fUseSafety, G4Track::GetDynamicParticle(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterialCutsCouple(), G4Step::GetPreStepPoint(), G4VMscModel::GetRange(), G4Track::GetStep(), G4VMscModel::GetTransportMeanFreePath(), G4VEmModel::SetCurrentCouple(), G4VMscModel::skin, G4InuclParticleNames::sp, and G4VMscModel::steppingAlgorithm.

00465 {
00466   tPathLength = currentMinimalStep;
00467   const G4DynamicParticle* dp = track.GetDynamicParticle();
00468   G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
00469   G4StepStatus stepStatus = sp->GetStepStatus();
00470   couple = track.GetMaterialCutsCouple();
00471   SetCurrentCouple(couple); 
00472   currentMaterialIndex = couple->GetIndex();
00473   currentKinEnergy = dp->GetKineticEnergy();
00474   currentRange = GetRange(particle,currentKinEnergy,couple);
00475   lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy);
00476 
00477   // stop here if small range particle
00478   if(inside || tPathLength < tlimitminfix) {
00479     return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00480   }
00481   
00482   if(tPathLength > currentRange) tPathLength = currentRange;
00483 
00484   presafety = sp->GetSafety();
00485 
00486   //    G4cout << "G4UrbanMscModel92::ComputeTruePathLengthLimit tPathLength= " 
00487   //       <<tPathLength<<" safety= " << presafety
00488   //       << " range= " <<currentRange<<G4endl;
00489 
00490   // far from geometry boundary
00491   if(currentRange < presafety)
00492     {
00493       inside = true;
00494       return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
00495     }
00496 
00497   // standard  version
00498   //
00499   if (steppingAlgorithm == fUseDistanceToBoundary)
00500     {
00501       //compute geomlimit and presafety 
00502       geomlimit = ComputeGeomLimit(track, presafety, currentRange);
00503 
00504       // is it far from boundary ?
00505       if(currentRange < presafety)
00506         {
00507           inside = true;
00508           return ConvertTrueToGeom(tPathLength, currentMinimalStep);   
00509         }
00510 
00511       smallstep += 1.;
00512       insideskin = false;
00513 
00514       if(firstStep || stepStatus == fGeomBoundary)
00515         {
00516           rangeinit = currentRange;
00517           if(firstStep) smallstep = 1.e10;
00518           else  smallstep = 1.;
00519 
00520           // constraint from the geometry 
00521           if((geomlimit < geombig) && (geomlimit > geommin))
00522             {
00523               if(stepStatus == fGeomBoundary)  
00524                 tgeom = geomlimit/facgeom;
00525               else
00526                 tgeom = 2.*geomlimit/facgeom;
00527             }
00528             else
00529               tgeom = geombig;
00530 
00531           //define stepmin here (it depends on lambda!)
00532           //rough estimation of lambda_elastic/lambda_transport
00533           G4double rat = currentKinEnergy/MeV ;
00534           rat = 1.e-3/(rat*(10.+rat)) ;
00535           //stepmin ~ lambda_elastic
00536           stepmin = rat*lambda0;
00537           skindepth = skin*stepmin;
00538 
00539           //define tlimitmin
00540           tlimitmin = 10.*stepmin;
00541           if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
00542 
00543         }
00544 
00545       //step limit 
00546       tlimit = facrange*rangeinit;              
00547       if(tlimit < facsafety*presafety)
00548         tlimit = facsafety*presafety; 
00549 
00550       //lower limit for tlimit
00551       if(tlimit < tlimitmin) tlimit = tlimitmin;
00552 
00553       if(tlimit > tgeom) tlimit = tgeom;
00554 
00555       //  G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit  
00556       //     << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
00557 
00558       // shortcut
00559       if((tPathLength < tlimit) && (tPathLength < presafety) &&
00560          (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
00561         return ConvertTrueToGeom(tPathLength, currentMinimalStep);   
00562 
00563       // step reduction near to boundary
00564       if(smallstep < skin)
00565         {
00566           tlimit = stepmin;
00567           insideskin = true;
00568         }
00569       else if(geomlimit < geombig)
00570         {
00571           if(geomlimit > skindepth)
00572             {
00573               if(tlimit > geomlimit-0.999*skindepth)
00574                 tlimit = geomlimit-0.999*skindepth;
00575             }
00576           else
00577             {
00578               insideskin = true;
00579               if(tlimit > stepmin) tlimit = stepmin;
00580             }
00581         }
00582 
00583       if(tlimit < stepmin) tlimit = stepmin;
00584 
00585       if(tPathLength > tlimit) tPathLength = tlimit  ; 
00586 
00587     }
00588     // for 'normal' simulation with or without magnetic field 
00589     //  there no small step/single scattering at boundaries
00590   else if(steppingAlgorithm == fUseSafety)
00591     {
00592       // compute presafety again if presafety <= 0 and no boundary
00593       // i.e. when it is needed for optimization purposes
00594       if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 
00595         presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
00596 
00597       // is far from boundary
00598       if(currentRange < presafety)
00599         {
00600           inside = true;
00601           return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
00602         }
00603 
00604       if(firstStep || stepStatus == fGeomBoundary)
00605       {
00606         rangeinit = currentRange;
00607         fr = facrange;
00608         // 9.1 like stepping for e+/e- only (not for muons,hadrons)
00609         if(mass < masslimite) 
00610         {
00611           if(lambda0 > currentRange)
00612             rangeinit = lambda0;
00613           if(lambda0 > lambdalimit)
00614             fr *= 0.75+0.25*lambda0/lambdalimit;
00615         }
00616 
00617         //lower limit for tlimit
00618         G4double rat = currentKinEnergy/MeV ;
00619         rat = 1.e-3/(rat*(10.+rat)) ;
00620         tlimitmin = 10.*lambda0*rat;
00621         if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
00622       }
00623       //step limit
00624       tlimit = fr*rangeinit;               
00625 
00626       if(tlimit < facsafety*presafety)
00627         tlimit = facsafety*presafety;
00628 
00629       //lower limit for tlimit
00630       if(tlimit < tlimitmin) tlimit = tlimitmin;
00631 
00632       if(tPathLength > tlimit) tPathLength = tlimit;
00633     }
00634   
00635   // version similar to 7.1 (needed for some experiments)
00636   else
00637     {
00638       if (stepStatus == fGeomBoundary)
00639         {
00640           if (currentRange > lambda0) tlimit = facrange*currentRange;
00641           else                        tlimit = facrange*lambda0;
00642 
00643           if(tlimit < tlimitmin) tlimit = tlimitmin;
00644           if(tPathLength > tlimit) tPathLength = tlimit;
00645         }
00646     }
00647   // G4cout << "tPathLength= " << tPathLength << "  geomlimit= " << geomlimit 
00648   //     << " currentMinimalStep= " << currentMinimalStep << G4endl;
00649   return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00650 }

G4double G4UrbanMscModel92::ComputeTrueStepLength ( G4double  geomStepLength  )  [virtual]

Reimplemented from G4VMscModel.

Definition at line 740 of file G4UrbanMscModel92.cc.

00741 {
00742   // step defined other than transportation 
00743   if(geomStepLength == zPathLength && tPathLength <= currentRange)
00744     return tPathLength;
00745 
00746   // t = z for very small step
00747   zPathLength = geomStepLength;
00748   tPathLength = geomStepLength;
00749   if(geomStepLength < tlimitminfix) return tPathLength;
00750   
00751   // recalculation
00752   if((geomStepLength > lambda0*tausmall) && !insideskin)
00753   {
00754     if(par1 <  0.)
00755       tPathLength = -lambda0*log(1.-geomStepLength/lambda0) ;
00756     else 
00757     {
00758       if(par1*par3*geomStepLength < 1.)
00759         tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
00760       else 
00761         tPathLength = currentRange;
00762     }  
00763   }
00764   if(tPathLength < geomStepLength) tPathLength = geomStepLength;
00765 
00766   return tPathLength;
00767 }

void G4UrbanMscModel92::Initialise ( const G4ParticleDefinition *  ,
const G4DataVector &   
) [virtual]

Implements G4VEmModel.

Definition at line 190 of file G4UrbanMscModel92.cc.

References G4cout, G4endl, G4VMscModel::GetParticleChangeForMSC(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGMass(), and G4VMscModel::skin.

00192 {
00193   skindepth = skin*stepmin;
00194   // set values of some data members
00195   SetParticle(p);
00196 
00197   if(p->GetPDGMass() > MeV) {
00198     G4cout << "### WARNING: G4UrbanMscModel92 model is used for " 
00199            << p->GetParticleName() << " !!! " << G4endl;
00200     G4cout << "###          This model should be used only for e+-" 
00201            << G4endl;
00202   }
00203 
00204   fParticleChange = GetParticleChangeForMSC(p);
00205 }

G4ThreeVector & G4UrbanMscModel92::SampleScattering ( const G4ThreeVector &  ,
G4double  safety 
) [virtual]

Reimplemented from G4VMscModel.

Definition at line 795 of file G4UrbanMscModel92.cc.

References G4VEmModel::CurrentCouple(), G4VMscModel::dtrl, G4VMscModel::fDisplacement, G4endl, G4Exception(), G4UniformRand, G4VMscModel::GetDEDX(), G4VMscModel::GetEnergy(), G4MaterialCutsCouple::GetMaterial(), G4Material::GetName(), G4ParticleDefinition::GetParticleName(), JustWarning, G4VMscModel::latDisplasment, and G4ParticleChangeForMSC::ProposeMomentumDirection().

00797 {
00798   fDisplacement.set(0.0,0.0,0.0);
00799   G4double kineticEnergy = currentKinEnergy;
00800   if (tPathLength > currentRange*dtrl) {
00801     kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
00802   } else {
00803     kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple);
00804   }
00805   if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix) ||
00806      (tPathLength/tausmall < lambda0)) return fDisplacement;
00807 
00808   G4double cth  = SampleCosineTheta(tPathLength,kineticEnergy);
00809 
00810   // protection against 'bad' cth values
00811   if(std::abs(cth) > 1.) return fDisplacement;
00812 
00813   // extra protection agaist high energy particles backscattered 
00814   if(cth < 1.0 - 1000*tPathLength/lambda0 && kineticEnergy > 20*MeV) { 
00815     //G4cout << "Warning: large scattering E(MeV)= " << kineticEnergy 
00816     //     << " s(mm)= " << tPathLength/mm
00817     //     << " 1-cosTheta= " << 1.0 - cth << G4endl;
00818     // do Gaussian central scattering
00819     if(kineticEnergy > GeV && cth < 0.0) {
00820       G4ExceptionDescription ed;
00821       ed << particle->GetParticleName()
00822          << " E(MeV)= " << kineticEnergy/MeV
00823          << " Step(mm)= " << tPathLength/mm
00824          << " in " << CurrentCouple()->GetMaterial()->GetName()
00825          << " CosTheta= " << cth 
00826          << " is too big - the angle is resampled" << G4endl;
00827       G4Exception("G4UrbanMscModel92::SampleScattering","em0004",
00828                   JustWarning, ed,"");
00829     }
00830     do {
00831       cth = 1.0 + 2*log(G4UniformRand())*tPathLength/lambda0;
00832     } while(cth < -1.0);
00833   }
00834 
00835   G4double sth  = sqrt((1.0 - cth)*(1.0 + cth));
00836   G4double phi  = twopi*G4UniformRand();
00837   G4double dirx = sth*cos(phi);
00838   G4double diry = sth*sin(phi);
00839 
00840   // G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
00841   G4ThreeVector newDirection(dirx,diry,cth);
00842   newDirection.rotateUz(oldDirection);
00843   fParticleChange->ProposeMomentumDirection(newDirection);
00844 
00845   if (latDisplasment && safety > tlimitminfix) {
00846 
00847     G4double r = SampleDisplacement();
00848     /*
00849     G4cout << "G4UrbanMscModel92::SampleSecondaries: e(MeV)= " << kineticEnergy
00850            << " sinTheta= " << sth << " r(mm)= " << r
00851            << " trueStep(mm)= " << tPathLength
00852            << " geomStep(mm)= " << zPathLength
00853            << G4endl;
00854     */
00855     if(r > 0.)
00856       {
00857         G4double latcorr = LatCorrelation();
00858         if(latcorr > r) latcorr = r;
00859 
00860         // sample direction of lateral displacement
00861         // compute it from the lateral correlation
00862         G4double Phi = 0.;
00863         if(std::abs(r*sth) < latcorr)
00864           Phi  = twopi*G4UniformRand();
00865         else
00866         {
00867           G4double psi = std::acos(latcorr/(r*sth));
00868           if(G4UniformRand() < 0.5)
00869             Phi = phi+psi;
00870           else
00871             Phi = phi-psi;
00872         }
00873 
00874         dirx = r*std::cos(Phi);
00875         diry = r*std::sin(Phi);
00876 
00877         fDisplacement.set(dirx,diry,0.0);
00878         fDisplacement.rotateUz(oldDirection);
00879       }
00880   }
00881   return fDisplacement;
00882 }

void G4UrbanMscModel92::StartTracking ( G4Track *   )  [virtual]

Reimplemented from G4VEmModel.

Definition at line 451 of file G4UrbanMscModel92.cc.

References G4DynamicParticle::GetDefinition(), and G4Track::GetDynamicParticle().

00452 {
00453   SetParticle(track->GetDynamicParticle()->GetDefinition());
00454   firstStep = true; 
00455   inside = false;
00456   insideskin = false;
00457   tlimit = geombig;
00458 }

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