| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 sca
37 // H.W.Lewis Phys Rev 78 (1950) 526 and others
38 // In its present form the model can be used
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........oo
58
59 using namespace std;
60
61 //....oooOO0OOooo........oooOO0OOooo........oo
62
63 G4UrbanAdjointMscModel::G4UrbanAdjointMscModel
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 = lambda
120 zPathLength = par1 = par2 = par3 = 0.;
121
122 currentMaterialIndex = -1;
123 fParticleChange = nullptr;
124 couple = nullptr;
125 }
126
127 //....oooOO0OOooo........oooOO0OOooo........oo
128 G4UrbanAdjointMscModel::~G4UrbanAdjointMscMode
129
130 //....oooOO0OOooo........oooOO0OOooo........oo
131 void G4UrbanAdjointMscModel::Initialise(const
132 const
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........oo
147
148 G4double G4UrbanAdjointMscModel::ComputeCrossS
149 const G4ParticleDefinition* part, G4double K
150 G4double AtomicNumber, G4double, G4double, G
151 {
152 static constexpr G4double epsmin = 1.e-4;
153 static constexpr G4double epsmax = 1.e10;
154
155 static constexpr G4double Zdat[15] = { 4.,
156 47.,
157
158 // corr. factors for e-/e+ lambda for T <= T
159 static constexpr G4double celectron[15][22]
160 { 1.125, 1.072, 1.051, 1.047, 1.047, 1.050
161 1.054, 1.057, 1.062, 1.069, 1.075, 1.090
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
164 1.052, 1.053, 1.058, 1.065, 1.072, 1.087
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
167 1.136, 1.114, 1.106, 1.106, 1.109, 1.119
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
170 1.190, 1.133, 1.107, 1.099, 1.098, 1.103
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
173 1.203, 1.122, 1.080, 1.065, 1.061, 1.063
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
176 1.256, 1.155, 1.099, 1.077, 1.070, 1.068
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
179 1.300, 1.180, 1.112, 1.082, 1.073, 1.066
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
182 1.339, 1.195, 1.108, 1.068, 1.053, 1.040
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
185 1.498, 1.301, 1.171, 1.105, 1.077, 1.048
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
188 1.528, 1.319, 1.178, 1.106, 1.075, 1.040
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
191 1.569, 1.342, 1.186, 1.102, 1.065, 1.022
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
194 1.677, 1.410, 1.224, 1.121, 1.073, 1.014
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
197 1.837, 1.512, 1.283, 1.153, 1.091, 1.010
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
200 1.939, 1.579, 1.325, 1.178, 1.108, 1.014
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
203 1.985, 1.609, 1.343, 1.188, 1.113, 1.013
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
209 1.097, 1.083, 1.080, 1.086, 1.092, 1.108
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
212 1.122, 1.096, 1.089, 1.092, 1.098, 1.114
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
215 1.357, 1.246, 1.194, 1.179, 1.178, 1.188
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
218 1.553, 1.353, 1.253, 1.219, 1.211, 1.214
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
221 1.624, 1.382, 1.259, 1.214, 1.202, 1.202
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
224 1.759, 1.465, 1.311, 1.252, 1.234, 1.228
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
227 1.880, 1.535, 1.353, 1.281, 1.258, 1.247
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
230 2.015, 1.602, 1.385, 1.297, 1.268, 1.251
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
233 2.407, 1.873, 1.564, 1.425, 1.374, 1.330
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
236 2.490, 1.925, 1.596, 1.447, 1.391, 1.342
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
239 2.602, 1.995, 1.641, 1.477, 1.414, 1.356
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
242 2.926, 2.188, 1.763, 1.563, 1.484, 1.405
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
245 3.417, 2.478, 1.944, 1.692, 1.589, 1.480
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
248 3.641, 2.648, 2.064, 1.779, 1.661, 1.531
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
251 3.752, 2.724, 2.116, 1.817, 1.692, 1.554
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
257 79.
258 41.
259 -7.
260
261 G4double sigma;
262 SetParticle(part);
263
264 Z23 = G4Pow::GetInstance()->Z23(G4lrint(Atom
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.) /
276 G4double w = c - 2.;
277 G4double tau = 0.5 * (w + sqrt(w * w + 4.
278 eKineticEnergy = electron_mass_c2 * tau;
279 }
280
281 G4double eTotalEnergy = eKineticEnergy + ele
282 G4double beta2 = eKineticEnergy * (eT
283 (eTotalEnergy * eTotalEnerg
284 G4double bg2 = eKineticEnergy * (eTotalEnerg
285 (electron_mass_c2 * electron_
286
287 static constexpr G4double epsfactor =
288 2. * CLHEP::electron_mass_c2 * CLHEP::elec
289 CLHEP::Bohr_radius * CLHEP::Bohr_radius /
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 /
296 else
297 sigma = G4Log(2. * eps) - 1. + 1. / eps;
298
299 sigma *= ChargeSquare * AtomicNumber * Atomi
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
317
318 static constexpr G4double Tlim = 10. * CLHEP
319 static constexpr G4double sigmafactor =
320 CLHEP::twopi * CLHEP::classic_electr_radiu
321 static const G4double beta2lim =
322 Tlim * (Tlim + 2. * CLHEP::electron_mass_c
323 ((Tlim + CLHEP::electron_mass_c2) * (Tlim
324 static const G4double bg2lim =
325 Tlim * (Tlim + 2. * CLHEP::electron_mass_c
326 (CLHEP::electron_mass_c2 * CLHEP::electron
327
328 static constexpr G4double sig0[15] = {
329 0.2672 * CLHEP::barn, 0.5922 * CLHEP::barn
330 6.235 * CLHEP::barn, 11.69 * CLHEP::barn,
331 16.12 * CLHEP::barn, 23.00 * CLHEP::barn,
332 39.95 * CLHEP::barn, 50.85 * CLHEP::barn,
333 91.15 * CLHEP::barn, 104.4 * CLHEP::barn,
334 };
335
336 static constexpr G4double Tdat[22] = {
337 100 * CLHEP::eV, 200 * CLHEP::eV, 400 *
338 1 * CLHEP::keV, 2 * CLHEP::keV, 4 * CL
339 10 * CLHEP::keV, 20 * CLHEP::keV, 40 * C
340 100 * CLHEP::keV, 200 * CLHEP::keV, 400 *
341 1 * CLHEP::MeV, 2 * CLHEP::MeV, 4 * CL
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] >= eKineticEn
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_ma
358 G4double b2small = T * (E + electron_mass_
359
360 T = Tdat[iT + 1];
361 E = T + electron_mass_c2;
362 G4double b2big = T * (E + electron_mass_c2
363 G4double ratb2 = (beta2 - b2small) / (b2bi
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]
397 c2 =
398 bg2lim * sig0[iZ + 1] * (1. + hecorr[iZ
399 if((AtomicNumber >= ZZ1) && (AtomicNumber
400 sigma = c1 + ratZ * (c2 - c1);
401 else if(AtomicNumber < ZZ1)
402 sigma = AtomicNumber * AtomicNumber * c1
403 else if(AtomicNumber > ZZ2)
404 sigma = AtomicNumber * AtomicNumber * c2
405 }
406 return sigma;
407 }
408
409 //....oooOO0OOooo........oooOO0OOooo........oo
410 void G4UrbanAdjointMscModel::StartTracking(G4T
411 {
412 SetParticle(track->GetDynamicParticle()->Get
413 firstStep = true;
414 insideskin = false;
415 fr = facrange;
416 tlimit = tgeom = rangeinit = rangecut = geom
417 smallstep = 1.e1
418 stepmin = tlim
419 tlimitmin = 10.
420 rndmEngineMod = G4Ra
421 }
422
423 //....oooOO0OOooo........oooOO0OOooo........oo
424 G4double G4UrbanAdjointMscModel::ComputeTruePa
425 const G4Track& track, G4double& currentMinim
426 {
427 tPathLength = currentMinimal
428 const G4DynamicParticle* dp = track.GetDynam
429
430 G4StepPoint* sp = track.GetStep()->G
431 G4StepStatus stepStatus = sp->GetStepStatus(
432 couple = track.GetMaterialC
433 SetCurrentCouple(couple);
434 currentMaterialIndex = couple->GetIndex();
435 currentKinEnergy = dp->GetKineticEnergy(
436
437 currentRange = GetRange(particle, currentKin
438 lambda0 = GetTransportMeanFreePath(part
439 tPathLength = min(tPathLength, currentRange
440
441 // set flag to default values
442 Zeff = couple->GetMaterial()->GetIonisation(
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, curr
452 }
453
454 // upper limit for the straight line distanc
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.22
465 }
466 presafety = sp->GetSafety();
467
468 // far from geometry boundary
469 if(distance < presafety)
470 {
471 latDisplasment = false;
472 return ConvertTrueToGeom(tPathLength, curr
473 }
474
475 latDisplasment = latDispla
476 static constexpr G4double invmev = 1.0 / CLH
477
478 // standard version
479 if(steppingAlgorithm == fUseDistanceToBounda
480 {
481 // compute geomlimit and presafety
482 geomlimit = ComputeGeomLimit(track, presaf
483
484 // is it far from boundary ?
485 if(distance < presafety)
486 {
487 latDisplasment = false;
488 return ConvertTrueToGeom(tPathLength, cu
489 }
490
491 smallstep += 1.;
492 insideskin = false;
493
494 // initialisation at firs step and at the
495 if(firstStep || (stepStatus == fGeomBounda
496 {
497 rangeinit = currentRange;
498 if(!firstStep)
499 {
500 smallstep = 1.;
501 }
502
503 // define stepmin here (it depends on la
504 // rough estimation of lambda_elastic/la
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, tlimitminf
511
512 // constraint from the geometry
513 if((geomlimit < geombig) && (geomlimit >
514 {
515 // geomlimit is a geometrical step len
516 // transform it to true path length (e
517 if((1. - geomlimit / lambda0) > 0.)
518 geomlimit = -lambda0 * G4Log(1. - ge
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
538 (smallstep > skin) && (tPathLength < ge
539 {
540 return ConvertTrueToGeom(tPathLength, cu
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
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
565 if((tlimit < tPathLength) && (smallstep >
566 {
567 tPathLength = min(tPathLength, Randomize
568 }
569 else
570 {
571 tPathLength = min(tPathLength, tlimit);
572 }
573 }
574 // for 'normal' simulation with or without m
575 // there no small step/single scattering at
576 else if(steppingAlgorithm == fUseSafety)
577 {
578 if(stepStatus != fGeomBoundary)
579 {
580 presafety = ComputeSafety(sp->GetPositio
581 }
582 // is far from boundary
583 if(distance < presafety)
584 {
585 latDisplasment = false;
586 return ConvertTrueToGeom(tPathLength, cu
587 }
588
589 if(firstStep || (stepStatus == fGeomBounda
590 {
591 rangeinit = currentRange;
592 fr = facrange;
593 // Geant4 version 9.1 like stepping for
594 if(mass < masslimite)
595 {
596 rangeinit = max(rangeinit, lambda0);
597 if(lambda0 > lambdalimit)
598 {
599 fr *= (0.75 + 0.25 * lambda0 / lambd
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, tlimitm
607 }
608
609 // step limit
610 tlimit = max(fr * rangeinit, facsafety * p
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, Randomize
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->GetPositio
631 }
632
633 // is far from boundary
634 if(distance < presafety)
635 {
636 latDisplasment = false;
637 return ConvertTrueToGeom(tPathLength, cu
638 }
639
640 if(firstStep || (stepStatus == fGeomBounda
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()
651 if(lambda0 > lambdalimit)
652 {
653 fr *= (0.84 + 0.16 * lambda0 / lambd
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, tlimitm
661 }
662 // step limit
663 tlimit = max(fr * rangeinit, facsafety * p
664
665 // lower limit for tlimit
666 tlimit = max(tlimit, tlimitmin);
667
668 // condition for tPathLength from drr and
669 if(currentRange > finalr)
670 {
671 G4double tmax =
672 drr * currentRange + finalr * (1. - dr
673 tPathLength = min(tPathLength, tmax);
674 }
675
676 // condition safety
677 if(currentRange > rangecut)
678 {
679 if(firstStep)
680 {
681 tPathLength = min(tPathLength, facsafe
682 }
683 else if(stepStatus != fGeomBoundary && p
684 {
685 tPathLength = min(tPathLength, presafe
686 }
687 }
688
689 // randomise if step determined by msc
690 if(tPathLength < tlimit)
691 {
692 tPathLength = min(tPathLength, Randomize
693 }
694 else
695 {
696 tPathLength = min(tPathLength, tlimit);
697 }
698 }
699
700 // version similar to 7.1 (needed for some e
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, Randomize
720 }
721 else
722 {
723 tPathLength = min(tPathLength, tlimit);
724 }
725 }
726 firstStep = false;
727 return ConvertTrueToGeom(tPathLength, curren
728 }
729
730 //....oooOO0OOooo........oooOO0OOooo........oo
731 G4double G4UrbanAdjointMscModel::ComputeGeomPa
732 {
733 lambdaeff = lambda0;
734 par1 = -1.;
735 par2 = par3 = 0.;
736
737 // this correction needed to run MSC with eI
738 // and makes no harm for a normal run
739 tPathLength = std::min(tPathLength, currentR
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 *
758 else
759 zPathLength = lambda0 * (1. - G4Exp(-tau
760 }
761 else if(currentKinEnergy < mass || tPathLeng
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. - tPathLen
770 (par1 * par3);
771 }
772 else
773 {
774 zPathLength = 1. / (par1 * par3);
775 }
776 }
777 else
778 {
779 G4double rfin = max(currentRange - tPat
780 G4double T1 = GetEnergy(particle, rfi
781 G4double lambda1 = GetTransportMeanFreePat
782
783 par1 = (lambda0 - lambda1) / (lambd
784 par2 = 1. / (par1 * lambda0);
785 par3 = 1. + par2;
786 zPathLength = (1. - G4Exp(par3 * G4Log(lam
787 }
788
789 zPathLength = min(zPathLength, lambda0);
790 return zPathLength;
791 }
792
793 //....oooOO0OOooo........oooOO0OOooo........oo
794 G4double G4UrbanAdjointMscModel::ComputeTrueSt
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) &
815 {
816 if(par1 < 0.)
817 {
818 tlength = -lambda0 * G4Log(1. - geomSt
819 }
820 else
821 {
822 if(par1 * par3 * geomStepLength < 1.)
823 {
824 tlength =
825 (1. - G4Exp(G4Log(1. - par1 * 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........oo
850 G4ThreeVector& G4UrbanAdjointMscModel::SampleS
851 const G4ThreeVector& oldDirection, G4double
852 {
853 fDisplacement.set(0.0, 0.0, 0.0);
854 G4double kineticEnergy = currentKinEnergy;
855 if(tPathLength > currentRange * dtrl)
856 {
857 kineticEnergy = GetEnergy(particle, curren
858 }
859 else
860 {
861 kineticEnergy -= tPathLength * GetDEDX(par
862 }
863
864 if((kineticEnergy <= eV) || (tPathLength <=
865 (tPathLength < tausmall * lambda0))
866 {
867 return fDisplacement;
868 }
869
870 G4double cth = SampleCosineTheta(tPathLength
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 + ct
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(ne
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........oo
904 G4double G4UrbanAdjointMscModel::SampleCosineT
905
906 {
907 G4double cth = 1.;
908 G4double tau = trueStepLength / lambda0;
909 currentTau = tau;
910 lambdaeff = lambda0;
911
912 G4double lambda1 = GetTransportMeanFreePath(
913 if(std::fabs(lambda1 - lambda0) > lambda0 *
914 {
915 // mean tau value
916 tau = trueStepLength * G4Log(lambda0 / lam
917 }
918
919 currentTau = tau;
920 lambdaeff = trueStepLength / currentT
921 currentRadLength = couple->GetMaterial()->Ge
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
935 }
936 else
937 {
938 xmeanth = G4Exp(-tau);
939 x2meanth = (1. + 2. * G4Exp(-2.5 * tau))
940 }
941
942 // too large step of low-energy particle
943 G4double relloss = 1. - Ki
944 static const G4double rellossmax = 0.50;
945 if(relloss > rellossmax)
946 {
947 return SimpleScattering(xmeanth, x2meant
948 }
949 // is step extreme small ?
950 G4bool extremesmallstep = false;
951 G4double tsmall = std::min(tlimitm
952 G4double theta0 = 0.;
953 if(trueStepLength > tsmall)
954 {
955 theta0 = ComputeTheta0(trueStepLength, K
956 }
957 else
958 {
959 theta0 =
960 sqrt(trueStepLength / tsmall) * Comput
961 extremesmallstep = true;
962 }
963
964 static constexpr G4double theta0max = CLHE
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, x2meant
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 / currentRa
992 G4double xsi = coeffc1 + u * (coeffc2 + co
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) *
1016 G4double x0 = 1. - xsi * x;
1017
1018 if(xmean1 <= 0.999 * xmeanth)
1019 {
1020 return SimpleScattering(xmeanth, x2mean
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)
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
1039
1040 if(rndmEngineMod->flat() < qprob)
1041 {
1042 G4double var = 0;
1043 if(rndmEngineMod->flat() < prob)
1044 {
1045 cth = 1. + G4Log(ea + rndmEngineMod->
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
1054 }
1055 else
1056 {
1057 cth = 1. + x * (c - xsi - c * G4Exp
1058 }
1059 }
1060 }
1061 else
1062 {
1063 cth = -1. + 2. * rndmEngineMod->flat();
1064 }
1065 }
1066 return cth;
1067 }
1068
1069 //....oooOO0OOooo........oooOO0OOooo........o
1070 G4double G4UrbanAdjointMscModel::ComputeTheta
1071
1072 {
1073 // for all particles take the width of the
1074 // from a parametrization similar to the H
1075 // ( Highland formula: Particle Physics Boo
1076 G4double invbetacp =
1077 std::sqrt((currentKinEnergy + mass) * (Ki
1078 (currentKinEnergy * (currentKin
1079 KineticEnergy * (KineticEnergy
1080 G4double y = trueStepLength / currentRadLen
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
1090 G4double x = std::sqrt(tau * (tau + 2.)
1091 G4double a = 0.994 - 4.08e-3 * Zeff;
1092 G4double b = 7.16 + (52.6 + 365. / 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 * Ze
1112 }
1113
1114 static constexpr G4double c_highland = 13.6
1115 G4double theta0 = c_highland * std::abs(cha
1116
1117 // correction factor from e- scattering dat
1118 theta0 *= (coeffth1 + coeffth2 * G4Log(y));
1119 return theta0;
1120 }
1121
1122 //....oooOO0OOooo........oooOO0OOooo........o
1123 void G4UrbanAdjointMscModel::SampleDisplaceme
1124 {
1125 G4double rmax =
1126 sqrt((tPathLength - zPathLength) * (tPath
1127
1128 static constexpr G4double third = 1. / 3.;
1129 G4double r = rmax * G4Exp(G4Log(rndmEngineM
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) && !insideski
1138 {
1139 if(currentTau < taulim)
1140 {
1141 latcorr = lambdaeff * kappa * current
1142 (1. - (kappa + 1.) * curren
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 / kappam
1154 latcorr *= 2. * lambdaeff * third;
1155 }
1156 }
1157 latcorr = std::min(latcorr, r);
1158
1159 // sample direction of lateral displaceme
1160 // compute it from the lateral correlatio
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 *
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 *
1179 }
1180 }
1181
1182 //....oooOO0OOooo........oooOO0OOooo........o
1183 void G4UrbanAdjointMscModel::SampleDisplaceme
1184 {
1185 // sample displacement r
1186
1187 G4double rmax =
1188 sqrt((tPathLength - zPathLength) * (tPath
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) * rndmEngin
1205 v = 1. - u;
1206 rej = rp0 * G4Exp(rp1 * G4Log(v) - rp2 *
1207 } while(rndmEngineMod->flat() > rej && ++co
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 c
1215 // simulation as
1216 // f(v) = 1.49e-2*exp(-v**2/(2*0.320))+2.
1217 // 1.96e-5*exp(8.42e-1*log(1.+1.45
1218 static constexpr G4double probv1 = 0.3055
1219 static constexpr G4double probv2 = 0.9551
1220 static constexpr G4double vhigh = 3.15;
1221 static const G4double w2v = 1. / G4Exp(30
1222 static const G4double w3v = 1. / G4Exp(-1
1223
1224 G4double Phi;
1225 G4double random = rndmEngineMod->flat();
1226 if(random < probv1)
1227 {
1228 do
1229 {
1230 v = G4RandGauss::shoot(rndmEngineMod,
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. - rndmEngine
1240 6.30e-2;
1241 }
1242 else
1243 {
1244 v = (-1. + 1. / G4Exp(G4Log(1. - rndm
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 *
1260 }
1261 }
1262
1263 //....oooOO0OOooo........oooOO0OOooo........o
1264