G4UrbanMscModel92.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 //
00027 // $Id: G4UrbanMscModel92.cc 66592 2012-12-23 09:34:55Z vnivanch $
00028 //
00029 // -------------------------------------------------------------------
00030 //
00031 // GEANT4 Class file
00032 //
00033 //
00034 // File name:   G4UrbanMscModel92
00035 //
00036 // Author:      Laszlo Urban
00037 //
00038 // Creation date: 06.03.2008
00039 //
00040 // Modifications:
00041 //
00042 // 06-03-2008 starting point : G4UrbanMscModel2 = G4UrbanMscModel 9.1 ref 02
00043 //
00044 // 13-03-08  Bug in SampleScattering (which caused lateral asymmetry) fixed
00045 //           (L.Urban)
00046 //
00047 // 14-03-08  Simplification of step limitation in ComputeTruePathLengthLimit,
00048 //           + tlimitmin is the same for UseDistancetoBoundary and
00049 //           UseSafety (L.Urban)           
00050 //
00051 // 16-03-08  Reorganization of SampleCosineTheta + new method SimpleScattering
00052 //           SimpleScattering is used if the relative energy loss is too big
00053 //           or theta0 is too big (see data members rellossmax, theta0max)
00054 //           (L.Urban)          
00055 //
00056 // 17-03-08  tuning of the correction factor in ComputeTheta0 (L.Urban)
00057 //
00058 // 19-03-08  exponent c of the 'tail' model function is not equal to 2 any more,
00059 //           value of c has been extracted from some e- scattering data (L.Urban) 
00060 //
00061 // 24-03-08  Step limitation in ComputeTruePathLengthLimit has been
00062 //           simplified further + some data members have been removed (L.Urban)
00063 //
00064 // 24-07-08  central part of scattering angle (theta0) has been tuned
00065 //           tail of the scattering angle distribution has been tuned
00066 //           using some e- and proton scattering data
00067 //
00068 // 05-08-08  bugfix in ComputeTruePathLengthLimit (L.Urban)
00069 //
00070 // 09-10-08  theta0 and tail have been retuned using some e-,mu,proton
00071 //           scattering data (L.Urban)
00072 //           + single scattering without path length correction for
00073 //           small steps (t < tlimitmin, for UseDistanceToBoundary only)
00074 //
00075 // 15-10-08  Moliere-Bethe screening in the single scattering part(L.Urban)          
00076 //
00077 // 17-10-08  stepping similar to that in model (9.1) for UseSafety case
00078 //           for e+/e- in order to speed up the code for calorimeters
00079 //
00080 // 23-10-08  bugfix in the screeningparameter of the single scattering part,
00081 //           some technical change in order to speed up the code (UpdateCache)
00082 //
00083 // 27-10-08  bugfix in ComputeTruePathLengthLimit (affects UseDistanceToBoundary
00084 //           stepping type only) (L.Urban)          
00085 //
00086 // 28-10-09  V.Ivanchenko moved G4UrbanMscModel to G4UrbanMscModel92, 
00087 //           now it is a frozen version of the Urban model corresponding 
00088 //           to g4 9.2
00089 
00090 // Class Description:
00091 //
00092 // Implementation of the model of multiple scattering based on
00093 // H.W.Lewis Phys Rev 78 (1950) 526 and others
00094 
00095 // -------------------------------------------------------------------
00096 // In its present form the model can be  used for simulation 
00097 //   of the e-/e+, muon and charged hadron multiple scattering
00098 //
00099 
00100 
00101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00103 
00104 #include "G4UrbanMscModel92.hh"
00105 #include "G4PhysicalConstants.hh"
00106 #include "G4SystemOfUnits.hh"
00107 #include "Randomize.hh"
00108 #include "G4Electron.hh"
00109 #include "G4LossTableManager.hh"
00110 #include "G4ParticleChangeForMSC.hh"
00111 
00112 #include "G4Poisson.hh"
00113 #include "globals.hh"
00114 
00115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00116 
00117 using namespace std;
00118 
00119 G4UrbanMscModel92::G4UrbanMscModel92(const G4String& nam)
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 }
00182 
00183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00184 
00185 G4UrbanMscModel92::~G4UrbanMscModel92()
00186 {}
00187 
00188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00189 
00190 void G4UrbanMscModel92::Initialise(const G4ParticleDefinition* p,
00191                                    const G4DataVector&)
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 }
00206 
00207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00208 
00209 G4double G4UrbanMscModel92::ComputeCrossSectionPerAtom( 
00210                              const G4ParticleDefinition* part,
00211                                    G4double KineticEnergy,
00212                                    G4double AtomicNumber,G4double,
00213                                    G4double, G4double)
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 }
00448 
00449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00450 
00451 void G4UrbanMscModel92::StartTracking(G4Track* track)
00452 {
00453   SetParticle(track->GetDynamicParticle()->GetDefinition());
00454   firstStep = true; 
00455   inside = false;
00456   insideskin = false;
00457   tlimit = geombig;
00458 }
00459 
00460 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00461 
00462 G4double G4UrbanMscModel92::ComputeTruePathLengthLimit(
00463                              const G4Track& track,
00464                              G4double& currentMinimalStep)
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 }
00651 
00652 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00653 
00654 G4double G4UrbanMscModel92::ComputeGeomPathLength(G4double)
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 }
00737 
00738 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00739 
00740 G4double G4UrbanMscModel92::ComputeTrueStepLength(G4double geomStepLength)
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 }
00768 
00769 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00770 
00771 G4double G4UrbanMscModel92::ComputeTheta0(G4double trueStepLength,
00772                                           G4double KineticEnergy)
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 }
00791 
00792 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00793 
00794 G4ThreeVector& 
00795 G4UrbanMscModel92::SampleScattering(const G4ThreeVector& oldDirection,
00796                                     G4double safety)
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 }
00883 
00884 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00885 
00886 G4double G4UrbanMscModel92::SampleCosineTheta(G4double trueStepLength,
00887                                               G4double KineticEnergy)
00888 {
00889   G4double cth = 1. ;
00890   G4double tau = trueStepLength/lambda0 ;
00891 
00892   Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
00893          couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ;
00894 
00895   if(Zold != Zeff)  
00896     UpdateCache();
00897 
00898   if(insideskin)
00899   {
00900     //no scattering, single or plural scattering
00901     G4double mean = trueStepLength/stepmin ;
00902 
00903     G4int n = G4Poisson(mean);
00904     if(n > 0)
00905     {
00906       //screening (Moliere-Bethe)
00907       G4double mom2 = KineticEnergy*(2.*mass+KineticEnergy);
00908       G4double beta2 = mom2/((KineticEnergy+mass)*(KineticEnergy+mass));
00909       G4double ascr = scr1/mom2;
00910       ascr *= 1.13+scr2/beta2;
00911       G4double ascr1 = 1.+2.*ascr;
00912       G4double bp1=ascr1+1.;
00913       G4double bm1=ascr1-1.;
00914 
00915       // single scattering from screened Rutherford x-section
00916       G4double ct,st,phi;
00917       G4double sx=0.,sy=0.,sz=0.;
00918       for(G4int i=1; i<=n; i++)
00919       {
00920         ct = ascr1-bp1*bm1/(2.*G4UniformRand()+bm1);
00921         if(ct < -1.) ct = -1.;
00922         if(ct >  1.) ct =  1.; 
00923         st = sqrt(1.-ct*ct);
00924         phi = twopi*G4UniformRand();
00925         sx += st*cos(phi);
00926         sy += st*sin(phi);
00927         sz += ct;
00928       }
00929       cth = sz/sqrt(sx*sx+sy*sy+sz*sz);
00930     }
00931   }
00932   else
00933   {
00934     if(trueStepLength >= currentRange*dtrl) 
00935     {
00936       if(par1*trueStepLength < 1.)
00937         tau = -par2*log(1.-par1*trueStepLength) ;
00938       // for the case if ioni/brems are inactivated
00939       // see the corresponding condition in ComputeGeomPathLength 
00940       else if(1.-KineticEnergy/currentKinEnergy > taulim)
00941         tau = taubig ;
00942     }
00943     currentTau = tau ;
00944     lambdaeff = trueStepLength/currentTau;
00945     currentRadLength = couple->GetMaterial()->GetRadlen();
00946 
00947     if (tau >= taubig) cth = -1.+2.*G4UniformRand();
00948     else if (tau >= tausmall)
00949     {
00950       G4double xsi = 3.0;  
00951       G4double x0 = 1.;
00952       G4double a = 1., ea = 0., eaa = 1.;
00953       G4double b=2.,b1=3.,bx=1.,eb1=3.,ebx=1.;
00954       G4double prob = 1. , qprob = 1. ;
00955       G4double xmean1 = 1., xmean2 = 0.;
00956       G4double xmeanth = exp(-tau);
00957       G4double x2meanth = (1.+2.*exp(-2.5*tau))/3.;
00958 
00959       G4double relloss = 1.-KineticEnergy/currentKinEnergy;
00960       if(relloss > rellossmax) 
00961         return SimpleScattering(xmeanth,x2meanth);
00962 
00963       G4double theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
00964 
00965       // protection for very small angles
00966       if(theta0*theta0 < tausmall) return cth;
00967     
00968       if(theta0 > theta0max)
00969         return SimpleScattering(xmeanth,x2meanth);
00970       G4double sth = sin(0.5*theta0);
00971       a = 0.25/(sth*sth);
00972 
00973       ea = exp(-xsi);
00974       eaa = 1.-ea ;
00975       xmean1 = 1.-(1.-(1.+xsi)*ea)/(a*eaa);
00976       x0 = 1.-xsi/a;
00977 
00978       if(xmean1 <= 0.999*xmeanth)
00979         return SimpleScattering(xmeanth,x2meanth);
00980 
00981       // from MUSCAT H,Be,Fe data 
00982       G4double c = coeffc1;                         
00983       if(y > -13.5) 
00984         c += coeffc2*exp(3.*log(y+13.5));
00985 
00986       if(abs(c-3.) < 0.001)  c = 3.001;      
00987       if(abs(c-2.) < 0.001)  c = 2.001;      
00988 
00989       G4double c1 = c-1.;
00990 
00991       //from continuity of derivatives
00992       b = 1.+(c-xsi)/a;
00993 
00994       b1 = b+1.;
00995       bx = c/a;
00996       eb1 = exp(c1*log(b1));
00997       ebx = exp(c1*log(bx));
00998 
00999       xmean2 = (x0*eb1+ebx-(eb1*bx-b1*ebx)/(c-2.))/(eb1-ebx);
01000       
01001       G4double f1x0 = a*ea/eaa;
01002       G4double f2x0 = c1*eb1/(bx*(eb1-ebx));
01003       prob = f2x0/(f1x0+f2x0);
01004 
01005       qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
01006 
01007       // sampling of costheta
01008       if(G4UniformRand() < qprob)
01009       {
01010         if(G4UniformRand() < prob)
01011           cth = 1.+log(ea+G4UniformRand()*eaa)/a ;
01012         else
01013           cth = b-b1*bx/exp(log(ebx+(eb1-ebx)*G4UniformRand())/c1) ;
01014       }
01015       else
01016         cth = -1.+2.*G4UniformRand();
01017     }
01018   }  
01019   return cth ;
01020 }
01021 
01022 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01023 
01024 G4double G4UrbanMscModel92::SimpleScattering(G4double xmeanth,G4double x2meanth)
01025 {
01026   // 'large angle scattering'
01027   // 2 model functions with correct xmean and x2mean
01028   G4double a = (2.*xmeanth+9.*x2meanth-3.)/(2.*xmeanth-3.*x2meanth+1.);
01029   G4double prob = (a+2.)*xmeanth/a;
01030 
01031   // sampling
01032   G4double cth = 1.;
01033   if(G4UniformRand() < prob)
01034     cth = -1.+2.*exp(log(G4UniformRand())/(a+1.));
01035   else
01036     cth = -1.+2.*G4UniformRand();
01037   return cth;
01038 }
01039 
01040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01041 
01042 G4double G4UrbanMscModel92::SampleDisplacement()
01043 {
01044   const G4double kappa = 2.5;
01045   const G4double kappapl1 = kappa+1.;
01046   const G4double kappami1 = kappa-1.;
01047   G4double rmean = 0.0;
01048   if ((currentTau >= tausmall) && !insideskin) {
01049     if (currentTau < taulim) {
01050       rmean = kappa*currentTau*currentTau*currentTau*
01051              (1.-kappapl1*currentTau*0.25)/6. ;
01052 
01053     } else {
01054       G4double etau = 0.0;
01055       if (currentTau<taubig) etau = exp(-currentTau);
01056       rmean = -kappa*currentTau;
01057       rmean = -exp(rmean)/(kappa*kappami1);
01058       rmean += currentTau-kappapl1/kappa+kappa*etau/kappami1;
01059     }
01060     if (rmean>0.) rmean = 2.*lambdaeff*sqrt(rmean/3.0);
01061     else          rmean = 0.;
01062   }
01063 
01064   // protection against z > t ...........................
01065   if(rmean > 0.) {
01066     G4double zt = (tPathLength-zPathLength)*(tPathLength+zPathLength);
01067     if(zt <= 0.)
01068       rmean = 0.;
01069     else if(rmean*rmean > zt)
01070       rmean = sqrt(zt);
01071   }
01072   return rmean;
01073 }
01074 
01075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
01076 
01077 G4double G4UrbanMscModel92::LatCorrelation()
01078 {
01079   const G4double kappa = 2.5;
01080   const G4double kappami1 = kappa-1.;
01081 
01082   G4double latcorr = 0.;
01083   if((currentTau >= tausmall) && !insideskin)
01084   {
01085     if(currentTau < taulim)
01086       latcorr = lambdaeff*kappa*currentTau*currentTau*
01087                 (1.-(kappa+1.)*currentTau/3.)/3.;
01088     else
01089     {
01090       G4double etau = 0.;
01091       if(currentTau < taubig) etau = exp(-currentTau);
01092       latcorr = -kappa*currentTau;
01093       latcorr = exp(latcorr)/kappami1;
01094       latcorr += 1.-kappa*etau/kappami1 ;
01095       latcorr *= 2.*lambdaeff/3. ;
01096     }
01097   }
01098 
01099   return latcorr;
01100 }
01101 
01102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

Generated on Mon May 27 17:50:08 2013 for Geant4 by  doxygen 1.4.7