ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UrbanMscModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UrbanMscModel.cc
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 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4UrbanMscModel
33 //
34 // Author: Laszlo Urban
35 //
36 // Creation date: 19.02.2013
37 //
38 // Created from G4UrbanMscModel96
39 //
40 // New parametrization for theta0
41 // Correction for very small step length
42 //
43 // Class Description:
44 //
45 // Implementation of the model of multiple scattering based on
46 // H.W.Lewis Phys Rev 78 (1950) 526 and others
47 
48 // -------------------------------------------------------------------
49 // In its present form the model can be used for simulation
50 // of the e-/e+ multiple scattering
51 //
52 
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 #include "G4UrbanMscModel.hh"
58 #include "G4PhysicalConstants.hh"
59 #include "G4SystemOfUnits.hh"
60 #include "Randomize.hh"
61 #include "G4Electron.hh"
62 #include "G4Positron.hh"
63 #include "G4LossTableManager.hh"
64 #include "G4EmParameters.hh"
65 #include "G4ParticleChangeForMSC.hh"
66 
67 #include "G4Poisson.hh"
68 #include "G4Pow.hh"
69 #include "globals.hh"
70 #include "G4Log.hh"
71 #include "G4Exp.hh"
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 
75 using namespace std;
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 
79 G4UrbanMscModel::G4UrbanMscModel(const G4String& nam)
80  : G4VMscModel(nam)
81 {
82  masslimite = 0.6*MeV;
83  fr = 0.02;
84  taubig = 8.0;
85  tausmall = 1.e-16;
86  taulim = 1.e-6;
87  currentTau = taulim;
88  tlimitminfix = 0.01*nm;
89  tlimitminfix2 = 1.*nm;
90  stepmin = tlimitminfix;
91  stepmina = 1.;
92  stepminb = 1.;
93  smallstep = 1.e10;
94  currentRange = 0. ;
95  rangeinit = 0.;
96  tlimit = 1.e10*mm;
97  tlimitmin = 10.*tlimitminfix;
98  tgeom = 1.e50*mm;
99  geombig = 1.e50*mm;
100  geommin = 1.e-3*mm;
101  geomlimit = geombig;
102  presafety = 0.*mm;
103 
104  Zold = 0.;
105  Zeff = 1.;
106  Z2 = 1.;
107  Z23 = 1.;
108  lnZ = 0.;
109  coeffth1 = 0.;
110  coeffth2 = 0.;
111  coeffc1 = 0.;
112  coeffc2 = 0.;
113  coeffc3 = 0.;
114  coeffc4 = 0.;
115  particle = 0;
116 
117  positron = G4Positron::Positron();
118  theManager = G4LossTableManager::Instance();
119  rndmEngineMod = G4Random::getTheEngine();
120 
121  firstStep = true;
122  insideskin = false;
123  latDisplasmentbackup = false;
124  dispAlg96 = true;
125 
126  rangecut = geombig;
127  drr = 0.35;
128  finalr = 10.*um;
129 
130  skindepth = skin*stepmin;
131 
132  mass = proton_mass_c2;
133  charge = ChargeSquare = 1.0;
134  currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength
135  = zPathLength = par1 = par2 = par3 = 0;
136  currentLogKinEnergy = LOG_EKIN_MIN;
137 
138  currentMaterialIndex = -1;
139  fParticleChange = nullptr;
140  couple = nullptr;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
145 G4UrbanMscModel::~G4UrbanMscModel()
146 {}
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149 
150 void G4UrbanMscModel::Initialise(const G4ParticleDefinition* p,
151  const G4DataVector&)
152 {
153  // set values of some data members
154  SetParticle(p);
155  fParticleChange = GetParticleChangeForMSC(p);
156  InitialiseParameters(p);
157 
158  latDisplasmentbackup = latDisplasment;
159  dispAlg96 = (G4EmParameters::Instance()->LateralDisplacementAlg96());
160  /*
161  G4cout << "### G4UrbanMscModel::Initialise done for "
162  << p->GetParticleName() << " type= " << steppingAlgorithm << G4endl;
163  G4cout << " RangeFact= " << facrange << " GeomFact= " << facgeom
164  << " SafetyFact= " << facsafety << " LambdaLim= " << lambdalimit
165  << G4endl;
166  */
167 }
168 
169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
170 
171 G4double G4UrbanMscModel::ComputeCrossSectionPerAtom(
172  const G4ParticleDefinition* part,
173  G4double KineticEnergy,
174  G4double AtomicNumber,G4double,
175  G4double, G4double)
176 {
177  static const G4double epsmin = 1.e-4 , epsmax = 1.e10;
178 
179  static const G4double Zdat[15] = { 4., 6., 13., 20., 26., 29., 32., 38.,47.,
180  50., 56., 64., 74., 79., 82. };
181 
182  // corr. factors for e-/e+ lambda for T <= Tlim
183  static const G4double celectron[15][22] =
184  {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
185  1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
186  1.112,1.108,1.100,1.093,1.089,1.087 },
187  {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
188  1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
189  1.109,1.105,1.097,1.090,1.086,1.082 },
190  {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
191  1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
192  1.131,1.124,1.113,1.104,1.099,1.098 },
193  {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
194  1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
195  1.112,1.105,1.096,1.089,1.085,1.098 },
196  {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
197  1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
198  1.073,1.070,1.064,1.059,1.056,1.056 },
199  {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
200  1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
201  1.074,1.070,1.063,1.059,1.056,1.052 },
202  {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
203  1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
204  1.068,1.064,1.059,1.054,1.051,1.050 },
205  {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
206  1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
207  1.039,1.037,1.034,1.031,1.030,1.036 },
208  {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
209  1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
210  1.031,1.028,1.024,1.022,1.021,1.024 },
211  {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
212  1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
213  1.020,1.017,1.015,1.013,1.013,1.020 },
214  {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
215  1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
216  0.995,0.993,0.993,0.993,0.993,1.011 },
217  {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
218  1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
219  0.974,0.972,0.973,0.974,0.975,0.987 },
220  {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
221  1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
222  0.950,0.947,0.949,0.952,0.954,0.963 },
223  {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
224  1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
225  0.941,0.938,0.940,0.944,0.946,0.954 },
226  {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
227  1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
228  0.933,0.930,0.933,0.936,0.939,0.949 }};
229 
230  static const G4double cpositron[15][22] = {
231  {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
232  1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
233  1.131,1.126,1.117,1.108,1.103,1.100 },
234  {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
235  1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
236  1.138,1.132,1.122,1.113,1.108,1.102 },
237  {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
238  1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
239  1.203,1.190,1.173,1.159,1.151,1.145 },
240  {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
241  1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
242  1.225,1.210,1.191,1.175,1.166,1.174 },
243  {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
244  1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
245  1.217,1.203,1.184,1.169,1.160,1.151 },
246  {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
247  1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
248  1.237,1.222,1.201,1.184,1.174,1.159 },
249  {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
250  1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
251  1.252,1.234,1.212,1.194,1.183,1.170 },
252  {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
253  2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
254  1.254,1.237,1.214,1.195,1.185,1.179 },
255  {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
256  2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
257  1.312,1.288,1.258,1.235,1.221,1.205 },
258  {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
259  2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
260  1.320,1.294,1.264,1.240,1.226,1.214 },
261  {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
262  2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
263  1.328,1.302,1.270,1.245,1.231,1.233 },
264  {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
265  2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
266  1.361,1.330,1.294,1.267,1.251,1.239 },
267  {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
268  3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
269  1.409,1.372,1.330,1.298,1.280,1.258 },
270  {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
271  3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
272  1.442,1.400,1.354,1.319,1.299,1.272 },
273  {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
274  3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
275  1.456,1.412,1.364,1.328,1.307,1.282 }};
276 
277  //data/corrections for T > Tlim
278 
279  static const G4double hecorr[15] = {
280  120.70, 117.50, 105.00, 92.92, 79.23, 74.510, 68.29,
281  57.39, 41.97, 36.14, 24.53, 10.21, -7.855, -16.84,
282  -22.30};
283 
284  G4double sigma;
285  SetParticle(part);
286 
287  Z23 = G4Pow::GetInstance()->Z23(G4lrint(AtomicNumber));
288 
289  // correction if particle .ne. e-/e+
290  // compute equivalent kinetic energy
291  // lambda depends on p*beta ....
292 
293  G4double eKineticEnergy = KineticEnergy;
294 
295  if(mass > electron_mass_c2)
296  {
297  G4double TAU = KineticEnergy/mass ;
298  G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
299  G4double w = c-2. ;
300  G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
301  eKineticEnergy = electron_mass_c2*tau ;
302  }
303 
304  G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
305  G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
306  /(eTotalEnergy*eTotalEnergy);
307  G4double bg2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
308  /(electron_mass_c2*electron_mass_c2);
309 
310  static const G4double epsfactor = 2.*CLHEP::electron_mass_c2*
311  CLHEP::electron_mass_c2*CLHEP::Bohr_radius*CLHEP::Bohr_radius
312  /(CLHEP::hbarc*CLHEP::hbarc);
313  G4double eps = epsfactor*bg2/Z23;
314 
315  if (eps<epsmin) sigma = 2.*eps*eps;
316  else if(eps<epsmax) sigma = G4Log(1.+2.*eps)-2.*eps/(1.+2.*eps);
317  else sigma = G4Log(2.*eps)-1.+1./eps;
318 
319  sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
320 
321  // interpolate in AtomicNumber and beta2
322  G4double c1,c2,cc1,cc2,corr;
323 
324  // get bin number in Z
325  G4int iZ = 14;
326  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
327  while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
328  if (iZ==14) iZ = 13;
329  if (iZ==-1) iZ = 0 ;
330 
331  G4double ZZ1 = Zdat[iZ];
332  G4double ZZ2 = Zdat[iZ+1];
333  G4double ratZ = (AtomicNumber-ZZ1)*(AtomicNumber+ZZ1)/
334  ((ZZ2-ZZ1)*(ZZ2+ZZ1));
335 
336  static const G4double Tlim = 10.*CLHEP::MeV;
337  static const G4double sigmafactor =
338  CLHEP::twopi*CLHEP::classic_electr_radius*CLHEP::classic_electr_radius;
339  static const G4double beta2lim = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
340  ((Tlim+CLHEP::electron_mass_c2)*(Tlim+CLHEP::electron_mass_c2));
341  static const G4double bg2lim = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
342  (CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
343 
344  static const G4double sig0[15] = {
345  0.2672*CLHEP::barn, 0.5922*CLHEP::barn, 2.653*CLHEP::barn, 6.235*CLHEP::barn,
346  11.69*CLHEP::barn , 13.24*CLHEP::barn , 16.12*CLHEP::barn, 23.00*CLHEP::barn,
347  35.13*CLHEP::barn , 39.95*CLHEP::barn , 50.85*CLHEP::barn, 67.19*CLHEP::barn,
348  91.15*CLHEP::barn , 104.4*CLHEP::barn , 113.1*CLHEP::barn};
349 
350  static const G4double Tdat[22] = {
351  100*CLHEP::eV, 200*CLHEP::eV, 400*CLHEP::eV, 700*CLHEP::eV,
352  1*CLHEP::keV, 2*CLHEP::keV, 4*CLHEP::keV, 7*CLHEP::keV,
353  10*CLHEP::keV, 20*CLHEP::keV, 40*CLHEP::keV, 70*CLHEP::keV,
354  100*CLHEP::keV, 200*CLHEP::keV, 400*CLHEP::keV, 700*CLHEP::keV,
355  1*CLHEP::MeV, 2*CLHEP::MeV, 4*CLHEP::MeV, 7*CLHEP::MeV,
356  10*CLHEP::MeV, 20*CLHEP::MeV};
357 
358  if(eKineticEnergy <= Tlim)
359  {
360  // get bin number in T (beta2)
361  G4int iT = 21;
362  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
363  while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
364  if(iT==21) iT = 20;
365  if(iT==-1) iT = 0 ;
366 
367  // calculate betasquare values
368  G4double T = Tdat[iT], E = T + electron_mass_c2;
369  G4double b2small = T*(E+electron_mass_c2)/(E*E);
370 
371  T = Tdat[iT+1]; E = T + electron_mass_c2;
372  G4double b2big = T*(E+electron_mass_c2)/(E*E);
373  G4double ratb2 = (beta2-b2small)/(b2big-b2small);
374 
375  if (charge < 0.)
376  {
377  c1 = celectron[iZ][iT];
378  c2 = celectron[iZ+1][iT];
379  cc1 = c1+ratZ*(c2-c1);
380 
381  c1 = celectron[iZ][iT+1];
382  c2 = celectron[iZ+1][iT+1];
383  cc2 = c1+ratZ*(c2-c1);
384 
385  corr = cc1+ratb2*(cc2-cc1);
386 
387  sigma *= sigmafactor/corr;
388  }
389  else
390  {
391  c1 = cpositron[iZ][iT];
392  c2 = cpositron[iZ+1][iT];
393  cc1 = c1+ratZ*(c2-c1);
394 
395  c1 = cpositron[iZ][iT+1];
396  c2 = cpositron[iZ+1][iT+1];
397  cc2 = c1+ratZ*(c2-c1);
398 
399  corr = cc1+ratb2*(cc2-cc1);
400 
401  sigma *= sigmafactor/corr;
402  }
403  }
404  else
405  {
406  c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
407  c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
408  if((AtomicNumber >= ZZ1) && (AtomicNumber <= ZZ2))
409  sigma = c1+ratZ*(c2-c1) ;
410  else if(AtomicNumber < ZZ1)
411  sigma = AtomicNumber*AtomicNumber*c1/(ZZ1*ZZ1);
412  else if(AtomicNumber > ZZ2)
413  sigma = AtomicNumber*AtomicNumber*c2/(ZZ2*ZZ2);
414  }
415  // low energy correction based on theory
416  sigma *= 1.+0.30/(1.+sqrt(1000.*eKineticEnergy));
417 
418  return sigma;
419 
420 }
421 
422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
423 
424 void G4UrbanMscModel::StartTracking(G4Track* track)
425 {
426  SetParticle(track->GetDynamicParticle()->GetDefinition());
427  firstStep = true;
428  insideskin = false;
429  fr = facrange;
430  tlimit = tgeom = rangeinit = rangecut = geombig;
431  smallstep = 1.e10;
432  stepmin = tlimitminfix;
433  tlimitmin = 10.*tlimitminfix;
434  rndmEngineMod = G4Random::getTheEngine();
435 }
436 
437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
438 
439 G4double G4UrbanMscModel::ComputeTruePathLengthLimit(
440  const G4Track& track,
441  G4double& currentMinimalStep)
442 {
443  tPathLength = currentMinimalStep;
444  const G4DynamicParticle* dp = track.GetDynamicParticle();
445 
446  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
447  G4StepStatus stepStatus = sp->GetStepStatus();
448  couple = track.GetMaterialCutsCouple();
449  SetCurrentCouple(couple);
450  currentMaterialIndex = couple->GetIndex();
451  currentKinEnergy = dp->GetKineticEnergy();
452  currentLogKinEnergy = dp->GetLogKineticEnergy();
453  currentRange = GetRange(particle,currentKinEnergy,couple,currentLogKinEnergy);
454  lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy,
455  currentLogKinEnergy);
456  tPathLength = min(tPathLength,currentRange);
457  /*
458  G4cout << "G4Urban::StepLimit tPathLength= " << tPathLength
459  << " range= " <<currentRange<< " lambda= "<<lambda0
460  <<G4endl;
461  */
462  // set flag to default values
463  Zeff = couple->GetMaterial()->GetIonisation()->GetZeffective();
464 
465  if(Zold != Zeff)
466  UpdateCache();
467 
468  // stop here if small step
469  if(tPathLength < tlimitminfix) {
470  latDisplasment = false;
471  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
472  }
473 
474  // upper limit for the straight line distance the particle can travel
475  // for electrons and positrons
476  G4double distance = (mass < masslimite)
477  ? currentRange*(1.20-Zeff*(1.62e-2-9.22e-5*Zeff))
478  // for muons, hadrons
479  : currentRange*(1.15-9.76e-4*Zeff);
480 
481  presafety = sp->GetSafety();
482  /*
483  G4cout << "G4Urban::StepLimit tPathLength= "
484  <<tPathLength<<" safety= " << presafety
485  << " range= " <<currentRange<< " lambda= "<<lambda0
486  << " Alg: " << steppingAlgorithm <<G4endl;
487  */
488  // far from geometry boundary
489  if(distance < presafety)
490  {
491  latDisplasment = false;
492  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
493  }
494 
495  latDisplasment = latDisplasmentbackup;
496  static const G4double invmev = 1.0/CLHEP::MeV;
497  // standard version
498  //
499  if (steppingAlgorithm == fUseDistanceToBoundary)
500  {
501  //compute geomlimit and presafety
502  geomlimit = ComputeGeomLimit(track, presafety, currentRange);
503  /*
504  G4cout << "G4Urban::Distance to boundary geomlimit= "
505  <<geomlimit<<" safety= " << presafety<<G4endl;
506  */
507 
508  // is it far from boundary ?
509  if(distance < presafety)
510  {
511  latDisplasment = false;
512  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
513  }
514 
515  smallstep += 1.;
516  insideskin = false;
517 
518  // initialisation at firs step and at the boundary
519  if(firstStep || (stepStatus == fGeomBoundary))
520  {
521  rangeinit = currentRange;
522  if(!firstStep) { smallstep = 1.; }
523 
524  //define stepmin here (it depends on lambda!)
525  //rough estimation of lambda_elastic/lambda_transport
526  G4double rat = currentKinEnergy*invmev;
527  rat = 1.e-3/(2.e-3+rat*(stepmina+stepminb*rat));
528  //stepmin ~ lambda_elastic
529  stepmin = rat*lambda0;
530  skindepth = skin*stepmin;
531  tlimitmin = max(0.7*sqrt(Zeff)*stepmin,tlimitminfix);
532  /*
533  G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
534  << " tlimitmin= " << tlimitmin << " geomlimit= "
535  << geomlimit <<G4endl;
536  */
537  // constraint from the geometry
538 
539  if((geomlimit < geombig) && (geomlimit > geommin))
540  {
541  // geomlimit is a geometrical step length
542  // transform it to true path length (estimation)
543  if(lambda0 > geomlimit) {
544  geomlimit = -lambda0*G4Log(1.-geomlimit/lambda0)+tlimitmin;
545  }
546  tgeom = (stepStatus == fGeomBoundary)
547  ? geomlimit/facgeom : 2.*geomlimit/facgeom;
548  }
549  else
550  {
551  tgeom = geombig;
552  }
553  }
554 
555  //step limit
556  tlimit = facrange*rangeinit;
557 
558  //lower limit for tlimit
559  tlimit = min(max(tlimit,tlimitmin), tgeom);
560  /*
561  G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
562  << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
563  */
564  // shortcut
565  if((tPathLength < tlimit) && (tPathLength < presafety) &&
566  (smallstep > skin) && (tPathLength < geomlimit-0.999*skindepth))
567  {
568  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
569  }
570 
571  // step reduction near to boundary
572  if(smallstep <= skin)
573  {
574  tlimit = stepmin;
575  insideskin = true;
576  }
577  else if(geomlimit < geombig)
578  {
579  if(geomlimit > skindepth)
580  {
581  tlimit = min(tlimit, geomlimit-0.999*skindepth);
582  }
583  else
584  {
585  insideskin = true;
586  tlimit = min(tlimit, stepmin);
587  }
588  }
589 
590  tlimit = max(tlimit, stepmin);
591 
592  // randomise if not 'small' step and step determined by msc
593  tPathLength = ((tlimit < tPathLength)&&(smallstep > skin)&& !insideskin)
594  ? std::min(tPathLength, Randomizetlimit())
595  : std::min(tPathLength, tlimit);
596  }
597  // for 'normal' simulation with or without magnetic field
598  // there no small step/single scattering at boundaries
599  else if(steppingAlgorithm == fUseSafety)
600  {
601  if(stepStatus != fGeomBoundary) {
602  presafety = ComputeSafety(sp->GetPosition(),tPathLength);
603  }
604  /*
605  G4cout << "presafety= " << presafety
606  << " firstStep= " << firstStep
607  << " stepStatus= " << stepStatus
608  << G4endl;
609  */
610  // is far from boundary
611  if(distance < presafety)
612  {
613  latDisplasment = false;
614  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
615  }
616 
617  if(firstStep || (stepStatus == fGeomBoundary)) {
618  rangeinit = currentRange;
619  fr = facrange;
620  // 9.1 like stepping for e+/e- only (not for muons,hadrons)
621  if(mass < masslimite)
622  {
623  rangeinit = max(rangeinit, lambda0);
624  if(lambda0 > lambdalimit) {
625  fr *= (0.75+0.25*lambda0/lambdalimit);
626  }
627  }
628  //lower limit for tlimit
629  G4double rat = currentKinEnergy*invmev;
630  stepmin = lambda0*1.e-3/(2.e-3+rat*(stepmina+stepminb*rat));
631  tlimitmin = max(0.7*sqrt(Zeff)*stepmin,tlimitminfix);
632  }
633 
634  //step limit
635  tlimit = (currentRange > presafety) ?
636  max(fr*rangeinit, facsafety*presafety) : currentRange;
637 
638  //lower limit for tlimit
639  tlimit = max(tlimit, tlimitmin);
640 
641  // randomise if step determined by msc
642  tPathLength = (tlimit < tPathLength) ?
643  std::min(tPathLength, Randomizetlimit()) : tPathLength;
644  }
645  // new stepping mode UseSafetyPlus
646  else if(steppingAlgorithm == fUseSafetyPlus)
647  {
648  if(stepStatus != fGeomBoundary) {
649  presafety = ComputeSafety(sp->GetPosition(),tPathLength);
650  }
651  /*
652  G4cout << "presafety= " << presafety
653  << " firstStep= " << firstStep
654  << " stepStatus= " << stepStatus
655  << G4endl;
656  */
657  // is far from boundary
658  if(distance < presafety)
659  {
660  latDisplasment = false;
661  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
662  }
663 
664  if(firstStep || (stepStatus == fGeomBoundary)) {
665  rangeinit = currentRange;
666  fr = facrange;
667  rangecut = geombig;
668  if(mass < masslimite)
669  {
670  G4int index = 1;
671  if(charge > 0.) index = 2;
672  rangecut = couple->GetProductionCuts()->GetProductionCut(index);
673  if(lambda0 > lambdalimit) {
674  fr *= (0.84+0.16*lambda0/lambdalimit);
675  }
676  }
677  //lower limit for tlimit
678  G4double rat = currentKinEnergy*invmev;
679  stepmin = lambda0*1.e-3/(2.e-3+rat*(stepmina+stepminb*rat));
680  tlimitmin = max(0.7*sqrt(Zeff)*stepmin,tlimitminfix);
681  }
682  //step limit
683  tlimit = (currentRange > presafety) ?
684  max(fr*rangeinit, facsafety*presafety) : currentRange;
685 
686  //lower limit for tlimit
687  tlimit = max(tlimit, tlimitmin);
688 
689  // condition for tPathLength from drr and finalr
690  if(currentRange > finalr) {
691  G4double tmax = drr*currentRange+
692  finalr*(1.-drr)*(2.-finalr/currentRange);
693  tPathLength = min(tPathLength,tmax);
694  }
695 
696  // condition safety
697  if(currentRange > rangecut) {
698  if(firstStep) {
699  tPathLength = min(tPathLength,facsafety*presafety);
700  } else if(stepStatus != fGeomBoundary && presafety > stepmin) {
701  tPathLength = min(tPathLength,presafety);
702  }
703  }
704 
705  // randomise if step determined by msc
706  tPathLength = (tlimit < tPathLength) ?
707  std::min(tPathLength, Randomizetlimit()) : tPathLength;
708  }
709 
710  // version similar to 7.1 (needed for some experiments)
711  else
712  {
713  if (stepStatus == fGeomBoundary)
714  {
715  tlimit = (currentRange > lambda0)
716  ? facrange*currentRange : facrange*lambda0;
717  tlimit = max(tlimit, tlimitmin);
718  }
719  // randomise if step determined by msc
720  tPathLength = (tlimit < tPathLength) ?
721  std::min(tPathLength, Randomizetlimit()) : tPathLength;
722  }
723  firstStep = false;
724  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
725 }
726 
727 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
728 
729 G4double G4UrbanMscModel::ComputeGeomPathLength(G4double)
730 {
731  lambdaeff = lambda0;
732  par1 = -1. ;
733  par2 = par3 = 0. ;
734 
735  // this correction needed to run MSC with eIoni and eBrem inactivated
736  // and makes no harm for a normal run
737  tPathLength = std::min(tPathLength,currentRange);
738 
739  // do the true -> geom transformation
740  zPathLength = tPathLength;
741 
742  // z = t for very small tPathLength
743  if(tPathLength < tlimitminfix2) return zPathLength;
744 
745  // VI: it is already checked
746  // if(tPathLength > currentRange)
747  // tPathLength = currentRange ;
748  /*
749  G4cout << "ComputeGeomPathLength: tpl= " << tPathLength
750  << " R= " << currentRange << " L0= " << lambda0
751  << " E= " << currentKinEnergy << " "
752  << particle->GetParticleName() << G4endl;
753  */
754  G4double tau = tPathLength/lambda0 ;
755 
756  if ((tau <= tausmall) || insideskin) {
757  zPathLength = min(tPathLength, lambda0);
758 
759  } else if (tPathLength < currentRange*dtrl) {
760  if(tau < taulim) zPathLength = tPathLength*(1.-0.5*tau) ;
761  else zPathLength = lambda0*(1.-G4Exp(-tau));
762 
763  } else if(currentKinEnergy < mass || tPathLength == currentRange) {
764  par1 = 1./currentRange ;
765  par2 = 1./(par1*lambda0) ;
766  par3 = 1.+par2 ;
767  if(tPathLength < currentRange) {
768  zPathLength =
769  (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3);
770  } else {
771  zPathLength = 1./(par1*par3);
772  }
773 
774  } else {
775  G4double rfin = max(currentRange-tPathLength, 0.01*currentRange);
776  G4double T1 = GetEnergy(particle,rfin,couple);
777  G4double lambda1 = GetTransportMeanFreePath(particle,T1);
778 
779  par1 = (lambda0-lambda1)/(lambda0*tPathLength);
780  //G4cout << "par1= " << par1 << " L1= " << lambda1 << G4endl;
781  par2 = 1./(par1*lambda0);
782  par3 = 1.+par2 ;
783  zPathLength = (1.-G4Exp(par3*G4Log(lambda1/lambda0)))/(par1*par3);
784  }
785 
786  zPathLength = min(zPathLength, lambda0);
787  //G4cout<< "zPathLength= "<< zPathLength<< " L0= " << lambda0 << G4endl;
788  return zPathLength;
789 }
790 
791 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
792 
793 G4double G4UrbanMscModel::ComputeTrueStepLength(G4double geomStepLength)
794 {
795  // step defined other than transportation
796  if(geomStepLength == zPathLength) {
797  //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength
798  // << " step= " << geomStepLength << " *** " << G4endl;
799  return tPathLength;
800  }
801 
802  zPathLength = geomStepLength;
803 
804  // t = z for very small step
805  if(geomStepLength < tlimitminfix2) {
806  tPathLength = geomStepLength;
807 
808  // recalculation
809  } else {
810 
811  G4double tlength = geomStepLength;
812  if((geomStepLength > lambda0*tausmall) && !insideskin) {
813 
814  if(par1 < 0.) {
815  tlength = -lambda0*G4Log(1.-geomStepLength/lambda0) ;
816  } else {
817  if(par1*par3*geomStepLength < 1.) {
818  tlength = (1.-G4Exp(G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
819  } else {
820  tlength = currentRange;
821  }
822  }
823 
824  if(tlength < geomStepLength) { tlength = geomStepLength; }
825  else if(tlength > tPathLength) { tlength = tPathLength; }
826  }
827  tPathLength = tlength;
828  }
829  //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength
830  // << " step= " << geomStepLength << " &&& " << G4endl;
831 
832  return tPathLength;
833 }
834 
835 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
836 
837 G4ThreeVector&
838 G4UrbanMscModel::SampleScattering(const G4ThreeVector& oldDirection,
839  G4double /*safety*/)
840 {
841  fDisplacement.set(0.0,0.0,0.0);
842  G4double kineticEnergy = currentKinEnergy;
843  if (tPathLength > currentRange*dtrl) {
844  kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
845  } else {
846  kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple,
847  currentLogKinEnergy);
848  }
849 
850  if((kineticEnergy <= eV) || (tPathLength <= tlimitminfix) ||
851  (tPathLength < tausmall*lambda0)) { return fDisplacement; }
852 
853  G4double cth = SampleCosineTheta(tPathLength,kineticEnergy);
854 
855  // protection against 'bad' cth values
856  if(std::abs(cth) >= 1.0) { return fDisplacement; }
857 
858  /*
859  if(cth < 1.0 - 1000*tPathLength/lambda0 && cth < 0.5 &&
860  kineticEnergy > 20*MeV) {
861  G4cout << "### G4UrbanMscModel::SampleScattering for "
862  << particle->GetParticleName()
863  << " E(MeV)= " << kineticEnergy/MeV
864  << " Step(mm)= " << tPathLength/mm
865  << " in " << CurrentCouple()->GetMaterial()->GetName()
866  << " CosTheta= " << cth << G4endl;
867  }
868  */
869  G4double sth = sqrt((1.0 - cth)*(1.0 + cth));
870  G4double phi = twopi*rndmEngineMod->flat();
871  G4double dirx = sth*cos(phi);
872  G4double diry = sth*sin(phi);
873 
874  G4ThreeVector newDirection(dirx,diry,cth);
875  newDirection.rotateUz(oldDirection);
876 
877  fParticleChange->ProposeMomentumDirection(newDirection);
878  /*
879  G4cout << "G4UrbanMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
880  << " sinTheta= " << sth << " safety(mm)= " << safety
881  << " trueStep(mm)= " << tPathLength
882  << " geomStep(mm)= " << zPathLength
883  << G4endl;
884  */
885 
886  if (latDisplasment && currentTau >= tausmall) {
887  if(dispAlg96) { SampleDisplacement(sth, phi); }
888  else { SampleDisplacementNew(cth, phi); }
889  fDisplacement.rotateUz(oldDirection);
890  }
891  return fDisplacement;
892 }
893 
894 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
895 
896 G4double G4UrbanMscModel::SampleCosineTheta(G4double trueStepLength,
897  G4double KineticEnergy)
898 {
899  G4double cth = 1. ;
900  G4double tau = trueStepLength/lambda0;
901  currentTau = tau;
902  lambdaeff = lambda0;
903 
904  G4double lambda1 = GetTransportMeanFreePath(particle,KineticEnergy);
905  if(std::abs(lambda1 - lambda0) > lambda0*0.01 && lambda1 > 0.)
906  {
907  // mean tau value
908  tau = trueStepLength*G4Log(lambda0/lambda1)/(lambda0-lambda1);
909  }
910 
911  currentTau = tau ;
912  lambdaeff = trueStepLength/currentTau;
913  currentRadLength = couple->GetMaterial()->GetRadlen();
914 
915  if (tau >= taubig) { cth = -1.+2.*rndmEngineMod->flat(); }
916  else if (tau >= tausmall) {
917  static const G4double numlim = 0.01;
918  static const G4double onethird = 1./3.;
919  G4double xmeanth, x2meanth;
920  if(tau < numlim) {
921  xmeanth = 1.0 - tau*(1.0 - 0.5*tau);
922  x2meanth= 1.0 - tau*(5.0 - 6.25*tau)*onethird;
923  } else {
924  xmeanth = G4Exp(-tau);
925  x2meanth = (1.+2.*G4Exp(-2.5*tau))*onethird;
926  }
927 
928  // too large step of low-energy particle
929  G4double relloss = 1. - KineticEnergy/currentKinEnergy;
930  static const G4double rellossmax= 0.50;
931  if(relloss > rellossmax) {
932  return SimpleScattering(xmeanth,x2meanth);
933  }
934  // is step extreme small ?
935  G4bool extremesmallstep = false;
936  G4double tsmall = std::min(tlimitmin,lambdalimit);
937  G4double theta0 = 0.;
938  if(trueStepLength > tsmall) {
939  theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
940  } else {
941  theta0 = sqrt(trueStepLength/tsmall)*ComputeTheta0(tsmall,KineticEnergy);
942  extremesmallstep = true ;
943  }
944 
945  static const G4double onesixth = 1./6.;
946  static const G4double theta0max = CLHEP::pi*onesixth;
947  //G4cout << "Theta0= " << theta0 << " theta0max= " << theta0max
948  // << " sqrt(tausmall)= " << sqrt(tausmall) << G4endl;
949 
950  // protection for very small angles
951  G4double theta2 = theta0*theta0;
952 
953  if(theta2 < tausmall) { return cth; }
954 
955  if(theta0 > theta0max) {
956  return SimpleScattering(xmeanth,x2meanth);
957  }
958 
959  G4double x = theta2*(1.0 - theta2/12.);
960  if(theta2 > numlim) {
961  G4double sth = 2*sin(0.5*theta0);
962  x = sth*sth;
963  }
964 
965  // parameter for tail
966  G4double ltau= G4Log(tau);
967  G4double u = extremesmallstep
968  ? G4Exp(G4Log(tsmall/lambda0)*onesixth)
969  : G4Exp(ltau*onesixth);
970  G4double xx = G4Log(lambdaeff/currentRadLength);
971  G4double xsi = coeffc1+u*(coeffc2+coeffc3*u)+coeffc4*xx;
972 
973  // tail should not be too big
974  xsi = std::max(xsi, 1.9);
975  /*
976  if(KineticEnergy > 20*MeV && xsi < 1.6) {
977  G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= "
978  << KineticEnergy/GeV
979  << " !!** c= " << xsi
980  << " **!! length(mm)= " << trueStepLength << " Zeff= " << Zeff
981  << " " << couple->GetMaterial()->GetName()
982  << " tau= " << tau << G4endl;
983  }
984  */
985 
986  G4double c = xsi;
987 
988  if(std::abs(c-3.) < 0.001) { c = 3.001; }
989  else if(std::abs(c-2.) < 0.001) { c = 2.001; }
990 
991  G4double c1 = c-1.;
992 
993  G4double ea = G4Exp(-xsi);
994  G4double eaa = 1.-ea ;
995  G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/eaa;
996  G4double x0 = 1. - xsi*x;
997 
998  // G4cout << " xmean1= " << xmean1 << " xmeanth= " << xmeanth << G4endl;
999 
1000  if(xmean1 <= 0.999*xmeanth) {
1001  return SimpleScattering(xmeanth,x2meanth);
1002  }
1003  //from continuity of derivatives
1004  G4double b = 1.+(c-xsi)*x;
1005 
1006  G4double b1 = b+1.;
1007  G4double bx = c*x;
1008 
1009  G4double eb1 = G4Exp(G4Log(b1)*c1);
1010  G4double ebx = G4Exp(G4Log(bx)*c1);
1011  G4double d = ebx/eb1;
1012 
1013  G4double xmean2 = (x0 + d - (bx - b1*d)/(c-2.))/(1. - d);
1014 
1015  G4double f1x0 = ea/eaa;
1016  G4double f2x0 = c1/(c*(1. - d));
1017  G4double prob = f2x0/(f1x0+f2x0);
1018 
1019  G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
1020 
1021  // sampling of costheta
1022  //G4cout << "c= " << c << " qprob= " << qprob << " eb1= " << eb1
1023  // << " c1= " << c1 << " b1= " << b1 << " bx= " << bx << " eb1= " << eb1
1024  // << G4endl;
1025  if(rndmEngineMod->flat() < qprob)
1026  {
1027  G4double var = 0;
1028  if(rndmEngineMod->flat() < prob) {
1029  cth = 1.+G4Log(ea+rndmEngineMod->flat()*eaa)*x;
1030  } else {
1031  var = (1.0 - d)*rndmEngineMod->flat();
1032  if(var < numlim*d) {
1033  var /= (d*c1);
1034  cth = -1.0 + var*(1.0 - 0.5*var*c)*(2. + (c - xsi)*x);
1035  } else {
1036  cth = 1. + x*(c - xsi - c*G4Exp(-G4Log(var + d)/c1));
1037  }
1038  }
1039  /*
1040  if(KineticEnergy > 5*GeV && cth < 0.9) {
1041  G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= "
1042  << KineticEnergy/GeV
1043  << " 1-cosT= " << 1 - cth
1044  << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
1045  << " tau= " << tau
1046  << " prob= " << prob << " var= " << var << G4endl;
1047  G4cout << " c= " << c << " qprob= " << qprob << " eb1= " << eb1
1048  << " ebx= " << ebx
1049  << " c1= " << c1 << " b= " << b << " b1= " << b1
1050  << " bx= " << bx << " d= " << d
1051  << " ea= " << ea << " eaa= " << eaa << G4endl;
1052  }
1053  */
1054  }
1055  else {
1056  cth = -1.+2.*rndmEngineMod->flat();
1057  /*
1058  if(KineticEnergy > 5*GeV) {
1059  G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= "
1060  << KineticEnergy/GeV
1061  << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
1062  << " qprob= " << qprob << G4endl;
1063  }
1064  */
1065  }
1066  }
1067  return cth ;
1068 }
1069 
1070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1071 
1072 G4double G4UrbanMscModel::ComputeTheta0(G4double trueStepLength,
1073  G4double KineticEnergy)
1074 {
1075  // for all particles take the width of the central part
1076  // from a parametrization similar to the Highland formula
1077  // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
1078  G4double invbetacp = std::sqrt((currentKinEnergy+mass)*(KineticEnergy+mass)/
1079  (currentKinEnergy*(currentKinEnergy+2.*mass)*
1080  KineticEnergy*(KineticEnergy+2.*mass)));
1081  G4double y = trueStepLength/currentRadLength;
1082 
1083  if(particle == positron)
1084  {
1085  static const G4double xl= 0.6;
1086  static const G4double xh= 0.9;
1087  static const G4double e = 113.0;
1088  G4double corr;
1089 
1090  G4double tau = std::sqrt(currentKinEnergy*KineticEnergy)/mass;
1091  G4double x = std::sqrt(tau*(tau+2.)/((tau+1.)*(tau+1.)));
1092  G4double a = 0.994-4.08e-3*Zeff;
1093  G4double b = 7.16+(52.6+365./Zeff)/Zeff;
1094  G4double c = 1.000-4.47e-3*Zeff;
1095  G4double d = 1.21e-3*Zeff;
1096  if(x < xl) {
1097  corr = a*(1.-G4Exp(-b*x));
1098  } else if(x > xh) {
1099  corr = c+d*G4Exp(e*(x-1.));
1100  } else {
1101  G4double yl = a*(1.-G4Exp(-b*xl));
1102  G4double yh = c+d*G4Exp(e*(xh-1.));
1103  G4double y0 = (yh-yl)/(xh-xl);
1104  G4double y1 = yl-y0*xl;
1105  corr = y0*x+y1;
1106  }
1107  //==================================================================
1108  y *= corr*(1.+Zeff*(1.84035e-4*Zeff-1.86427e-2)+0.41125);
1109  }
1110 
1111  static const G4double c_highland = 13.6*CLHEP::MeV;
1112  G4double theta0 = c_highland*std::abs(charge)*std::sqrt(y)*invbetacp;
1113 
1114  // correction factor from e- scattering data
1115  theta0 *= (coeffth1+coeffth2*G4Log(y));
1116  return theta0;
1117 }
1118 
1119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1120 
1121 void G4UrbanMscModel::SampleDisplacement(G4double, G4double phi)
1122 {
1123  // simple and fast sampling
1124  // based on single scattering results
1125  // u = r/rmax : mean value
1126 
1127  G4double rmax = sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
1128  if(rmax > 0.)
1129  {
1130  G4double r = 0.73*rmax;
1131 
1132  // simple distribution for v=Phi-phi=psi ~exp(-beta*v)
1133  // beta determined from the requirement that distribution should give
1134  // the same mean value than that obtained from the ss simulation
1135 
1136  static const G4double cbeta = 2.160;
1137  static const G4double cbeta1 = 1.-G4Exp(-cbeta*CLHEP::pi);
1138  G4double psi = -G4Log(1.-rndmEngineMod->flat()*cbeta1)/cbeta;
1139  G4double Phi = (rndmEngineMod->flat() < 0.5) ? phi+psi : phi-psi;
1140  fDisplacement.set(r*std::cos(Phi),r*std::sin(Phi),0.0);
1141  }
1142 }
1143 
1144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1145 
1146 void G4UrbanMscModel::SampleDisplacementNew(G4double, G4double phi)
1147 {
1148  // best sampling based on single scattering results
1149 
1150  G4double rmax = sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
1151  G4double r = 0.;
1152  G4double u = r/rmax;
1153  static const G4double reps = 5.e-3;
1154 
1155  if(rmax > 0.)
1156  {
1157  static const G4double umax = 0.855;
1158  static const G4double wlow = 0.750;
1159 
1160  static const G4double ralpha = 6.83e+0;
1161  static const G4double ra1 =-4.16179e+1;
1162  static const G4double ra2 = 1.12548e+2;
1163  static const G4double ra3 =-8.66665e+1;
1164  static const G4double ralpha1 = 0.751*ralpha;
1165  static const G4double ralpha2 =ralpha-ralpha1;
1166  static const G4double rwa1 = G4Exp(ralpha1*reps);
1167  static const G4double rwa2 = G4Exp(ralpha1*umax)-rwa1;
1168  static const G4double rejamax = 1.16456;
1169 
1170  static const G4double rbeta = 2.18e+1;
1171  static const G4double rb0 = 4.81382e+2;
1172  static const G4double rb1 =-1.12842e+4;
1173  static const G4double rb2 = 4.57745e+4;
1174  static const G4double rbeta1 = 0.732*rbeta;
1175  static const G4double rbeta2 = rbeta-rbeta1;
1176  static const G4double rwb1 = G4Exp(-rbeta1*umax);
1177  static const G4double rwb2 = rwb1-G4Exp(-rbeta1*(1.-reps));
1178  static const G4double rejbmax = 1.62651;
1179 
1180  G4int count = 0;
1181  G4double uc,rej;
1182 
1183  if(rndmEngineMod->flat() < wlow)
1184  {
1185  do {
1186  u = G4Log(rwa1+rwa2*rndmEngineMod->flat())/ralpha1;
1187  uc = umax-u;
1188  rej = G4Exp(-ralpha2*uc)*
1189  (1.+ralpha*uc+ra1*uc*uc+ra2*uc*uc*uc+ra3*uc*uc*uc*uc);
1190  } while (rejamax*rndmEngineMod->flat() > rej && ++count < 1000);
1191  }
1192  else
1193  {
1194  do {
1195  u = -G4Log(rwb1-rwb2*rndmEngineMod->flat())/rbeta1;
1196  uc = u-umax;
1197  rej = G4Exp(-rbeta2*uc)*
1198  (1.+rbeta*uc+rb0*uc*uc+rb1*uc*uc*uc+rb2*uc*uc*uc*uc);
1199  } while (rejbmax*rndmEngineMod->flat() > rej && ++count < 1000);
1200  }
1201  r = rmax*u;
1202  }
1203 
1204  if(r > 0.)
1205  {
1206  // sample Phi using lateral correlation
1207  // and r/rmax - (Phi-phi) correlation
1208  // v = Phi-phi = acos(latcorr/(r*sth))
1209  // from SS simulation f(v)*g(v)
1210  // f(v) ~ exp(-a1*v) normalized distribution
1211  // g(v) rejection function (0 < g(v) <= 1)
1212  G4double v, rej;
1213 
1214  static const G4double peps = 1.e-4;
1215  static const G4double Pi = CLHEP::pi;
1216  static const G4double palpha[10] = {2.300e+0,2.490e+0,2.610e+0,2.820e+0,2.710e+0,
1217  2.750e+0,2.910e+0,3.400e+0,4.150e+0,5.400e+0};
1218  static const G4double palpha1[10]= {4.600e-2,1.245e-1,2.610e-1,2.820e-1,2.710e-1,
1219  6.875e-1,1.019e+0,1.360e+0,1.660e+0,2.430e+0};
1220  static const G4double pejmax[10] = {3.513,1.968,1.479,1.239,1.116,
1221  1.081,1.064,1.073,1.103,1.158};
1222 
1223  static const G4double pa1[10] = { 3.218e+0, 2.412e+0, 2.715e+0, 2.787e+0, 2.541e+0,
1224  2.508e+0, 2.600e+0, 3.231e+0, 4.588e+0, 6.584e+0};
1225  static const G4double pa2[10] = {-5.528e-1, 2.523e+0, 1.738e+0, 2.082e+0, 1.423e+0,
1226  4.682e-1,-6.883e-1,-2.147e+0,-5.127e+0,-1.054e+1};
1227  static const G4double pa3[10] = { 3.618e+0, 2.032e+0, 2.341e+0, 2.172e+0, 7.205e-1,
1228  4.655e-1, 6.318e-1, 1.255e+0, 2.425e+0, 4.938e+0};
1229  static const G4double pa4[10] = { 2.437e+0, 9.450e-1, 4.349e-1, 2.221e-1, 1.130e-1,
1230  5.405e-2, 2.245e-2, 7.370e-3, 1.456e-3, 1.508e-4};
1231  static const G4double pw1[10] = {G4Exp(-palpha1[0]*peps),G4Exp(-palpha1[1]*peps),
1232  G4Exp(-palpha1[2]*peps),G4Exp(-palpha1[3]*peps),
1233  G4Exp(-palpha1[4]*peps),G4Exp(-palpha1[5]*peps),
1234  G4Exp(-palpha1[6]*peps),G4Exp(-palpha1[7]*peps),
1235  G4Exp(-palpha1[8]*peps),G4Exp(-palpha1[9]*peps)};
1236  static const G4double pw2[10] = {pw1[0]-G4Exp(-palpha1[0]*(Pi-peps)),
1237  pw1[1]-G4Exp(-palpha1[1]*(Pi-peps)),
1238  pw1[2]-G4Exp(-palpha1[2]*(Pi-peps)),
1239  pw1[3]-G4Exp(-palpha1[3]*(Pi-peps)),
1240  pw1[4]-G4Exp(-palpha1[4]*(Pi-peps)),
1241  pw1[5]-G4Exp(-palpha1[5]*(Pi-peps)),
1242  pw1[6]-G4Exp(-palpha1[6]*(Pi-peps)),
1243  pw1[7]-G4Exp(-palpha1[7]*(Pi-peps)),
1244  pw1[8]-G4Exp(-palpha1[8]*(Pi-peps)),
1245  pw1[9]-G4Exp(-palpha1[9]*(Pi-peps))};
1246 
1247  G4int iphi = u*10.;
1248  if(iphi < 0) { iphi = 0; }
1249  else if(iphi > 9) { iphi = 9; }
1250  G4int count = 0;
1251 
1252  do {
1253  v = -G4Log(pw1[iphi]-pw2[iphi]*rndmEngineMod->flat())/palpha1[iphi];
1254  rej = (G4Exp(-palpha[iphi]*v)*
1255  (1+pa1[iphi]*v+pa2[iphi]*v*v+pa3[iphi]*v*v*v)+pa4[iphi])/
1256  G4Exp(-pw1[iphi]*v);
1257  }
1258  // Loop checking, 5-March-2018, Vladimir Ivanchenko
1259  while (pejmax[iphi]*rndmEngineMod->flat() > rej && ++count < 1000);
1260 
1261  G4double Phi = (rndmEngineMod->flat() < 0.5) ? phi+v : phi-v;
1262  fDisplacement.set(r*std::cos(Phi),r*std::sin(Phi),0.0);
1263  }
1264 }
1265 
1266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......