80 CLHEP::twopi*CLHEP::classic_electr_radius*CLHEP::classic_electr_radius;
82 CLHEP::electron_mass_c2*CLHEP::Bohr_radius*CLHEP::Bohr_radius
83 /(CLHEP::hbarc*CLHEP::hbarc);
85 ((
Tlim+CLHEP::electron_mass_c2)*(
Tlim+CLHEP::electron_mass_c2));
87 (CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
171 mass = proton_mass_c2;
216 static const G4double epsmin = 1.e-4 , epsmax = 1.e10;
218 static const G4double Zdat[15] = { 4., 6., 13., 20., 26., 29., 32., 38.,47.,
219 50., 56., 64., 74., 79., 82. };
222 static const G4double celectron[15][22] =
223 {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
224 1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
225 1.112,1.108,1.100,1.093,1.089,1.087 },
226 {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
227 1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
228 1.109,1.105,1.097,1.090,1.086,1.082 },
229 {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
230 1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
231 1.131,1.124,1.113,1.104,1.099,1.098 },
232 {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
233 1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
234 1.112,1.105,1.096,1.089,1.085,1.098 },
235 {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
236 1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
237 1.073,1.070,1.064,1.059,1.056,1.056 },
238 {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
239 1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
240 1.074,1.070,1.063,1.059,1.056,1.052 },
241 {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
242 1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
243 1.068,1.064,1.059,1.054,1.051,1.050 },
244 {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
245 1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
246 1.039,1.037,1.034,1.031,1.030,1.036 },
247 {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
248 1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
249 1.031,1.028,1.024,1.022,1.021,1.024 },
250 {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
251 1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
252 1.020,1.017,1.015,1.013,1.013,1.020 },
253 {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
254 1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
255 0.995,0.993,0.993,0.993,0.993,1.011 },
256 {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
257 1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
258 0.974,0.972,0.973,0.974,0.975,0.987 },
259 {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
260 1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
261 0.950,0.947,0.949,0.952,0.954,0.963 },
262 {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
263 1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
264 0.941,0.938,0.940,0.944,0.946,0.954 },
265 {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
266 1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
267 0.933,0.930,0.933,0.936,0.939,0.949 }};
269 static const G4double cpositron[15][22] = {
270 {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
271 1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
272 1.131,1.126,1.117,1.108,1.103,1.100 },
273 {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
274 1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
275 1.138,1.132,1.122,1.113,1.108,1.102 },
276 {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
277 1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
278 1.203,1.190,1.173,1.159,1.151,1.145 },
279 {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
280 1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
281 1.225,1.210,1.191,1.175,1.166,1.174 },
282 {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
283 1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
284 1.217,1.203,1.184,1.169,1.160,1.151 },
285 {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
286 1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
287 1.237,1.222,1.201,1.184,1.174,1.159 },
288 {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
289 1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
290 1.252,1.234,1.212,1.194,1.183,1.170 },
291 {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
292 2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
293 1.254,1.237,1.214,1.195,1.185,1.179 },
294 {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
295 2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
296 1.312,1.288,1.258,1.235,1.221,1.205 },
297 {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
298 2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
299 1.320,1.294,1.264,1.240,1.226,1.214 },
300 {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
301 2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
302 1.328,1.302,1.270,1.245,1.231,1.233 },
303 {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
304 2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
305 1.361,1.330,1.294,1.267,1.251,1.239 },
306 {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
307 3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
308 1.409,1.372,1.330,1.298,1.280,1.258 },
309 {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
310 3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
311 1.442,1.400,1.354,1.319,1.299,1.272 },
312 {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
313 3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
314 1.456,1.412,1.364,1.328,1.307,1.282 }};
318 static const G4double hecorr[15] = {
319 120.70, 117.50, 105.00, 92.92, 79.23, 74.510, 68.29,
320 57.39, 41.97, 36.14, 24.53, 10.21, -7.855, -16.84,
332 G4double eKineticEnergy = KineticEnergy;
334 if(
mass > electron_mass_c2)
337 G4double c =
mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
339 G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
340 eKineticEnergy = electron_mass_c2*tau ;
343 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
344 G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
345 /(eTotalEnergy*eTotalEnergy);
346 G4double bg2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
347 /(electron_mass_c2*electron_mass_c2);
351 if (eps<epsmin) sigma = 2.*eps*
eps;
352 else if(eps<epsmax) sigma =
G4Log(1.+2.*eps)-2.*eps/(1.+2.*
eps);
353 else sigma =
G4Log(2.*eps)-1.+1./
eps;
355 sigma *=
ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
363 while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
369 G4double ratZ = (AtomicNumber-ZZ1)*(AtomicNumber+ZZ1)/
370 ((ZZ2-ZZ1)*(ZZ2+ZZ1));
372 if(eKineticEnergy <=
Tlim)
377 while ((iT>=0)&&(
Tdat[iT]>=eKineticEnergy)) iT -= 1;
383 G4double b2small = T*(E+electron_mass_c2)/(E*E);
385 T =
Tdat[iT+1]; E = T + electron_mass_c2;
386 G4double b2big = T*(E+electron_mass_c2)/(E*E);
387 G4double ratb2 = (beta2-b2small)/(b2big-b2small);
391 c1 = celectron[iZ][iT];
392 c2 = celectron[iZ+1][iT];
393 cc1 = c1+ratZ*(c2-
c1);
395 c1 = celectron[iZ][iT+1];
396 c2 = celectron[iZ+1][iT+1];
397 cc2 = c1+ratZ*(c2-
c1);
399 corr = cc1+ratb2*(cc2-cc1);
405 c1 = cpositron[iZ][iT];
406 c2 = cpositron[iZ+1][iT];
407 cc1 = c1+ratZ*(c2-
c1);
409 c1 = cpositron[iZ][iT+1];
410 c2 = cpositron[iZ+1][iT+1];
411 cc2 = c1+ratZ*(c2-
c1);
413 corr = cc1+ratb2*(cc2-cc1);
422 if((AtomicNumber >= ZZ1) && (AtomicNumber <= ZZ2))
423 sigma = c1+ratZ*(c2-
c1) ;
424 else if(AtomicNumber < ZZ1)
425 sigma = AtomicNumber*AtomicNumber*c1/(ZZ1*ZZ1);
426 else if(AtomicNumber > ZZ2)
427 sigma = AtomicNumber*AtomicNumber*c2/(ZZ2*ZZ2);
489 distance *= (1.15-9.76e-4*
Zeff);
491 distance *= (1.20-
Zeff*(1.62e-2-9.22e-5*
Zeff));
538 rat = 1.e-3/(rat*(10.+rat)) ;
648 rat = 1.e-3/(rat*(10 + rat)) ;
691 if(
charge > 0.) index = 2;
699 rat = 1.e-3/(rat*(10 + rat)) ;
849 if(tlength < geomStepLength) { tlength = geomStepLength; }
893 G4double sth = sqrt((1.0 - cth)*(1.0 + cth));
899 newDirection.rotateUz(oldDirection);
942 static const G4double numlim = 0.01;
945 xmeanth = 1.0 - tau*(1.0 - 0.5*tau);
946 x2meanth= 1.0 - tau*(5.0 - 6.25*tau)/3.;
948 xmeanth =
G4Exp(-tau);
949 x2meanth = (1.+2.*
G4Exp(-2.5*tau))/3.;
958 G4bool extremesmallstep = false ;
961 if(trueStepLength > tsmall) {
964 theta0 = sqrt(trueStepLength/tsmall)*
ComputeTheta0(tsmall,KineticEnergy);
965 extremesmallstep = true ;
974 if(theta2 <
tausmall) {
return cth; }
981 if(theta2 > numlim) {
1010 if(fabs(c-3.) < 0.001) { c = 3.001; }
1011 else if(fabs(c-2.) < 0.001) { c = 2.001; }
1017 G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/eaa;
1022 if(xmean1 <= 0.999*xmeanth) {
1035 G4double xmean2 = (x0 + d - (bx - b1*d)/(c-2.))/(1. - d);
1041 G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
1054 if(var < numlim*d) {
1056 cth = -1.0 + var*(1.0 - 0.5*var*c)*(2. + (c - xsi)*
x);
1058 cth = 1. + x*(c - xsi - c*
G4Exp(-
G4Log(var + d)/c1));
1102 KineticEnergy*(KineticEnergy+2.*
mass)));
1113 G4double x = std::sqrt(tau*(tau+2.)/((tau+1.)*(tau+1.)));
1119 corr = a*(1.-
G4Exp(-b*x));
1121 corr = c+d*
G4Exp(e*(x-1.));
1130 y *= corr*(1.+
Zeff*(1.84035e-4*
Zeff-1.86427e-2)+0.41125);
1156 static const G4double kappami1 = 1.5;
1168 latcorr =
G4Exp(latcorr)/kappami1;
1169 latcorr += 1.-kappa*etau/kappami1 ;
1178 if(std::abs(r*sth) < latcorr) {
1184 G4double psi = std::acos(latcorr/(r*sth));
1224 static const G4double probv1 = 0.305533;
1225 static const G4double probv2 = 0.955176;
1226 static const G4double vhigh = 3.15;
1232 if(random < probv1) {
1237 while (std::abs(v) >= vhigh);
1242 if(random < probv2) {
1249 if(random < 0.5) { Phi = phi+v; }
1250 else { Phi = phi-v; }
G4ParticleChangeForMSC * fParticleChange
static G4Pow * GetInstance()
ThreeVector shoot(const G4int Ap, const G4int Af)
G4IonisParamMat * GetIonisation() const
static c2_factory< G4double > c2
CLHEP::HepRandomEngine * rndmEngineMod
static G4LossTableManager * Instance()
G4UrbanMscModel(const G4String &nam="UrbanMsc")
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const G4DynamicParticle * GetDynamicParticle() const
virtual ~G4UrbanMscModel()
static const G4double Tdat[22]
G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety)
G4bool latDisplasmentbackup
G4double GetProductionCut(G4int index) const
G4StepStatus GetStepStatus() const
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
static const G4double rp2
const G4double w[NPOINTSGL]
static const G4double rp3
static const G4double eps
G4ParticleDefinition * GetDefinition() const
const G4ParticleDefinition * particle
static const G4double reps
G4double currentRadLength
const G4Step * GetStep() const
G4double ComputeSafety(const G4ThreeVector &position, G4double limit=DBL_MAX)
static const G4double Tlim
G4int currentMaterialIndex
G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy)
static const G4double bg2lim
G4double GetZeffective() const
static const G4double rellossmax
static const G4double invmev
const G4ParticleDefinition * positron
G4StepPoint * GetPreStepPoint() const
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
static const G4double rp4
void SampleDisplacementNew(G4double sinTheta, G4double phi)
const G4ThreeVector & GetPosition() const
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
G4double Randomizetlimit()
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
static const double twopi
G4double ComputeGeomPathLength(G4double truePathLength)
static const G4double beta2lim
const G4MaterialCutsCouple * couple
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
G4ThreeVector fDisplacement
G4double ComputeTrueStepLength(G4double geomStepLength)
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
G4double GetRadlen() const
static const G4double rp1
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
void SetParticle(const G4ParticleDefinition *)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4Positron * Positron()
static const G4double theta0max
static const G4double sigmafactor
G4LossTableManager * theManager
static const G4double epsfactor
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetDEDX(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
void SampleDisplacement(G4double sinTheta, G4double phi)
void SetCurrentCouple(const G4MaterialCutsCouple *)
static const G4double c_highland
const G4double x[NPOINTSGL]
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetSafety() const
G4double Z23(G4int Z) const
G4MscStepLimitType steppingAlgorithm
static const G4double third
static const G4double rp0
G4double currentKinEnergy
G4double ComputeTheta0(G4double truePathLength, G4double KineticEnergy)
static const G4double sig0[15]
G4ProductionCuts * GetProductionCuts() const
G4double SimpleScattering(G4double xmeanth, G4double x2meanth)
void StartTracking(G4Track *)
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=0., G4double emax=DBL_MAX)
const G4Material * GetMaterial() const
G4double ComputeTruePathLengthLimit(const G4Track &track, G4double ¤tMinimalStep)