Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4UrbanMscModel90.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// $Id$
27//
28// -------------------------------------------------------------------
29//
30// GEANT4 Class file
31//
32//
33// File name: G4UrbanMscModel90
34//
35// Author: V.Ivanchenko clone Laszlo Urban model
36//
37// Creation date: 07.12.2007
38//
39// Modifications:
40//
41//
42
43// Class Description:
44//
45// Implementation of the model of multiple scattering based on
46// H.W.Lewis Phys Rev 78 (1950) 526 and others
47
48// -------------------------------------------------------------------
49//
50
51
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54
55#include "G4UrbanMscModel90.hh"
57#include "G4SystemOfUnits.hh"
58#include "Randomize.hh"
59#include "G4Electron.hh"
60
61#include "G4LossTableManager.hh"
63
64#include "G4Poisson.hh"
65
66//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67
68using namespace std;
69
71 : G4VMscModel(nam)
72{
73 taubig = 8.0;
74 tausmall = 1.e-20;
75 taulim = 1.e-6;
76 currentTau = taulim;
77 tlimitminfix = 1.e-6*mm;
78 stepmin = tlimitminfix;
79 smallstep = 1.e10;
80 currentRange = 0. ;
81 frscaling2 = 0.25;
82 frscaling1 = 1.-frscaling2;
83 tlimit = 1.e10*mm;
84 tlimitmin = 10.*tlimitminfix;
85 nstepmax = 25.;
86 geombig = 1.e50*mm;
87 geommin = 1.e-3*mm;
88 geomlimit = geombig;
89 presafety = 0.*mm;
90 Zeff = 1.;
91 particle = 0;
92 theManager = G4LossTableManager::Instance();
93 firstStep = true;
94 inside = false;
95 insideskin = false;
96
97 skindepth = skin*stepmin;
98
99 mass = proton_mass_c2;
100 charge = 1.0;
101 currentKinEnergy = currentRange = currentRadLength = masslimite = masslimitmu
102 = lambda0 = lambdaeff = tPathLength = zPathLength = par1 = par2 = par3 = 0;
103
104 currentMaterialIndex = -1;
105 fParticleChange = 0;
106 couple = 0;
107 G4cout << "### G4UrbanMscModel90 is obsolete and will be removed for "
108 << "the next Geant4 version" << G4endl;
109}
110
111//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112
114{}
115
116//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117
119 const G4DataVector&)
120{
121 skindepth = skin*stepmin;
122
123 // set values of some data members
124 SetParticle(p);
125
126 if(p->GetPDGMass() < MeV) {
127 G4cout << "### WARNING: G4UrbanMscModel90 model is used for "
128 << p->GetParticleName() << " !!! " << G4endl;
129 G4cout << "### This model should be used only for heavy particles"
130 << G4endl;
131 }
132
133
134 fParticleChange = GetParticleChangeForMSC(p);
135}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
138
140 const G4ParticleDefinition* part,
141 G4double KineticEnergy,
142 G4double AtomicNumber,G4double,
144{
145 const G4double sigmafactor = twopi*classic_electr_radius*classic_electr_radius;
146 const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2*
147 Bohr_radius*Bohr_radius/(hbarc*hbarc);
148 const G4double epsmin = 1.e-4 , epsmax = 1.e10;
149
150 const G4double Zdat[15] = { 4., 6., 13., 20., 26., 29., 32., 38., 47.,
151 50., 56., 64., 74., 79., 82. };
152
153 const G4double Tdat[22] = { 100*eV, 200*eV, 400*eV, 700*eV,
154 1*keV, 2*keV, 4*keV, 7*keV,
155 10*keV, 20*keV, 40*keV, 70*keV,
156 100*keV, 200*keV, 400*keV, 700*keV,
157 1*MeV, 2*MeV, 4*MeV, 7*MeV,
158 10*MeV, 20*MeV};
159
160 // corr. factors for e-/e+ lambda for T <= Tlim
161 G4double celectron[15][22] =
162 {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
163 1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
164 1.112,1.108,1.100,1.093,1.089,1.087 },
165 {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
166 1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
167 1.109,1.105,1.097,1.090,1.086,1.082 },
168 {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
169 1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
170 1.131,1.124,1.113,1.104,1.099,1.098 },
171 {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
172 1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
173 1.112,1.105,1.096,1.089,1.085,1.098 },
174 {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
175 1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
176 1.073,1.070,1.064,1.059,1.056,1.056 },
177 {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
178 1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
179 1.074,1.070,1.063,1.059,1.056,1.052 },
180 {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
181 1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
182 1.068,1.064,1.059,1.054,1.051,1.050 },
183 {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
184 1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
185 1.039,1.037,1.034,1.031,1.030,1.036 },
186 {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
187 1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
188 1.031,1.028,1.024,1.022,1.021,1.024 },
189 {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
190 1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
191 1.020,1.017,1.015,1.013,1.013,1.020 },
192 {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
193 1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
194 0.995,0.993,0.993,0.993,0.993,1.011 },
195 {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
196 1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
197 0.974,0.972,0.973,0.974,0.975,0.987 },
198 {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
199 1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
200 0.950,0.947,0.949,0.952,0.954,0.963 },
201 {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
202 1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
203 0.941,0.938,0.940,0.944,0.946,0.954 },
204 {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
205 1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
206 0.933,0.930,0.933,0.936,0.939,0.949 }};
207
208 G4double cpositron[15][22] = {
209 {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
210 1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
211 1.131,1.126,1.117,1.108,1.103,1.100 },
212 {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
213 1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
214 1.138,1.132,1.122,1.113,1.108,1.102 },
215 {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
216 1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
217 1.203,1.190,1.173,1.159,1.151,1.145 },
218 {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
219 1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
220 1.225,1.210,1.191,1.175,1.166,1.174 },
221 {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
222 1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
223 1.217,1.203,1.184,1.169,1.160,1.151 },
224 {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
225 1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
226 1.237,1.222,1.201,1.184,1.174,1.159 },
227 {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
228 1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
229 1.252,1.234,1.212,1.194,1.183,1.170 },
230 {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
231 2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
232 1.254,1.237,1.214,1.195,1.185,1.179 },
233 {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
234 2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
235 1.312,1.288,1.258,1.235,1.221,1.205 },
236 {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
237 2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
238 1.320,1.294,1.264,1.240,1.226,1.214 },
239 {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
240 2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
241 1.328,1.302,1.270,1.245,1.231,1.233 },
242 {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
243 2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
244 1.361,1.330,1.294,1.267,1.251,1.239 },
245 {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
246 3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
247 1.409,1.372,1.330,1.298,1.280,1.258 },
248 {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
249 3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
250 1.442,1.400,1.354,1.319,1.299,1.272 },
251 {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
252 3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
253 1.456,1.412,1.364,1.328,1.307,1.282 }};
254
255 //data/corrections for T > Tlim
256 G4double Tlim = 10.*MeV;
257 G4double beta2lim = Tlim*(Tlim+2.*electron_mass_c2)/
258 ((Tlim+electron_mass_c2)*(Tlim+electron_mass_c2));
259 G4double bg2lim = Tlim*(Tlim+2.*electron_mass_c2)/
260 (electron_mass_c2*electron_mass_c2);
261
262 G4double sig0[15] = {0.2672*barn, 0.5922*barn, 2.653*barn, 6.235*barn,
263 11.69*barn , 13.24*barn , 16.12*barn, 23.00*barn ,
264 35.13*barn , 39.95*barn , 50.85*barn, 67.19*barn ,
265 91.15*barn , 104.4*barn , 113.1*barn};
266
267 G4double hecorr[15] = {120.70, 117.50, 105.00, 92.92, 79.23, 74.510, 68.29,
268 57.39, 41.97, 36.14, 24.53, 10.21, -7.855, -16.84,
269 -22.30};
270
271 G4double sigma;
272 SetParticle(part);
273
274 G4double Z23 = 2.*log(AtomicNumber)/3.; Z23 = exp(Z23);
275
276 // correction if particle .ne. e-/e+
277 // compute equivalent kinetic energy
278 // lambda depends on p*beta ....
279
280 G4double eKineticEnergy = KineticEnergy;
281
282 if(mass > electron_mass_c2)
283 {
284 G4double TAU = KineticEnergy/mass ;
285 G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
286 G4double w = c-2. ;
287 G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
288 eKineticEnergy = electron_mass_c2*tau ;
289 }
290
291 G4double ChargeSquare = charge*charge;
292
293 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
294 G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
295 /(eTotalEnergy*eTotalEnergy);
296 G4double bg2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
297 /(electron_mass_c2*electron_mass_c2);
298
299 G4double eps = epsfactor*bg2/Z23;
300
301 if (eps<epsmin) sigma = 2.*eps*eps;
302 else if(eps<epsmax) sigma = log(1.+2.*eps)-2.*eps/(1.+2.*eps);
303 else sigma = log(2.*eps)-1.+1./eps;
304
305 sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
306
307 // interpolate in AtomicNumber and beta2
308 G4double c1,c2,cc1,cc2,corr;
309
310 // get bin number in Z
311 G4int iZ = 14;
312 while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
313 if (iZ==14) iZ = 13;
314 if (iZ==-1) iZ = 0 ;
315
316 G4double Z1 = Zdat[iZ];
317 G4double Z2 = Zdat[iZ+1];
318 G4double ratZ = (AtomicNumber-Z1)*(AtomicNumber+Z1)/
319 ((Z2-Z1)*(Z2+Z1));
320
321 if(eKineticEnergy <= Tlim)
322 {
323 // get bin number in T (beta2)
324 G4int iT = 21;
325 while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
326 if(iT==21) iT = 20;
327 if(iT==-1) iT = 0 ;
328
329 // calculate betasquare values
330 G4double T = Tdat[iT], E = T + electron_mass_c2;
331 G4double b2small = T*(E+electron_mass_c2)/(E*E);
332
333 T = Tdat[iT+1]; E = T + electron_mass_c2;
334 G4double b2big = T*(E+electron_mass_c2)/(E*E);
335 G4double ratb2 = (beta2-b2small)/(b2big-b2small);
336
337 if (charge < 0.)
338 {
339 c1 = celectron[iZ][iT];
340 c2 = celectron[iZ+1][iT];
341 cc1 = c1+ratZ*(c2-c1);
342
343 c1 = celectron[iZ][iT+1];
344 c2 = celectron[iZ+1][iT+1];
345 cc2 = c1+ratZ*(c2-c1);
346
347 corr = cc1+ratb2*(cc2-cc1);
348
349 sigma *= sigmafactor/corr;
350 }
351 else
352 {
353 c1 = cpositron[iZ][iT];
354 c2 = cpositron[iZ+1][iT];
355 cc1 = c1+ratZ*(c2-c1);
356
357 c1 = cpositron[iZ][iT+1];
358 c2 = cpositron[iZ+1][iT+1];
359 cc2 = c1+ratZ*(c2-c1);
360
361 corr = cc1+ratb2*(cc2-cc1);
362
363 sigma *= sigmafactor/corr;
364 }
365 }
366 else
367 {
368 c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
369 c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
370 if((AtomicNumber >= Z1) && (AtomicNumber <= Z2))
371 sigma = c1+ratZ*(c2-c1) ;
372 else if(AtomicNumber < Z1)
373 sigma = AtomicNumber*AtomicNumber*c1/(Z1*Z1);
374 else if(AtomicNumber > Z2)
375 sigma = AtomicNumber*AtomicNumber*c2/(Z2*Z2);
376 }
377 return sigma;
378
379}
380
381//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
382
384{
385 SetParticle(track->GetDynamicParticle()->GetDefinition());
386 firstStep = true;
387 inside = false;
388 insideskin = false;
389 tlimit = geombig;
390}
391
392//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
393
395 const G4Track& track,
396 G4double& currentMinimalStep)
397{
398 tPathLength = currentMinimalStep;
399 const G4DynamicParticle* dp = track.GetDynamicParticle();
400 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
401 G4StepStatus stepStatus = sp->GetStepStatus();
402 couple = track.GetMaterialCutsCouple();
403 SetCurrentCouple(couple);
404 currentMaterialIndex = couple->GetIndex();
405 currentKinEnergy = dp->GetKineticEnergy();
406 currentRange = GetRange(particle,currentKinEnergy,couple);
407 lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy);
408
409 // stop here if small range particle
410 if(inside || tPathLength < tlimitminfix) {
411 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
412 }
413
414 if(tPathLength > currentRange) { tPathLength = currentRange; }
415
416 presafety = sp->GetSafety();
417 /*
418 G4cout << "G4UrbanMscModel90::ComputeTruePathLengthLimit tPathLength= "
419 <<tPathLength<<" safety= " << presafety
420 << " range= " <<currentRange<<G4endl;
421 */
422 // far from geometry boundary
423 if(currentRange < presafety)
424 {
425 inside = true;
426 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
427 }
428
429 // standard version
430 //
432 {
433 //compute geomlimit and presafety
434 geomlimit = ComputeGeomLimit(track, presafety, currentRange);
435
436 // is far from boundary
437 if(currentRange <= presafety)
438 {
439 inside = true;
440 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
441 }
442
443 smallstep += 1.;
444 insideskin = false;
445
446 if(firstStep || stepStatus == fGeomBoundary)
447 {
448
449 if(firstStep) smallstep = 1.e10;
450 else smallstep = 1.;
451
452 // facrange scaling in lambda
453 // not so strong step restriction above lambdalimit
454 G4double facr = facrange;
455 if(lambda0 > lambdalimit)
456 facr *= frscaling1+frscaling2*lambda0/lambdalimit;
457
458 // constraint from the physics
459 if (currentRange > lambda0) tlimit = facr*currentRange;
460 else tlimit = facr*lambda0;
461
462 // constraint from the geometry (if tlimit above is too big)
463 G4double tgeom = geombig;
464 if(geomlimit > geommin)
465 {
466 if(stepStatus == fGeomBoundary)
467 tgeom = geomlimit/facgeom;
468 else
469 tgeom = 2.*geomlimit/facgeom;
470 }
471
472 //define stepmin here (it depends on lambda!)
473 //rough estimation of lambda_elastic/lambda_transport
474 G4double rat = currentKinEnergy/MeV ;
475 rat = 1.e-3/(rat*(10.+rat)) ;
476 //stepmin ~ lambda_elastic
477 stepmin = rat*lambda0;
478 skindepth = skin*stepmin;
479
480 //define tlimitmin
481 tlimitmin = lambda0/nstepmax;
482 if(tlimitmin < stepmin) tlimitmin = 1.01*stepmin;
483 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
484
485 //lower limit for tlimit
486 if(tlimit < tlimitmin) tlimit = tlimitmin;
487
488 //check against geometry limit
489 if(tlimit > tgeom) tlimit = tgeom;
490 }
491
492 //if track starts far from boundaries increase tlimit!
493 if(tlimit < facsafety*presafety) tlimit = facsafety*presafety ;
494
495 // G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
496 // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
497
498 // shortcut
499 if((tPathLength < tlimit) && (tPathLength < presafety))
500 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
501
502 G4double tnow = tlimit;
503 // optimization ...
504 if(geomlimit < geombig) tnow = max(tlimit,facsafety*geomlimit);
505
506 // step reduction near to boundary
507 if(smallstep < skin)
508 {
509 tnow = stepmin;
510 insideskin = true;
511 }
512 else if(geomlimit < geombig)
513 {
514 if(geomlimit > skindepth)
515 {
516 if(tnow > geomlimit-0.999*skindepth)
517 tnow = geomlimit-0.999*skindepth;
518 }
519 else
520 {
521 insideskin = true;
522 if(tnow > stepmin) tnow = stepmin;
523 }
524 }
525
526 if(tnow < stepmin) tnow = stepmin;
527
528 if(tPathLength > tnow) tPathLength = tnow ;
529 }
530 // for 'normal' simulation with or without magnetic field
531 // there no small step/single scattering at boundaries
532 else if(steppingAlgorithm == fUseSafety)
533 {
534 // compute presafety again if presafety <= 0 and no boundary
535 // i.e. when it is needed for optimization purposes
536 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
537 presafety = ComputeSafety(sp->GetPosition(),tPathLength);
538
539 // is far from boundary
540 if(currentRange < presafety)
541 {
542 inside = true;
543 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
544 }
545
546 if(firstStep || stepStatus == fGeomBoundary)
547 {
548 // facrange scaling in lambda
549 // not so strong step restriction above lambdalimit
550 G4double facr = facrange;
551 if(lambda0 > lambdalimit)
552 facr *= frscaling1+frscaling2*lambda0/lambdalimit;
553
554 // constraint from the physics
555 if (currentRange > lambda0) tlimit = facr*currentRange;
556 else tlimit = facr*lambda0;
557
558 //lower limit for tlimit
559 tlimitmin = std::max(tlimitminfix,lambda0/nstepmax);
560 if(tlimit < tlimitmin) tlimit = tlimitmin;
561 }
562
563 //if track starts far from boundaries increase tlimit!
564 if(tlimit < facsafety*presafety) tlimit = facsafety*presafety ;
565
566 if(tPathLength > tlimit) tPathLength = tlimit;
567 }
568
569 // version similar to 7.1 (needed for some experiments)
570 else
571 {
572 if (stepStatus == fGeomBoundary)
573 {
574 if (currentRange > lambda0) tlimit = facrange*currentRange;
575 else tlimit = facrange*lambda0;
576
577 if(tlimit < tlimitmin) tlimit = tlimitmin;
578 if(tPathLength > tlimit) tPathLength = tlimit;
579 }
580 }
581 // G4cout << "tPathLength= " << tPathLength << " geomlimit= " << geomlimit
582 // << " currentMinimalStep= " << currentMinimalStep << G4endl;
583 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
584}
585
586//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
587
589{
590 firstStep = false;
591 lambdaeff = lambda0;
592 par1 = -1. ;
593 par2 = par3 = 0. ;
594
595 // do the true -> geom transformation
596 zPathLength = tPathLength;
597
598 // z = t for very small tPathLength
599 if(tPathLength < tlimitminfix) return zPathLength;
600
601 // this correction needed to run MSC with eIoni and eBrem inactivated
602 // and makes no harm for a normal run
603 if(tPathLength > currentRange)
604 tPathLength = currentRange ;
605
606 G4double tau = tPathLength/lambda0 ;
607
608 if ((tau <= tausmall) || insideskin) {
609 zPathLength = tPathLength;
610 if(zPathLength > lambda0) zPathLength = lambda0;
611 return zPathLength;
612 }
613
614 G4double zmean = tPathLength;
615 if (tPathLength < currentRange*dtrl) {
616 if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
617 else zmean = lambda0*(1.-exp(-tau));
618 } else if(currentKinEnergy < mass || tPathLength == currentRange) {
619 par1 = 1./currentRange ;
620 par2 = 1./(par1*lambda0) ;
621 par3 = 1.+par2 ;
622 if(tPathLength < currentRange)
623 zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
624 else
625 zmean = 1./(par1*par3) ;
626 } else {
627 G4double T1 = GetEnergy(particle,currentRange-tPathLength,couple);
628 G4double lambda1 = GetTransportMeanFreePath(particle,T1);
629
630 par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
631 par2 = 1./(par1*lambda0) ;
632 par3 = 1.+par2 ;
633 zmean = (1.-exp(par3*log(lambda1/lambda0)))/(par1*par3) ;
634 }
635
636 zPathLength = zmean ;
637
638 // sample z
639 if(samplez)
640 {
641 const G4double ztmax = 0.99, onethird = 1./3. ;
642 G4double zt = zmean/tPathLength ;
643
644 if (tPathLength > stepmin && zt < ztmax)
645 {
646 G4double u,cz1;
647 if(zt >= onethird)
648 {
649 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
650 cz1 = 1.+cz ;
651 G4double u0 = cz/cz1 ;
652 G4double grej ;
653 do {
654 u = exp(log(G4UniformRand())/cz1) ;
655 grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
656 } while (grej < G4UniformRand()) ;
657 }
658 else
659 {
660 cz1 = 1./zt-1.;
661 u = 1.-exp(log(G4UniformRand())/cz1) ;
662 }
663 zPathLength = tPathLength*u ;
664 }
665 }
666
667 if(zPathLength > lambda0) zPathLength = lambda0;
668
669 return zPathLength;
670}
671
672//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
673
675{
676 // step defined other than transportation
677 if(geomStepLength == zPathLength && tPathLength <= currentRange)
678 { return tPathLength; }
679
680 // t = z for very small step
681 zPathLength = geomStepLength;
682 tPathLength = geomStepLength;
683 if(geomStepLength < tlimitminfix) return tPathLength;
684
685 // recalculation
686 if((geomStepLength > lambda0*tausmall) && !insideskin)
687 {
688 if(par1 < 0.)
689 tPathLength = -lambda0*log(1.-geomStepLength/lambda0) ;
690 else
691 {
692 if(par1*par3*geomStepLength < 1.)
693 tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
694 else
695 tPathLength = currentRange;
696 }
697 }
698 if(tPathLength < geomStepLength) tPathLength = geomStepLength;
699
700 return tPathLength;
701}
702
703//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
704
706 G4double KineticEnergy)
707{
708 // for all particles take the width of the central part
709 // from a parametrization similar to the Highland formula
710 // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
711 const G4double c_highland = 13.6*MeV ;
712 G4double betacp = sqrt(currentKinEnergy*(currentKinEnergy+2.*mass)*
713 KineticEnergy*(KineticEnergy+2.*mass)/
714 ((currentKinEnergy+mass)*(KineticEnergy+mass)));
715 G4double y = trueStepLength/currentRadLength;
716 G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp;
717 y = log(y);
718 theta0 *= sqrt(1.+y*(0.105+0.0035*y));
719
720 //correction for small Zeff (based on high energy
721 // proton scattering data)
722 // see G.Shen at al. Phys.Rev.D20(1979) p.1584
723 theta0 *= 1.-0.24/(Zeff*(Zeff+1.));
724
725 return theta0;
726
727}
728
729//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
730
733 G4double safety)
734{
735 fDisplacement.set(0.0,0.0,0.0);
736 G4double kineticEnergy = currentKinEnergy;
737 if (tPathLength > currentRange*dtrl) {
738 kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
739 } else {
740 kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple);
741 }
742 if((kineticEnergy <= eV) || (tPathLength <= tlimitminfix) ||
743 (tPathLength/tausmall < lambda0) ) { return fDisplacement; }
744
745 G4double cth = SampleCosineTheta(tPathLength,kineticEnergy);
746 // protection against 'bad' cth values
747 if(std::fabs(cth) > 1.) { return fDisplacement; }
748
749 const G4double checkEnergy = GeV;
750 if(kineticEnergy > checkEnergy && cth < 0.0
751 && tPathLength < taulim*lambda0) {
753 ed << dynParticle->GetDefinition()->GetParticleName()
754 << " E(MeV)= " << kineticEnergy/MeV
755 << " Step(mm)= " << tPathLength/mm
756 << " in " << CurrentCouple()->GetMaterial()->GetName()
757 << " CosTheta= " << cth
758 << " is too big - the angle is set to zero" << G4endl;
759 G4Exception("G4UrbanMscModel90::SampleScattering","em0004",JustWarning,
760 ed,"");
761 return fDisplacement;
762 }
763
764 G4double sth = sqrt((1.0 - cth)*(1.0 + cth));
765 G4double phi = twopi*G4UniformRand();
766 G4double dirx = sth*cos(phi);
767 G4double diry = sth*sin(phi);
768
769 G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
770 G4ThreeVector newDirection(dirx,diry,cth);
771 newDirection.rotateUz(oldDirection);
772 fParticleChange->ProposeMomentumDirection(newDirection);
773
774 if (latDisplasment && safety > tlimitminfix) {
775
776 G4double r = SampleDisplacement();
777/*
778 G4cout << "G4UrbanMscModel90::SampleSecondaries: e(MeV)= " << kineticEnergy
779 << " sinTheta= " << sth << " r(mm)= " << r
780 << " trueStep(mm)= " << truestep
781 << " geomStep(mm)= " << zPathLength
782 << G4endl;
783*/
784 if(r > 0.)
785 {
786 G4double latcorr = LatCorrelation();
787 if(latcorr > r) latcorr = r;
788
789 // sample direction of lateral displacement
790 // compute it from the lateral correlation
791 G4double Phi = 0.;
792 if(std::abs(r*sth) < latcorr) {
793 Phi = twopi*G4UniformRand();
794 } else {
795 G4double psi = std::acos(latcorr/(r*sth));
796 if(G4UniformRand() < 0.5) Phi = phi+psi;
797 else Phi = phi-psi;
798 }
799
800 dirx = std::cos(Phi);
801 diry = std::sin(Phi);
802
803 fDisplacement.set(r*dirx,r*diry,0.0);
804 fDisplacement.rotateUz(oldDirection);
805 }
806 }
807 return fDisplacement;
808}
809
810//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
811
812G4double G4UrbanMscModel90::SampleCosineTheta(G4double trueStepLength,
813 G4double KineticEnergy)
814{
815 G4double cth = 1. ;
816 G4double tau = trueStepLength/lambda0 ;
817
818 Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
820
821 if(insideskin)
822 {
823 //no scattering, single or plural scattering
824 G4double mean = trueStepLength/stepmin ;
825
826 G4int n = G4Poisson(mean);
827 if(n > 0)
828 {
829 G4double tm = KineticEnergy/electron_mass_c2;
830 // ascr - screening parameter
831 G4double ascr = exp(log(Zeff)/3.)/(137.*sqrt(tm*(tm+2.)));
832 G4double ascr1 = 1.+0.5*ascr*ascr;
833 G4double bp1=ascr1+1.;
834 G4double bm1=ascr1-1.;
835 // single scattering from screened Rutherford x-section
836 G4double ct,st,phi;
837 G4double sx=0.,sy=0.,sz=0.;
838 for(G4int i=1; i<=n; i++)
839 {
840 ct = ascr1-bp1*bm1/(2.*G4UniformRand()+bm1);
841 if(ct < -1.) ct = -1.;
842 if(ct > 1.) ct = 1.;
843 st = sqrt(1.-ct*ct);
844 phi = twopi*G4UniformRand();
845 sx += st*cos(phi);
846 sy += st*sin(phi);
847 sz += ct;
848 }
849 cth = sz/sqrt(sx*sx+sy*sy+sz*sz);
850 }
851 }
852 else
853 {
854 if(trueStepLength >= currentRange*dtrl) {
855 if(par1*trueStepLength < 1.)
856 tau = -par2*log(1.-par1*trueStepLength) ;
857 // for the case if ioni/brems are inactivated
858 // see the corresponding condition in ComputeGeomPathLength
859 else if(1.-KineticEnergy/currentKinEnergy > taulim)
860 tau = taubig ;
861 }
862 currentTau = tau ;
863 lambdaeff = trueStepLength/currentTau;
864 currentRadLength = couple->GetMaterial()->GetRadlen();
865
866 if (tau >= taubig) cth = -1.+2.*G4UniformRand();
867 else if (tau >= tausmall)
868 {
869 G4double b,bx,b1,ebx,eb1;
870 G4double prob = 0., qprob = 1. ;
871 G4double a = 1., ea = 0., eaa = 1.;
872 G4double xmean1 = 1., xmean2 = 0.;
873 G4double xsi = 3.;
874
875 G4double theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
876
877 // protexction for very small angles
878 if(theta0*theta0 < tausmall) return cth;
879
880 G4double sth = sin(0.5*theta0);
881 a = 0.25/(sth*sth);
882
883 G4double xmeanth = exp(-tau);
884
885 G4double c = 3. ;
886 G4double c1 = c-1.;
887
888 G4double x0 = 1.-xsi/a ;
889 if(x0 < 0.)
890 {
891 // 1 model function
892 b = exp(tau);
893 bx = b-1.;
894 b1 = b+1.;
895 ebx=exp((c1)*log(bx)) ;
896 eb1=exp((c1)*log(b1)) ;
897 }
898 else
899 {
900 //empirical tail parameter
901 // based some exp. data
902 c = 2.40-0.027*exp(2.*log(Zeff)/3.);
903
904 if(c == 2.) c = 2.+taulim ;
905 if(c <= 1.) c = 1.+taulim ;
906 c1 = c-1.;
907
908 ea = exp(-xsi) ;
909 eaa = 1.-ea ;
910 xmean1 = 1.-(1.-(1.+xsi)*ea)/(eaa*a) ;
911
912 // from the continuity of the 1st derivative at x=x0
913 b = 1.+(c-xsi)/a ;
914
915 b1 = b+1. ;
916 bx = c/a ;
917 eb1=exp((c1)*log(b1)) ;
918 ebx=exp((c1)*log(bx)) ;
919 xmean2 = (x0*eb1+ebx-(eb1*bx-b1*ebx)/(c-2.))/(eb1-ebx) ;
920
921 G4double f1x0 = a*ea/eaa ;
922 G4double f2x0 = c1*eb1*ebx/(eb1-ebx)/exp(c*log(bx)) ;
923
924 // from continuity at x=x0
925 prob = f2x0/(f1x0+f2x0) ;
926
927 // from xmean = xmeanth
928 qprob = (f1x0+f2x0)*xmeanth/(f2x0*xmean1+f1x0*xmean2) ;
929 }
930
931 // sampling of costheta
932 if (G4UniformRand() < qprob)
933 {
934 if (G4UniformRand() < prob)
935 cth = 1.+log(ea+G4UniformRand()*eaa)/a ;
936 else
937 cth = b-b1*bx/exp(log(ebx-G4UniformRand()*(ebx-eb1))/c1) ;
938 }
939 else
940 {
941 cth = -1.+2.*G4UniformRand();
942 }
943 }
944 }
945
946 return cth ;
947}
948//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
949
950G4double G4UrbanMscModel90::SampleDisplacement()
951{
952 const G4double kappa = 2.5;
953 const G4double kappapl1 = kappa+1.;
954 const G4double kappami1 = kappa-1.;
955 G4double rmean = 0.0;
956 if ((currentTau >= tausmall) && !insideskin) {
957 if (currentTau < taulim) {
958 rmean = kappa*currentTau*currentTau*currentTau*
959 (1.-kappapl1*currentTau*0.25)/6. ;
960
961 } else {
962 G4double etau = 0.0;
963 if (currentTau<taubig) etau = exp(-currentTau);
964 rmean = -kappa*currentTau;
965 rmean = -exp(rmean)/(kappa*kappami1);
966 rmean += currentTau-kappapl1/kappa+kappa*etau/kappami1;
967 }
968 if (rmean>0.) rmean = 2.*lambdaeff*sqrt(rmean/3.0);
969 else rmean = 0.;
970 }
971
972 // protection against z > t ...........................
973 if(rmean > 0.) {
974 G4double zt = (tPathLength-zPathLength)*(tPathLength+zPathLength);
975 if(zt <= 0.)
976 rmean = 0.;
977 else if(rmean*rmean > zt)
978 rmean = sqrt(zt);
979 }
980 return rmean;
981}
982
983//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
984
985G4double G4UrbanMscModel90::LatCorrelation()
986{
987 const G4double kappa = 2.5;
988 const G4double kappami1 = kappa-1.;
989
990 G4double latcorr = 0.;
991 if((currentTau >= tausmall) && !insideskin)
992 {
993 if(currentTau < taulim)
994 latcorr = lambdaeff*kappa*currentTau*currentTau*
995 (1.-(kappa+1.)*currentTau/3.)/3.;
996 else
997 {
998 G4double etau = 0.;
999 if(currentTau < taubig) etau = exp(-currentTau);
1000 latcorr = -kappa*currentTau;
1001 latcorr = exp(latcorr)/kappami1;
1002 latcorr += 1.-kappa*etau/kappami1 ;
1003 latcorr *= 2.*lambdaeff/3. ;
1004 }
1005 }
1006
1007 return latcorr;
1008}
1009
1010//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1011
1012void G4UrbanMscModel90::SampleSecondaries(std::vector<G4DynamicParticle*>*,
1013 const G4MaterialCutsCouple*,
1014 const G4DynamicParticle*,
1015 G4double,
1016 G4double)
1017{}
1018
1019//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ JustWarning
@ fUseSafety
@ fUseDistanceToBoundary
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:50
G4StepStatus
Definition: G4StepStatus.hh:51
@ fGeomBoundary
Definition: G4StepStatus.hh:54
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
void set(double x, double y, double z)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4LossTableManager * Instance()
const G4Material * GetMaterial() const
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:208
G4double GetTotNbOfElectPerVolume() const
Definition: G4Material.hh:211
G4double GetRadlen() const
Definition: G4Material.hh:219
const G4String & GetName() const
Definition: G4Material.hh:177
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
const G4String & GetParticleName() const
G4StepPoint * GetPreStepPoint() const
const G4DynamicParticle * GetDynamicParticle() const
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
const G4Step * GetStep() const
G4double ComputeTheta0(G4double truePathLength, G4double KineticEnergy)
G4double ComputeGeomPathLength(G4double truePathLength)
void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4UrbanMscModel90(const G4String &nam="UrbanMsc90")
G4ThreeVector & SampleScattering(const G4DynamicParticle *, G4double safety)
G4double ComputeTrueStepLength(G4double geomStepLength)
G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep)
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double)
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=0., G4double emax=DBL_MAX)
void StartTracking(G4Track *)
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:370
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:377
G4double dtrl
Definition: G4VMscModel.hh:180
G4double GetDEDX(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:273
G4double facrange
Definition: G4VMscModel.hh:176
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:256
G4bool samplez
Definition: G4VMscModel.hh:188
G4double skin
Definition: G4VMscModel.hh:179
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:332
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:304
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:288
G4double lambdalimit
Definition: G4VMscModel.hh:181
G4double ComputeSafety(const G4ThreeVector &position, G4double limit)
Definition: G4VMscModel.hh:238
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:186
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:89
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:246
G4bool latDisplasment
Definition: G4VMscModel.hh:189
G4double facsafety
Definition: G4VMscModel.hh:178
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:185
G4double facgeom
Definition: G4VMscModel.hh:177
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76