Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4UrbanAdjointMscModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// -------------------------------------------------------------------
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
59using namespace std;
60
61//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62
63G4UrbanAdjointMscModel::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......
128G4UrbanAdjointMscModel::~G4UrbanAdjointMscModel() {}
129
130//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131void 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
148G4double 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......
410void 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......
424G4double 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......
731G4double 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......
794G4double 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......
850G4ThreeVector& 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......
904G4double 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......
1070G4double 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......
1123void 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......
1183void 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......
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
fUseSafety
@ fUseSafety
Definition: G4MscStepLimitType.hh:49
fUseSafetyPlus
@ fUseSafetyPlus
Definition: G4MscStepLimitType.hh:50
fUseDistanceToBoundary
@ fUseDistanceToBoundary
Definition: G4MscStepLimitType.hh:51
G4StepStatus
G4StepStatus
Definition: G4StepStatus.hh:40
fGeomBoundary
@ fGeomBoundary
Definition: G4StepStatus.hh:43
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:33
CLHEP::HepRandomEngine::flat
virtual double flat()=0
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4IonisParamMat::GetZeffective
G4double GetZeffective() const
Definition: G4IonisParamMat.hh:139
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:87
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
G4MaterialCutsCouple::GetProductionCuts
G4ProductionCuts * GetProductionCuts() const
Definition: G4MaterialCutsCouple.hh:141
G4Material::GetIonisation
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:221
G4Material::GetRadlen
G4double GetRadlen() const
Definition: G4Material.hh:215
G4ParticleChangeForMSC::ProposeMomentumDirection
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
Definition: G4ParticleChangeForMSC.hh:78
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:109
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:93
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4Pow::Z23
G4double Z23(G4int Z) const
Definition: G4Pow.hh:125
G4ProductionCuts::GetProductionCut
G4double GetProductionCut(G4int index) const
Definition: G4ProductionCuts.cc:116
G4Step::GetPreStepPoint
G4StepPoint * GetPreStepPoint() const
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4Track::GetStep
const G4Step * GetStep() const
G4UrbanAdjointMscModel::ComputeCrossSectionPerAtom
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=0., G4double emax=DBL_MAX) override
Definition: G4UrbanAdjointMscModel.cc:148
G4UrbanAdjointMscModel::G4UrbanAdjointMscModel
G4UrbanAdjointMscModel(const G4String &nam="UrbanMsc")
Definition: G4UrbanAdjointMscModel.cc:63
G4UrbanAdjointMscModel::ComputeTrueStepLength
G4double ComputeTrueStepLength(G4double geomStepLength) override
Definition: G4UrbanAdjointMscModel.cc:794
G4UrbanAdjointMscModel::StartTracking
void StartTracking(G4Track *) override
Definition: G4UrbanAdjointMscModel.cc:410
G4UrbanAdjointMscModel::ComputeTheta0
G4double ComputeTheta0(G4double truePathLength, G4double KineticEnergy)
Definition: G4UrbanAdjointMscModel.cc:1070
G4UrbanAdjointMscModel::ComputeGeomPathLength
G4double ComputeGeomPathLength(G4double truePathLength) override
Definition: G4UrbanAdjointMscModel.cc:731
G4UrbanAdjointMscModel::SampleScattering
G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety) override
Definition: G4UrbanAdjointMscModel.cc:850
G4UrbanAdjointMscModel::ComputeTruePathLengthLimit
G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep) override
Definition: G4UrbanAdjointMscModel.cc:424
G4UrbanAdjointMscModel::Initialise
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: G4UrbanAdjointMscModel.cc:131
G4VEmModel::SetCurrentCouple
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:468
G4VMscModel::dtrl
G4double dtrl
Definition: G4VMscModel.hh:203
G4VMscModel::GetDEDX
G4double GetDEDX(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.cc:158
G4VMscModel::facrange
G4double facrange
Definition: G4VMscModel.hh:199
G4VMscModel::ComputeGeomLimit
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:296
G4VMscModel::skin
G4double skin
Definition: G4VMscModel.hh:202
G4VMscModel::GetTransportMeanFreePath
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:325
G4VMscModel::GetParticleChangeForMSC
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=nullptr)
Definition: G4VMscModel.cc:77
G4VMscModel::GetEnergy
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.cc:223
G4VMscModel::GetRange
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.cc:188
G4VMscModel::steppingAlgorithm
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:209
G4VMscModel::ConvertTrueToGeom
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:286
G4VMscModel::latDisplasment
G4bool latDisplasment
Definition: G4VMscModel.hh:212
G4VMscModel::ComputeSafety
G4double ComputeSafety(const G4ThreeVector &position, G4double limit=DBL_MAX)
Definition: G4VMscModel.hh:278
G4VMscModel::facsafety
G4double facsafety
Definition: G4VMscModel.hh:201
G4VMscModel::fDisplacement
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:208
G4VMscModel::facgeom
G4double facgeom
Definition: G4VMscModel.hh:200
field_utils::Value1D::KineticEnergy
@ KineticEnergy
G4lrint
int G4lrint(double ad)
Definition: templates.hh:134