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