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