Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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: G4UrbanMscModel90.cc 66592 2012-12-23 09:34:55Z vnivanch $
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"
56 #include "G4PhysicalConstants.hh"
57 #include "G4SystemOfUnits.hh"
58 #include "Randomize.hh"
59 #include "G4Electron.hh"
60 
61 #include "G4LossTableManager.hh"
62 #include "G4ParticleChangeForMSC.hh"
63 
64 #include "G4Poisson.hh"
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 
68 using namespace std;
69 
70 G4UrbanMscModel90::G4UrbanMscModel90(const G4String& nam)
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 
113 G4UrbanMscModel90::~G4UrbanMscModel90()
114 {}
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
118 void G4UrbanMscModel90::Initialise(const G4ParticleDefinition* p,
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 
139 G4double G4UrbanMscModel90::ComputeCrossSectionPerAtom(
140  const G4ParticleDefinition* part,
141  G4double KineticEnergy,
142  G4double AtomicNumber,G4double,
143  G4double, 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 
383 void G4UrbanMscModel90::StartTracking(G4Track* track)
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 
394 G4double G4UrbanMscModel90::ComputeTruePathLengthLimit(
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  //
431  if (steppingAlgorithm == fUseDistanceToBoundary)
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 
588 G4double G4UrbanMscModel90::ComputeGeomPathLength(G4double)
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 
674 G4double G4UrbanMscModel90::ComputeTrueStepLength(G4double geomStepLength)
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 
705 G4double G4UrbanMscModel90::ComputeTheta0(G4double trueStepLength,
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 
731 G4ThreeVector&
732 G4UrbanMscModel90::SampleScattering(const G4ThreeVector& oldDirection,
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) {
752  G4ExceptionDescription ed;
753  ed << particle->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 
812 G4double G4UrbanMscModel90::SampleCosineTheta(G4double trueStepLength,
813  G4double KineticEnergy)
814 {
815  G4double cth = 1. ;
816  G4double tau = trueStepLength/lambda0 ;
817 
818  Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
819  couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ;
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 
950 G4double 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 
985 G4double 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 
1012 void G4UrbanMscModel90::SampleSecondaries(std::vector<G4DynamicParticle*>*,
1013  const G4MaterialCutsCouple*,
1014  const G4DynamicParticle*,
1015  G4double,
1016  G4double)
1017 {}
1018 
1019 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:33:38 for Geant4 by   doxygen 1.8.4