Geant4  10.02.p02
G4UrbanMscModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: $
27 // GEANT4 tag $Name: $
28 //
29 // -------------------------------------------------------------------
30 //
31 // GEANT4 Class file
32 //
33 //
34 // File name: G4UrbanMscModel
35 //
36 // Author: Laszlo Urban
37 //
38 // Creation date: 19.02.2013
39 //
40 // Created from G4UrbanMscModel96
41 //
42 // New parametrization for theta0
43 // Correction for very small step length
44 //
45 // Class Description:
46 //
47 // Implementation of the model of multiple scattering based on
48 // H.W.Lewis Phys Rev 78 (1950) 526 and others
49 
50 // -------------------------------------------------------------------
51 // In its present form the model can be used for simulation
52 // of the e-/e+ multiple scattering
53 //
54 
55 
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
59 #include "G4UrbanMscModel.hh"
60 #include "G4PhysicalConstants.hh"
61 #include "G4SystemOfUnits.hh"
62 #include "Randomize.hh"
63 #include "G4Electron.hh"
64 #include "G4Positron.hh"
65 #include "G4LossTableManager.hh"
66 #include "G4ParticleChangeForMSC.hh"
67 
68 #include "G4Poisson.hh"
69 #include "G4Pow.hh"
70 #include "globals.hh"
71 #include "G4Log.hh"
72 #include "G4Exp.hh"
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 
76 using namespace std;
77 
78 static const G4double Tlim = 10.*CLHEP::MeV;
79 static const G4double sigmafactor =
80  CLHEP::twopi*CLHEP::classic_electr_radius*CLHEP::classic_electr_radius;
81 static const G4double epsfactor = 2.*CLHEP::electron_mass_c2*
82  CLHEP::electron_mass_c2*CLHEP::Bohr_radius*CLHEP::Bohr_radius
83  /(CLHEP::hbarc*CLHEP::hbarc);
84 static const G4double beta2lim = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
85  ((Tlim+CLHEP::electron_mass_c2)*(Tlim+CLHEP::electron_mass_c2));
86 static const G4double bg2lim = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
87  (CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
88 
89 static const G4double sig0[15] = {
90  0.2672*CLHEP::barn, 0.5922*CLHEP::barn, 2.653*CLHEP::barn, 6.235*CLHEP::barn,
91  11.69*CLHEP::barn , 13.24*CLHEP::barn , 16.12*CLHEP::barn, 23.00*CLHEP::barn,
92  35.13*CLHEP::barn , 39.95*CLHEP::barn , 50.85*CLHEP::barn, 67.19*CLHEP::barn,
93  91.15*CLHEP::barn , 104.4*CLHEP::barn , 113.1*CLHEP::barn};
94 
95 static const G4double Tdat[22] = {
96  100*CLHEP::eV, 200*CLHEP::eV, 400*CLHEP::eV, 700*CLHEP::eV,
97  1*CLHEP::keV, 2*CLHEP::keV, 4*CLHEP::keV, 7*CLHEP::keV,
98  10*CLHEP::keV, 20*CLHEP::keV, 40*CLHEP::keV, 70*CLHEP::keV,
99  100*CLHEP::keV, 200*CLHEP::keV, 400*CLHEP::keV, 700*CLHEP::keV,
100  1*CLHEP::MeV, 2*CLHEP::MeV, 4*CLHEP::MeV, 7*CLHEP::MeV,
101  10*CLHEP::MeV, 20*CLHEP::MeV};
102 
103 static const G4double reps = 1.e-6;
104 static const G4double rp0 = 2.2747e+4;
105 static const G4double rp1 = 4.5980e+0;
106 static const G4double rp2 = 1.5580e+1;
107 static const G4double rp3 = 7.1287e-1;
108 static const G4double rp4 =-5.7069e-1;
109 
110 static const G4double theta0max = CLHEP::pi/6.;
111 static const G4double rellossmax= 0.50;
112 static const G4double invmev = 1.0/CLHEP::MeV;
113 static const G4double third = 1./3.;
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116 
117 G4UrbanMscModel::G4UrbanMscModel(const G4String& nam)
118  : G4VMscModel(nam)
119 {
120  masslimite = 0.6*MeV;
121  lambdalimit = 1.*mm;
122  fr = 0.02;
123  taubig = 8.0;
124  tausmall = 1.e-16;
125  taulim = 1.e-6;
126  currentTau = taulim;
127  tlimitminfix = 0.01*nm;
128  tlimitminfix2 = 1.*nm;
129  stepmin = tlimitminfix;
130  smallstep = 1.e10;
131  currentRange = 0. ;
132  rangeinit = 0.;
133  tlimit = 1.e10*mm;
134  tlimitmin = 10.*tlimitminfix;
135  tgeom = 1.e50*mm;
136  geombig = 1.e50*mm;
137  geommin = 1.e-3*mm;
138  geomlimit = geombig;
139  presafety = 0.*mm;
140 
141  facsafety = 0.6;
142 
143  Zold = 0.;
144  Zeff = 1.;
145  Z2 = 1.;
146  Z23 = 1.;
147  lnZ = 0.;
148  coeffth1 = 0.;
149  coeffth2 = 0.;
150  coeffc1 = 0.;
151  coeffc2 = 0.;
152  coeffc3 = 0.;
153  coeffc4 = 0.;
154  particle = 0;
155 
156  positron = G4Positron::Positron();
157  theManager = G4LossTableManager::Instance();
158  rndmEngineMod = G4Random::getTheEngine();
159 
160  firstStep = true;
161  insideskin = false;
162  latDisplasmentbackup = false;
163  displacementFlag = true;
164 
165  rangecut = geombig;
166  drr = 0.35 ;
167  finalr = 10.*um ;
168 
169  skindepth = skin*stepmin;
170 
171  mass = proton_mass_c2;
172  charge = ChargeSquare = 1.0;
173  currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength
174  = zPathLength = par1 = par2 = par3 = 0;
175 
176  currentMaterialIndex = -1;
177  fParticleChange = 0;
178  couple = 0;
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 
183 G4UrbanMscModel::~G4UrbanMscModel()
184 {}
185 
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
187 
188 void G4UrbanMscModel::Initialise(const G4ParticleDefinition* p,
189  const G4DataVector&)
190 {
191  // set values of some data members
192  SetParticle(p);
193  /*
194  if(p->GetPDGMass() > MeV) {
195  G4cout << "### WARNING: G4UrbanMscModel model is used for "
196  << p->GetParticleName() << " !!! " << G4endl;
197  G4cout << "### This model should be used only for e+-"
198  << G4endl;
199  }
200  */
201  fParticleChange = GetParticleChangeForMSC(p);
202 
203  latDisplasmentbackup = latDisplasment;
204 
205  //G4cout << "### G4UrbanMscModel::Initialise done!" << G4endl;
206 }
207 
208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209 
210 G4double G4UrbanMscModel::ComputeCrossSectionPerAtom(
211  const G4ParticleDefinition* part,
212  G4double KineticEnergy,
213  G4double AtomicNumber,G4double,
214  G4double, G4double)
215 {
216  static const G4double epsmin = 1.e-4 , epsmax = 1.e10;
217 
218  static const G4double Zdat[15] = { 4., 6., 13., 20., 26., 29., 32., 38.,47.,
219  50., 56., 64., 74., 79., 82. };
220 
221  // corr. factors for e-/e+ lambda for T <= Tlim
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 }};
268 
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 }};
315 
316  //data/corrections for T > Tlim
317 
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,
321  -22.30};
322 
323  G4double sigma;
324  SetParticle(part);
325 
326  Z23 = G4Pow::GetInstance()->Z23(G4lrint(AtomicNumber));
327 
328  // correction if particle .ne. e-/e+
329  // compute equivalent kinetic energy
330  // lambda depends on p*beta ....
331 
332  G4double eKineticEnergy = KineticEnergy;
333 
334  if(mass > electron_mass_c2)
335  {
336  G4double TAU = KineticEnergy/mass ;
337  G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
338  G4double w = c-2. ;
339  G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
340  eKineticEnergy = electron_mass_c2*tau ;
341  }
342 
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);
348 
349  G4double eps = epsfactor*bg2/Z23;
350 
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;
354 
355  sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
356 
357  // interpolate in AtomicNumber and beta2
358  G4double c1,c2,cc1,cc2,corr;
359 
360  // get bin number in Z
361  G4int iZ = 14;
362  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
363  while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
364  if (iZ==14) iZ = 13;
365  if (iZ==-1) iZ = 0 ;
366 
367  G4double ZZ1 = Zdat[iZ];
368  G4double ZZ2 = Zdat[iZ+1];
369  G4double ratZ = (AtomicNumber-ZZ1)*(AtomicNumber+ZZ1)/
370  ((ZZ2-ZZ1)*(ZZ2+ZZ1));
371 
372  if(eKineticEnergy <= Tlim)
373  {
374  // get bin number in T (beta2)
375  G4int iT = 21;
376  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
377  while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
378  if(iT==21) iT = 20;
379  if(iT==-1) iT = 0 ;
380 
381  // calculate betasquare values
382  G4double T = Tdat[iT], E = T + electron_mass_c2;
383  G4double b2small = T*(E+electron_mass_c2)/(E*E);
384 
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);
388 
389  if (charge < 0.)
390  {
391  c1 = celectron[iZ][iT];
392  c2 = celectron[iZ+1][iT];
393  cc1 = c1+ratZ*(c2-c1);
394 
395  c1 = celectron[iZ][iT+1];
396  c2 = celectron[iZ+1][iT+1];
397  cc2 = c1+ratZ*(c2-c1);
398 
399  corr = cc1+ratb2*(cc2-cc1);
400 
401  sigma *= sigmafactor/corr;
402  }
403  else
404  {
405  c1 = cpositron[iZ][iT];
406  c2 = cpositron[iZ+1][iT];
407  cc1 = c1+ratZ*(c2-c1);
408 
409  c1 = cpositron[iZ][iT+1];
410  c2 = cpositron[iZ+1][iT+1];
411  cc2 = c1+ratZ*(c2-c1);
412 
413  corr = cc1+ratb2*(cc2-cc1);
414 
415  sigma *= sigmafactor/corr;
416  }
417  }
418  else
419  {
420  c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
421  c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
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);
428  }
429  return sigma;
430 
431 }
432 
433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
434 
435 void G4UrbanMscModel::StartTracking(G4Track* track)
436 {
437  SetParticle(track->GetDynamicParticle()->GetDefinition());
438  firstStep = true;
439  insideskin = false;
440  fr = facrange;
441  tlimit = tgeom = rangeinit = rangecut = geombig;
442  smallstep = 1.e10;
443  stepmin = tlimitminfix;
444  tlimitmin = 10.*tlimitminfix;
445  rndmEngineMod = G4Random::getTheEngine();
446 }
447 
448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
449 
450 G4double G4UrbanMscModel::ComputeTruePathLengthLimit(
451  const G4Track& track,
452  G4double& currentMinimalStep)
453 {
454  tPathLength = currentMinimalStep;
455  const G4DynamicParticle* dp = track.GetDynamicParticle();
456 
457  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
458  G4StepStatus stepStatus = sp->GetStepStatus();
459  couple = track.GetMaterialCutsCouple();
460  SetCurrentCouple(couple);
461  currentMaterialIndex = couple->GetIndex();
462  currentKinEnergy = dp->GetKineticEnergy();
463  currentRange = GetRange(particle,currentKinEnergy,couple);
464  lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy);
465  tPathLength = min(tPathLength,currentRange);
466  /*
467  G4cout << "G4Urban::StepLimit tPathLength= " << tPathLength
468  << " range= " <<currentRange<< " lambda= "<<lambda0
469  <<G4endl;
470  */
471  // set flag to default values
472  Zeff = couple->GetMaterial()->GetIonisation()->GetZeffective();
473  // couple->GetMaterial()->GetTotNbOfAtomsPerVolume();
474 
475  if(Zold != Zeff)
476  UpdateCache();
477 
478  // stop here if small step
479  if(tPathLength < tlimitminfix) {
480  latDisplasment = false;
481  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
482  }
483 
484  // upper limit for the straight line distance the particle can travel
485  // for electrons and positrons
486  G4double distance = currentRange;
487  // for muons, hadrons
488  if(mass > masslimite) {
489  distance *= (1.15-9.76e-4*Zeff);
490  } else {
491  distance *= (1.20-Zeff*(1.62e-2-9.22e-5*Zeff));
492  }
493  presafety = sp->GetSafety();
494  /*
495  G4cout << "G4Urban::StepLimit tPathLength= "
496  <<tPathLength<<" safety= " << presafety
497  << " range= " <<currentRange<< " lambda= "<<lambda0
498  << " Alg: " << steppingAlgorithm <<G4endl;
499  */
500  // far from geometry boundary
501  if(distance < presafety)
502  {
503  latDisplasment = false;
504  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
505  }
506 
507  latDisplasment = latDisplasmentbackup;
508  // standard version
509  //
510  if (steppingAlgorithm == fUseDistanceToBoundary)
511  {
512  //compute geomlimit and presafety
513  geomlimit = ComputeGeomLimit(track, presafety, currentRange);
514  /*
515  G4cout << "G4Urban::Distance to boundary geomlimit= "
516  <<geomlimit<<" safety= " << presafety<<G4endl;
517  */
518 
519  // is it far from boundary ?
520  if(distance < presafety)
521  {
522  latDisplasment = false;
523  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
524  }
525 
526  smallstep += 1.;
527  insideskin = false;
528 
529  // initialisation at firs step and at the boundary
530  if(firstStep || (stepStatus == fGeomBoundary))
531  {
532  rangeinit = currentRange;
533  if(!firstStep) { smallstep = 1.; }
534 
535  //define stepmin here (it depends on lambda!)
536  //rough estimation of lambda_elastic/lambda_transport
537  G4double rat = currentKinEnergy*invmev;
538  rat = 1.e-3/(rat*(10.+rat)) ;
539  //stepmin ~ lambda_elastic
540  stepmin = rat*lambda0;
541  skindepth = skin*stepmin;
542  tlimitmin = max(10*stepmin,tlimitminfix);
543  /*
544  G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
545  << " tlimitmin= " << tlimitmin << " geomlimit= "
546  << geomlimit <<G4endl;
547  */
548  // constraint from the geometry
549 
550  if((geomlimit < geombig) && (geomlimit > geommin))
551  {
552  // geomlimit is a geometrical step length
553  // transform it to true path length (estimation)
554  if((1.-geomlimit/lambda0) > 0.)
555  geomlimit = -lambda0*G4Log(1.-geomlimit/lambda0)+tlimitmin ;
556 
557  if(stepStatus == fGeomBoundary)
558  tgeom = geomlimit/facgeom;
559  else
560  tgeom = 2.*geomlimit/facgeom;
561  }
562  else
563  tgeom = geombig;
564  }
565 
566  //step limit
567  tlimit = facrange*rangeinit;
568 
569  //lower limit for tlimit
570  tlimit = max(tlimit,tlimitmin);
571  tlimit = min(tlimit,tgeom);
572  /*
573  G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
574  << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
575  */
576  // shortcut
577  if((tPathLength < tlimit) && (tPathLength < presafety) &&
578  (smallstep > skin) && (tPathLength < geomlimit-0.999*skindepth))
579  {
580  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
581  }
582 
583  // step reduction near to boundary
584  if(smallstep <= skin)
585  {
586  tlimit = stepmin;
587  insideskin = true;
588  }
589  else if(geomlimit < geombig)
590  {
591  if(geomlimit > skindepth)
592  {
593  tlimit = min(tlimit, geomlimit-0.999*skindepth);
594  }
595  else
596  {
597  insideskin = true;
598  tlimit = min(tlimit, stepmin);
599  }
600  }
601 
602  tlimit = max(tlimit, stepmin);
603 
604  // randomize 1st step or 1st 'normal' step in volume
605  if(firstStep || ((smallstep == skin+1) && !insideskin))
606  {
607  tPathLength = min(tPathLength, Randomizetlimit());
608  }
609  else
610  {
611  tPathLength = min(tPathLength, tlimit);
612  }
613 
614  }
615  // for 'normal' simulation with or without magnetic field
616  // there no small step/single scattering at boundaries
617  else if(steppingAlgorithm == fUseSafety)
618  {
619  if(stepStatus != fGeomBoundary) {
620  presafety = ComputeSafety(sp->GetPosition(),tPathLength);
621  }
622  /*
623  G4cout << "presafety= " << presafety
624  << " firstStep= " << firstStep
625  << " stepStatus= " << stepStatus
626  << G4endl;
627  */
628  // is far from boundary
629  if(distance < presafety)
630  {
631  latDisplasment = false;
632  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
633  }
634 
635  if(firstStep || (stepStatus == fGeomBoundary)) {
636  rangeinit = currentRange;
637  fr = facrange;
638  // 9.1 like stepping for e+/e- only (not for muons,hadrons)
639  if(mass < masslimite)
640  {
641  rangeinit = max(rangeinit, lambda0);
642  if(lambda0 > lambdalimit) {
643  fr *= (0.75+0.25*lambda0/lambdalimit);
644  }
645  }
646  //lower limit for tlimit
647  G4double rat = currentKinEnergy*invmev;
648  rat = 1.e-3/(rat*(10 + rat)) ;
649  stepmin = lambda0*rat;
650  tlimitmin = max(10*stepmin, tlimitminfix);
651  }
652 
653  //step limit
654  tlimit = max(fr*rangeinit, facsafety*presafety);
655 
656  //lower limit for tlimit
657  tlimit = max(tlimit, tlimitmin);
658 
659  if(firstStep || stepStatus == fGeomBoundary)
660  {
661  tPathLength = min(tPathLength, Randomizetlimit());
662  }
663  else { tPathLength = min(tPathLength, tlimit); }
664  }
665  // new stepping mode UseSafetyPlus
666  else if(steppingAlgorithm == fUseSafetyPlus)
667  {
668  if(stepStatus != fGeomBoundary) {
669  presafety = ComputeSafety(sp->GetPosition(),tPathLength);
670  }
671  /*
672  G4cout << "presafety= " << presafety
673  << " firstStep= " << firstStep
674  << " stepStatus= " << stepStatus
675  << G4endl;
676  */
677  // is far from boundary
678  if(distance < presafety)
679  {
680  latDisplasment = false;
681  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
682  }
683 
684  if(firstStep || (stepStatus == fGeomBoundary)) {
685  rangeinit = currentRange;
686  fr = facrange;
687  rangecut = geombig;
688  if(mass < masslimite)
689  {
690  G4int index = 1;
691  if(charge > 0.) index = 2;
692  rangecut = couple->GetProductionCuts()->GetProductionCut(index);
693  if(lambda0 > lambdalimit) {
694  fr *= (0.84+0.16*lambda0/lambdalimit);
695  }
696  }
697  //lower limit for tlimit
698  G4double rat = currentKinEnergy*invmev;
699  rat = 1.e-3/(rat*(10 + rat)) ;
700  stepmin = lambda0*rat;
701  tlimitmin = max(10*stepmin, tlimitminfix);
702  }
703  //step limit
704  tlimit = max(fr*rangeinit, facsafety*presafety);
705 
706  //lower limit for tlimit
707  tlimit = max(tlimit, tlimitmin);
708 
709  // condition for tPathLength from drr and finalr
710  if(currentRange > finalr) {
711  G4double tmax = drr*currentRange+
712  finalr*(1.-drr)*(2.-finalr/currentRange);
713  tPathLength = min(tPathLength,tmax);
714  }
715 
716  // condition safety
717  if(currentRange > rangecut) {
718  if(firstStep) {
719  tPathLength = min(tPathLength,facsafety*presafety);
720  } else if(stepStatus != fGeomBoundary && presafety > stepmin) {
721  tPathLength = min(tPathLength,presafety);
722  }
723  }
724 
725  if(firstStep || stepStatus == fGeomBoundary)
726  {
727  tPathLength = min(tPathLength, Randomizetlimit());
728  }
729  else { tPathLength = min(tPathLength, tlimit); }
730  }
731 
732  // version similar to 7.1 (needed for some experiments)
733  else
734  {
735  if (stepStatus == fGeomBoundary)
736  {
737  if (currentRange > lambda0) { tlimit = facrange*currentRange; }
738  else { tlimit = facrange*lambda0; }
739 
740  tlimit = max(tlimit, tlimitmin);
741  }
742  if(firstStep || stepStatus == fGeomBoundary)
743  {
744  tPathLength = min(tPathLength, Randomizetlimit());
745  }
746  else { tPathLength = min(tPathLength, tlimit); }
747  }
748  firstStep = false;
749  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
750 }
751 
752 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
753 
754 G4double G4UrbanMscModel::ComputeGeomPathLength(G4double)
755 {
756  lambdaeff = lambda0;
757  par1 = -1. ;
758  par2 = par3 = 0. ;
759 
760  // this correction needed to run MSC with eIoni and eBrem inactivated
761  // and makes no harm for a normal run
762  tPathLength = std::min(tPathLength,currentRange);
763 
764  // do the true -> geom transformation
765  zPathLength = tPathLength;
766 
767  // z = t for very small tPathLength
768  if(tPathLength < tlimitminfix2) return zPathLength;
769 
770  // VI: it is already checked
771  // if(tPathLength > currentRange)
772  // tPathLength = currentRange ;
773  /*
774  G4cout << "ComputeGeomPathLength: tpl= " << tPathLength
775  << " R= " << currentRange << " L0= " << lambda0
776  << " E= " << currentKinEnergy << " "
777  << particle->GetParticleName() << G4endl;
778  */
779  G4double tau = tPathLength/lambda0 ;
780 
781  if ((tau <= tausmall) || insideskin) {
782  zPathLength = min(tPathLength, lambda0);
783 
784  } else if (tPathLength < currentRange*dtrl) {
785  if(tau < taulim) zPathLength = tPathLength*(1.-0.5*tau) ;
786  else zPathLength = lambda0*(1.-G4Exp(-tau));
787 
788  } else if(currentKinEnergy < mass || tPathLength == currentRange) {
789  par1 = 1./currentRange ;
790  par2 = 1./(par1*lambda0) ;
791  par3 = 1.+par2 ;
792  if(tPathLength < currentRange) {
793  zPathLength =
794  (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3);
795  } else {
796  zPathLength = 1./(par1*par3);
797  }
798 
799  } else {
800  G4double rfin = max(currentRange-tPathLength, 0.01*currentRange);
801  G4double T1 = GetEnergy(particle,rfin,couple);
802  G4double lambda1 = GetTransportMeanFreePath(particle,T1);
803 
804  par1 = (lambda0-lambda1)/(lambda0*tPathLength);
805  //G4cout << "par1= " << par1 << " L1= " << lambda1 << G4endl;
806  par2 = 1./(par1*lambda0);
807  par3 = 1.+par2 ;
808  zPathLength = (1.-G4Exp(par3*G4Log(lambda1/lambda0)))/(par1*par3);
809  }
810 
811  zPathLength = min(zPathLength, lambda0);
812  //G4cout<< "zPathLength= "<< zPathLength<< " L0= " << lambda0 << G4endl;
813  return zPathLength;
814 }
815 
816 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
817 
818 G4double G4UrbanMscModel::ComputeTrueStepLength(G4double geomStepLength)
819 {
820  // step defined other than transportation
821  if(geomStepLength == zPathLength) {
822  //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength
823  // << " step= " << geomStepLength << " *** " << G4endl;
824  return tPathLength;
825  }
826 
827  zPathLength = geomStepLength;
828 
829  // t = z for very small step
830  if(geomStepLength < tlimitminfix2) {
831  tPathLength = geomStepLength;
832 
833  // recalculation
834  } else {
835 
836  G4double tlength = geomStepLength;
837  if((geomStepLength > lambda0*tausmall) && !insideskin) {
838 
839  if(par1 < 0.) {
840  tlength = -lambda0*G4Log(1.-geomStepLength/lambda0) ;
841  } else {
842  if(par1*par3*geomStepLength < 1.) {
843  tlength = (1.-G4Exp(G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
844  } else {
845  tlength = currentRange;
846  }
847  }
848 
849  if(tlength < geomStepLength) { tlength = geomStepLength; }
850  else if(tlength > tPathLength) { tlength = tPathLength; }
851  }
852  tPathLength = tlength;
853  }
854  //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength
855  // << " step= " << geomStepLength << " &&& " << G4endl;
856 
857  return tPathLength;
858 }
859 
860 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
861 
862 G4ThreeVector&
863 G4UrbanMscModel::SampleScattering(const G4ThreeVector& oldDirection,
864  G4double /*safety*/)
865 {
866  fDisplacement.set(0.0,0.0,0.0);
867  G4double kineticEnergy = currentKinEnergy;
868  if (tPathLength > currentRange*dtrl) {
869  kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
870  } else {
871  kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple);
872  }
873 
874  if((kineticEnergy <= eV) || (tPathLength <= tlimitminfix) ||
875  (tPathLength < tausmall*lambda0)) { return fDisplacement; }
876 
877  G4double cth = SampleCosineTheta(tPathLength,kineticEnergy);
878 
879  // protection against 'bad' cth values
880  if(std::fabs(cth) >= 1.0) { return fDisplacement; }
881 
882  /*
883  if(cth < 1.0 - 1000*tPathLength/lambda0 && cth < 0.5 &&
884  kineticEnergy > 20*MeV) {
885  G4cout << "### G4UrbanMscModel::SampleScattering for "
886  << particle->GetParticleName()
887  << " E(MeV)= " << kineticEnergy/MeV
888  << " Step(mm)= " << tPathLength/mm
889  << " in " << CurrentCouple()->GetMaterial()->GetName()
890  << " CosTheta= " << cth << G4endl;
891  }
892  */
893  G4double sth = sqrt((1.0 - cth)*(1.0 + cth));
894  G4double phi = twopi*rndmEngineMod->flat();
895  G4double dirx = sth*cos(phi);
896  G4double diry = sth*sin(phi);
897 
898  G4ThreeVector newDirection(dirx,diry,cth);
899  newDirection.rotateUz(oldDirection);
900 
901  fParticleChange->ProposeMomentumDirection(newDirection);
902  /*
903  G4cout << "G4UrbanMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
904  << " sinTheta= " << sth << " safety(mm)= " << safety
905  << " trueStep(mm)= " << tPathLength
906  << " geomStep(mm)= " << zPathLength
907  << G4endl;
908  */
909 
910 
911  if (latDisplasment && currentTau >= tausmall) {
912  if(displacementFlag) { SampleDisplacementNew(cth, phi); }
913  else { SampleDisplacement(sth, phi); }
914  fDisplacement.rotateUz(oldDirection);
915  }
916  return fDisplacement;
917 }
918 
919 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
920 
921 G4double G4UrbanMscModel::SampleCosineTheta(G4double trueStepLength,
922  G4double KineticEnergy)
923 {
924  G4double cth = 1. ;
925  G4double tau = trueStepLength/lambda0;
926  currentTau = tau;
927  lambdaeff = lambda0;
928 
929  G4double lambda1 = GetTransportMeanFreePath(particle,KineticEnergy);
930  if(std::fabs(lambda1 - lambda0) > lambda0*0.01 && lambda1 > 0.)
931  {
932  // mean tau value
933  tau = trueStepLength*G4Log(lambda0/lambda1)/(lambda0-lambda1);
934  }
935 
936  currentTau = tau ;
937  lambdaeff = trueStepLength/currentTau;
938  currentRadLength = couple->GetMaterial()->GetRadlen();
939 
940  if (tau >= taubig) { cth = -1.+2.*rndmEngineMod->flat(); }
941  else if (tau >= tausmall) {
942  static const G4double numlim = 0.01;
943  G4double xmeanth, x2meanth;
944  if(tau < numlim) {
945  xmeanth = 1.0 - tau*(1.0 - 0.5*tau);
946  x2meanth= 1.0 - tau*(5.0 - 6.25*tau)/3.;
947  } else {
948  xmeanth = G4Exp(-tau);
949  x2meanth = (1.+2.*G4Exp(-2.5*tau))/3.;
950  }
951 
952  // too large step of low-energy particle
953  G4double relloss = 1. - KineticEnergy/currentKinEnergy;
954  if(relloss > rellossmax) {
955  return SimpleScattering(xmeanth,x2meanth);
956  }
957  // is step extreme small ?
958  G4bool extremesmallstep = false ;
959  G4double tsmall = std::min(tlimitmin,lambdalimit);
960  G4double theta0 = 0.;
961  if(trueStepLength > tsmall) {
962  theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
963  } else {
964  theta0 = sqrt(trueStepLength/tsmall)*ComputeTheta0(tsmall,KineticEnergy);
965  extremesmallstep = true ;
966  }
967 
968  //G4cout << "Theta0= " << theta0 << " theta0max= " << theta0max
969  // << " sqrt(tausmall)= " << sqrt(tausmall) << G4endl;
970 
971  // protection for very small angles
972  G4double theta2 = theta0*theta0;
973 
974  if(theta2 < tausmall) { return cth; }
975 
976  if(theta0 > theta0max) {
977  return SimpleScattering(xmeanth,x2meanth);
978  }
979 
980  G4double x = theta2*(1.0 - theta2/12.);
981  if(theta2 > numlim) {
982  G4double sth = 2*sin(0.5*theta0);
983  x = sth*sth;
984  }
985 
986  // parameter for tail
987  G4double ltau= G4Log(tau);
988  G4double u = G4Exp(ltau/6.);
989  if(extremesmallstep) u = G4Exp(G4Log(tsmall/lambda0)/6.);
990  G4double xx = G4Log(lambdaeff/currentRadLength);
991  G4double xsi = coeffc1+u*(coeffc2+coeffc3*u)+coeffc4*xx;
992 
993  // tail should not be too big
994  if(xsi < 1.9) {
995  /*
996  if(KineticEnergy > 20*MeV && xsi < 1.6) {
997  G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= "
998  << KineticEnergy/GeV
999  << " !!** c= " << xsi
1000  << " **!! length(mm)= " << trueStepLength << " Zeff= " << Zeff
1001  << " " << couple->GetMaterial()->GetName()
1002  << " tau= " << tau << G4endl;
1003  }
1004  */
1005  xsi = 1.9;
1006  }
1007 
1008  G4double c = xsi;
1009 
1010  if(fabs(c-3.) < 0.001) { c = 3.001; }
1011  else if(fabs(c-2.) < 0.001) { c = 2.001; }
1012 
1013  G4double c1 = c-1.;
1014 
1015  G4double ea = G4Exp(-xsi);
1016  G4double eaa = 1.-ea ;
1017  G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/eaa;
1018  G4double x0 = 1. - xsi*x;
1019 
1020  // G4cout << " xmean1= " << xmean1 << " xmeanth= " << xmeanth << G4endl;
1021 
1022  if(xmean1 <= 0.999*xmeanth) {
1023  return SimpleScattering(xmeanth,x2meanth);
1024  }
1025  //from continuity of derivatives
1026  G4double b = 1.+(c-xsi)*x;
1027 
1028  G4double b1 = b+1.;
1029  G4double bx = c*x;
1030 
1031  G4double eb1 = G4Exp(G4Log(b1)*c1);
1032  G4double ebx = G4Exp(G4Log(bx)*c1);
1033  G4double d = ebx/eb1;
1034 
1035  G4double xmean2 = (x0 + d - (bx - b1*d)/(c-2.))/(1. - d);
1036 
1037  G4double f1x0 = ea/eaa;
1038  G4double f2x0 = c1/(c*(1. - d));
1039  G4double prob = f2x0/(f1x0+f2x0);
1040 
1041  G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
1042 
1043  // sampling of costheta
1044  //G4cout << "c= " << c << " qprob= " << qprob << " eb1= " << eb1
1045  // << " c1= " << c1 << " b1= " << b1 << " bx= " << bx << " eb1= " << eb1
1046  // << G4endl;
1047  if(rndmEngineMod->flat() < qprob)
1048  {
1049  G4double var = 0;
1050  if(rndmEngineMod->flat() < prob) {
1051  cth = 1.+G4Log(ea+rndmEngineMod->flat()*eaa)*x;
1052  } else {
1053  var = (1.0 - d)*rndmEngineMod->flat();
1054  if(var < numlim*d) {
1055  var /= (d*c1);
1056  cth = -1.0 + var*(1.0 - 0.5*var*c)*(2. + (c - xsi)*x);
1057  } else {
1058  cth = 1. + x*(c - xsi - c*G4Exp(-G4Log(var + d)/c1));
1059  }
1060  }
1061  /*
1062  if(KineticEnergy > 5*GeV && cth < 0.9) {
1063  G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= "
1064  << KineticEnergy/GeV
1065  << " 1-cosT= " << 1 - cth
1066  << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
1067  << " tau= " << tau
1068  << " prob= " << prob << " var= " << var << G4endl;
1069  G4cout << " c= " << c << " qprob= " << qprob << " eb1= " << eb1
1070  << " ebx= " << ebx
1071  << " c1= " << c1 << " b= " << b << " b1= " << b1
1072  << " bx= " << bx << " d= " << d
1073  << " ea= " << ea << " eaa= " << eaa << G4endl;
1074  }
1075  */
1076  }
1077  else {
1078  cth = -1.+2.*rndmEngineMod->flat();
1079  /*
1080  if(KineticEnergy > 5*GeV) {
1081  G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= "
1082  << KineticEnergy/GeV
1083  << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
1084  << " qprob= " << qprob << G4endl;
1085  }
1086  */
1087  }
1088  }
1089  return cth ;
1090 }
1091 
1092 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1093 
1094 G4double G4UrbanMscModel::ComputeTheta0(G4double trueStepLength,
1095  G4double KineticEnergy)
1096 {
1097  // for all particles take the width of the central part
1098  // from a parametrization similar to the Highland formula
1099  // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
1100  G4double invbetacp = std::sqrt((currentKinEnergy+mass)*(KineticEnergy+mass)/
1101  (currentKinEnergy*(currentKinEnergy+2.*mass)*
1102  KineticEnergy*(KineticEnergy+2.*mass)));
1103  G4double y = trueStepLength/currentRadLength;
1104 
1105  if(particle == positron)
1106  {
1107  const G4double xl= 0.6;
1108  const G4double xh= 0.9;
1109  const G4double e = 113.0;
1110  G4double corr;
1111 
1112  G4double tau = std::sqrt(currentKinEnergy*KineticEnergy)/mass;
1113  G4double x = std::sqrt(tau*(tau+2.)/((tau+1.)*(tau+1.)));
1114  G4double a = 0.994-4.08e-3*Zeff;
1115  G4double b = 7.16+(52.6+365./Zeff)/Zeff;
1116  G4double c = 1.000-4.47e-3*Zeff;
1117  G4double d = 1.21e-3*Zeff;
1118  if(x < xl) {
1119  corr = a*(1.-G4Exp(-b*x));
1120  } else if(x > xh) {
1121  corr = c+d*G4Exp(e*(x-1.));
1122  } else {
1123  G4double yl = a*(1.-G4Exp(-b*xl));
1124  G4double yh = c+d*G4Exp(e*(xh-1.));
1125  G4double y0 = (yh-yl)/(xh-xl);
1126  G4double y1 = yl-y0*xl;
1127  corr = y0*x+y1;
1128  }
1129  //==================================================================
1130  y *= corr*(1.+Zeff*(1.84035e-4*Zeff-1.86427e-2)+0.41125);
1131  }
1132 
1133  G4double theta0 = c_highland*std::abs(charge)*std::sqrt(y)*invbetacp;
1134 
1135  // correction factor from e- scattering data
1136  theta0 *= (coeffth1+coeffth2*G4Log(y));
1137  return theta0;
1138 }
1139 
1140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1141 
1142 void G4UrbanMscModel::SampleDisplacement(G4double sth, G4double phi)
1143 {
1144  G4double rmax = sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
1145  G4double r = rmax*G4Exp(G4Log(rndmEngineMod->flat())*third);
1146  /*
1147  G4cout << "G4UrbanMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
1148  << " sinTheta= " << sth << " r(mm)= " << r
1149  << " trueStep(mm)= " << tPathLength
1150  << " geomStep(mm)= " << zPathLength
1151  << G4endl;
1152  */
1153 
1154  if(r > 0.) {
1155  static const G4double kappa = 2.5;
1156  static const G4double kappami1 = 1.5;
1157 
1158  G4double latcorr = 0.;
1159  if((currentTau >= tausmall) && !insideskin) {
1160  if(currentTau < taulim) {
1161  latcorr = lambdaeff*kappa*currentTau*currentTau*
1162  (1.-(kappa+1.)*currentTau*third)*third;
1163 
1164  } else {
1165  G4double etau = 0.;
1166  if(currentTau < taubig) { etau = G4Exp(-currentTau); }
1167  latcorr = -kappa*currentTau;
1168  latcorr = G4Exp(latcorr)/kappami1;
1169  latcorr += 1.-kappa*etau/kappami1 ;
1170  latcorr *= 2.*lambdaeff*third;
1171  }
1172  }
1173  latcorr = std::min(latcorr, r);
1174 
1175  // sample direction of lateral displacement
1176  // compute it from the lateral correlation
1177  G4double Phi = 0.;
1178  if(std::abs(r*sth) < latcorr) {
1179  Phi = twopi*rndmEngineMod->flat();
1180 
1181  } else {
1182  //G4cout << "latcorr= " << latcorr << " r*sth= " << r*sth
1183  // << " ratio= " << latcorr/(r*sth) << G4endl;
1184  G4double psi = std::acos(latcorr/(r*sth));
1185  if(rndmEngineMod->flat() < 0.5) {
1186  Phi = phi+psi;
1187  } else {
1188  Phi = phi-psi;
1189  }
1190  }
1191  fDisplacement.set(r*std::cos(Phi),r*std::sin(Phi),0.0);
1192  }
1193 }
1194 
1195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1196 
1197 void G4UrbanMscModel::SampleDisplacementNew(G4double , G4double phi)
1198 {
1199  //sample displacement r
1200 
1201  G4double rmax = sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
1202  // u = (r/rmax)**2 , v=1-u
1203  // paramerization from ss simulation
1204  // f(u) = p0*exp(p1*log(v)-p2*v)+v*(p3+p4*v)
1205  G4double u ,v , rej;
1206  G4int count = 0;
1207  do {
1208  u = reps+(1.-2.*reps)*rndmEngineMod->flat();
1209  v = 1.-u ;
1210  rej = rp0*G4Exp(rp1*G4Log(v)-rp2*v) + v*(rp3+rp4*v);
1211  }
1212  // Loop checking, 15-Sept-2015, Vladimir Ivanchenko
1213  while (rndmEngineMod->flat() > rej && ++count < 1000);
1214  G4double r = rmax*sqrt(u);
1215 
1216  if(r > 0.)
1217  {
1218  // sample Phi using lateral correlation
1219  // v = Phi-phi = acos(latcorr/(r*sth))
1220  // v has a universal distribution which can be parametrized from ss
1221  // simulation as
1222  // f(v) = 1.49e-2*exp(-v**2/(2*0.320))+2.50e-2*exp(-31.0*log(1.+6.30e-2*v))+
1223  // 1.96e-5*exp(8.42e-1*log(1.+1.45e1*v))
1224  static const G4double probv1 = 0.305533;
1225  static const G4double probv2 = 0.955176;
1226  static const G4double vhigh = 3.15;
1227  static const G4double w2v = 1./G4Exp(30.*G4Log(1. + 6.30e-2*vhigh));
1228  static const G4double w3v = 1./G4Exp(-1.842*G4Log(1. + 1.45e1*vhigh));
1229 
1230  G4double Phi;
1231  G4double random = rndmEngineMod->flat();
1232  if(random < probv1) {
1233  do {
1234  v = G4RandGauss::shoot(rndmEngineMod,0.,0.320);
1235  }
1236  // Loop checking, 15-Sept-2015, Vladimir Ivanchenko
1237  while (std::abs(v) >= vhigh);
1238  Phi = phi + v;
1239 
1240  } else {
1241 
1242  if(random < probv2) {
1243  v = (-1.+1./G4Exp(G4Log(1.-rndmEngineMod->flat()*(1.-w2v))/30.))/6.30e-2;
1244  } else {
1245  v = (-1.+1./G4Exp(G4Log(1.-rndmEngineMod->flat()*(1.-w3v))/-1.842))/1.45e1;
1246  }
1247 
1248  random = rndmEngineMod->flat();
1249  if(random < 0.5) { Phi = phi+v; }
1250  else { Phi = phi-v; }
1251  }
1252  fDisplacement.set(r*std::cos(Phi),r*std::sin(Phi),0.0);
1253  }
1254 }
1255 
1256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4Exp.hh
G4UrbanMscModel::fParticleChange
G4ParticleChangeForMSC * fParticleChange
Definition: G4UrbanMscModel.hh:131
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4VMscModel::facgeom
G4double facgeom
Definition: G4VMscModel.hh:178
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:68
G4Material::GetIonisation
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
c2
static c2_factory< G4double > c2
Definition: G4ScreenedNuclearRecoil.cc:126
G4UrbanMscModel::rndmEngineMod
CLHEP::HepRandomEngine * rndmEngineMod
Definition: G4UrbanMscModel.hh:127
MeV
static const double MeV
Definition: G4SIunits.hh:211
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:112
G4UrbanMscModel::G4UrbanMscModel
G4UrbanMscModel(const G4String &nam="UrbanMsc")
Definition: G4UrbanMscModel.cc:117
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4VMscModel::dtrl
G4double dtrl
Definition: G4VMscModel.hh:181
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4UrbanMscModel::~G4UrbanMscModel
virtual ~G4UrbanMscModel()
Definition: G4UrbanMscModel.cc:183
G4VMscModel::facrange
G4double facrange
Definition: G4VMscModel.hh:177
Tdat
static const G4double Tdat[22]
Definition: G4UrbanMscModel.cc:95
G4UrbanMscModel::SampleScattering
G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety)
Definition: G4UrbanMscModel.cc:863
G4ProductionCuts::GetProductionCut
G4double GetProductionCut(G4int index) const
Definition: G4ProductionCuts.hh:157
G4StepPoint::GetStepStatus
G4StepStatus GetStepStatus() const
G4VMscModel::ConvertTrueToGeom
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:247
G4VMscModel::skin
G4double skin
Definition: G4VMscModel.hh:180
G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
rp2
static const G4double rp2
Definition: G4UrbanMscModel.cc:106
w
const G4double w[NPOINTSGL]
Definition: G4PreCompoundFragment.cc:44
fUseSafety
Definition: G4MscStepLimitType.hh:50
rp3
static const G4double rp3
Definition: G4UrbanMscModel.cc:107
G4VMscModel::latDisplasment
G4bool latDisplasment
Definition: G4VMscModel.hh:190
G4LossTableManager.hh
eps
static const G4double eps
Definition: G4ElectroNuclearCrossSection.cc:99
fUseDistanceToBoundary
Definition: G4MscStepLimitType.hh:52
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
fUseSafetyPlus
Definition: G4MscStepLimitType.hh:51
a
G4double a
Definition: TRTMaterials.hh:39
G4UrbanMscModel::particle
const G4ParticleDefinition * particle
Definition: G4UrbanMscModel.hh:129
reps
static const G4double reps
Definition: G4UrbanMscModel.cc:103
G4UrbanMscModel::currentRadLength
G4double currentRadLength
Definition: G4UrbanMscModel.hh:168
G4Track::GetStep
const G4Step * GetStep() const
G4VMscModel::ComputeSafety
G4double ComputeSafety(const G4ThreeVector &position, G4double limit=DBL_MAX)
Definition: G4VMscModel.hh:239
G4int
int G4int
Definition: G4Types.hh:78
Tlim
static const G4double Tlim
Definition: G4UrbanMscModel.cc:78
G4UrbanMscModel::SampleCosineTheta
G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy)
Definition: G4UrbanMscModel.cc:921
G4StepStatus
G4StepStatus
Definition: G4StepStatus.hh:49
bg2lim
static const G4double bg2lim
Definition: G4UrbanMscModel.cc:86
G4IonisParamMat::GetZeffective
G4double GetZeffective() const
Definition: G4IonisParamMat.hh:127
rellossmax
static const G4double rellossmax
Definition: G4UrbanMscModel.cc:111
invmev
static const G4double invmev
Definition: G4UrbanMscModel.cc:112
G4UrbanMscModel::positron
const G4ParticleDefinition * positron
Definition: G4UrbanMscModel.hh:130
G4Step::GetPreStepPoint
G4StepPoint * GetPreStepPoint() const
G4UrbanMscModel.hh
G4VMscModel::GetEnergy
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:315
rp4
static const G4double rp4
Definition: G4UrbanMscModel.cc:108
G4UrbanMscModel::SampleDisplacementNew
void SampleDisplacementNew(G4double sinTheta, G4double phi)
Definition: G4UrbanMscModel.cc:1197
G4StepPoint::GetPosition
const G4ThreeVector & GetPosition() const
G4VMscModel::GetRange
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:295
G4UrbanMscModel::Randomizetlimit
G4double Randomizetlimit()
Definition: G4UrbanMscModel.hh:211
G4VMscModel::GetParticleChangeForMSC
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:90
G4bool
bool G4bool
Definition: G4Types.hh:79
nm
static const double nm
Definition: G4SIunits.hh:111
G4Poisson.hh
twopi
static const double twopi
Definition: G4SIunits.hh:75
G4UrbanMscModel::ComputeGeomPathLength
G4double ComputeGeomPathLength(G4double truePathLength)
Definition: G4UrbanMscModel.cc:754
G4Positron.hh
beta2lim
static const G4double beta2lim
Definition: G4UrbanMscModel.cc:84
G4UrbanMscModel::couple
const G4MaterialCutsCouple * couple
Definition: G4UrbanMscModel.hh:133
Randomize.hh
c1
static const G4double c1
Definition: G4BetheHeitlerModel.cc:80
G4Log.hh
G4VMscModel::ComputeGeomLimit
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:257
G4VMscModel::fDisplacement
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:186
G4UrbanMscModel::ComputeTrueStepLength
G4double ComputeTrueStepLength(G4double geomStepLength)
Definition: G4UrbanMscModel.cc:818
G4VMscModel::GetTransportMeanFreePath
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:352
G4Material::GetRadlen
G4double GetRadlen() const
Definition: G4Material.hh:220
e1
static const G4double e1
Definition: G4KleinNishinaCompton.cc:68
G4PhysicalConstants.hh
rp1
static const G4double rp1
Definition: G4UrbanMscModel.cc:105
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4UrbanMscModel::SetParticle
void SetParticle(const G4ParticleDefinition *)
Definition: G4UrbanMscModel.hh:199
fGeomBoundary
Definition: G4StepStatus.hh:54
G4ParticleChangeForMSC::ProposeMomentumDirection
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4UrbanMscModel::Initialise
void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4UrbanMscModel.cc:188
eV
static const double eV
Definition: G4SIunits.hh:212
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4Electron.hh
theta0max
static const G4double theta0max
Definition: G4UrbanMscModel.cc:110
sigmafactor
static const G4double sigmafactor
Definition: G4UrbanMscModel.cc:79
G4UrbanMscModel::theManager
G4LossTableManager * theManager
Definition: G4UrbanMscModel.hh:134
pi
static const double pi
Definition: G4SIunits.hh:74
epsfactor
static const G4double epsfactor
Definition: G4UrbanMscModel.cc:81
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4VMscModel::facsafety
G4double facsafety
Definition: G4VMscModel.hh:179
G4VMscModel::GetDEDX
G4double GetDEDX(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:280
G4UrbanMscModel::SampleDisplacement
void SampleDisplacement(G4double sinTheta, G4double phi)
Definition: G4UrbanMscModel.cc:1142
G4VEmModel::SetCurrentCouple
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:445
c_highland
static const G4double c_highland
Definition: G4UrbanMscModel.hh:70
x
const G4double x[NPOINTSGL]
Definition: G4PreCompoundFragment.cc:56
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4StepPoint::GetSafety
G4double GetSafety() const
b1
static const G4double b1
Definition: G4BetheHeitlerModel.cc:76
G4Pow::Z23
G4double Z23(G4int Z) const
Definition: G4Pow.hh:154
um
static const double um
Definition: G4SIunits.hh:112
G4VMscModel::steppingAlgorithm
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:187
third
static const G4double third
Definition: G4UrbanMscModel.cc:113
keV
static const double keV
Definition: G4SIunits.hh:213
barn
static const double barn
Definition: G4SIunits.hh:104
rp0
static const G4double rp0
Definition: G4UrbanMscModel.cc:104
G4UrbanMscModel::currentKinEnergy
G4double currentKinEnergy
Definition: G4UrbanMscModel.hh:165
G4double
double G4double
Definition: G4Types.hh:76
G4UrbanMscModel::ComputeTheta0
G4double ComputeTheta0(G4double truePathLength, G4double KineticEnergy)
Definition: G4UrbanMscModel.cc:1094
G4SystemOfUnits.hh
sig0
static const G4double sig0[15]
Definition: G4UrbanMscModel.cc:89
G4Pow.hh
G4MaterialCutsCouple::GetProductionCuts
G4ProductionCuts * GetProductionCuts() const
Definition: G4MaterialCutsCouple.hh:127
G4UrbanMscModel::SimpleScattering
G4double SimpleScattering(G4double xmeanth, G4double x2meanth)
Definition: G4UrbanMscModel.hh:254
mm
static const double mm
Definition: G4SIunits.hh:114
G4UrbanMscModel::StartTracking
void StartTracking(G4Track *)
Definition: G4UrbanMscModel.cc:435
G4ParticleChangeForMSC.hh
G4UrbanMscModel::ComputeCrossSectionPerAtom
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=0., G4double emax=DBL_MAX)
Definition: G4UrbanMscModel.cc:210
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4UrbanMscModel::ComputeTruePathLengthLimit
G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep)
Definition: G4UrbanMscModel.cc:450