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