CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
UrbanMscModel93 Class Reference

#include <UrbanMscModel93.h>

Inheritance diagram for UrbanMscModel93:

Public Member Functions

G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=0., G4double emax=DBL_MAX)
 
G4double ComputeGeomPathLength (G4double truePathLength)
 
G4double ComputeTheta0 (G4double truePathLength, G4double KineticEnergy)
 
G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep)
 
G4double ComputeTrueStepLength (G4double geomStepLength)
 
void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
G4ThreeVector & SampleScattering (const G4ThreeVector &, G4double safety)
 
void StartTracking (G4Track *)
 
 UrbanMscModel93 (const G4String &nam="UrbanMsc93")
 
virtual ~UrbanMscModel93 ()
 

Private Member Functions

G4double LatCorrelation ()
 
UrbanMscModel93operator= (const UrbanMscModel93 &right)
 
G4double SampleCosineTheta (G4double trueStepLength, G4double KineticEnergy)
 
G4double SampleDisplacement ()
 
void SetParticle (const G4ParticleDefinition *)
 
G4double SimpleScattering (G4double xmeanth, G4double x2meanth)
 
void UpdateCache ()
 
 UrbanMscModel93 (const UrbanMscModel93 &)
 

Private Attributes

G4double charge
 
G4double ChargeSquare
 
G4double coeffc1
 
G4double coeffc2
 
G4double coeffth1
 
G4double coeffth2
 
const G4MaterialCutsCouple * couple
 
G4double currentKinEnergy
 
G4int currentMaterialIndex
 
G4double currentRadLength
 
G4double currentRange
 
G4double currentTau
 
G4double ea
 
G4double eaa
 
G4bool firstStep
 
G4ParticleChangeForMSC * fParticleChange
 
G4double fr
 
G4double geombig
 
G4double geomlimit
 
G4double geommin
 
G4bool inside
 
G4bool insideskin
 
G4double lambda0
 
G4double lambdaeff
 
G4double lambdalimit
 
G4double lnZ
 
G4double mass
 
G4double masslimite
 
G4double numlim
 
G4double par1
 
G4double par2
 
G4double par3
 
const G4ParticleDefinition * particle
 
G4double presafety
 
G4double rangeinit
 
G4double rellossmax
 
G4double scr1
 
G4double scr1ini
 
G4double scr2
 
G4double scr2ini
 
G4double skindepth
 
G4double smallstep
 
G4double stepmin
 
G4double taubig
 
G4double taulim
 
G4double tausmall
 
G4double tgeom
 
G4LossTableManager * theManager
 
G4double theta0max
 
G4double third
 
G4double tlimit
 
G4double tlimitmin
 
G4double tlimitminfix
 
G4double tPathLength
 
G4double xsi
 
G4double y
 
G4double Z2
 
G4double Z23
 
G4double Zeff
 
G4double Zold
 
G4double zPathLength
 

Detailed Description

Definition at line 69 of file UrbanMscModel93.h.

Constructor & Destructor Documentation

UrbanMscModel93::UrbanMscModel93 ( const G4String &  nam = "UrbanMsc93")

Definition at line 76 of file UrbanMscModel93.cc.

References charge, ChargeSquare, coeffc1, coeffc2, coeffth1, coeffth2, couple, currentKinEnergy, currentMaterialIndex, currentRadLength, currentRange, currentTau, ea, eaa, firstStep, fParticleChange, fr, geombig, geomlimit, geommin, inside, insideskin, lambda0, lambdaeff, lambdalimit, lnZ, mass, masslimite, MeV, numlim, par1, par2, par3, particle, pi, presafety, rangeinit, rellossmax, scr1, scr1ini, scr2, scr2ini, skindepth, smallstep, stepmin, taubig, taulim, tausmall, tgeom, theManager, theta0max, third, tlimit, tlimitmin, tlimitminfix, tPathLength, xsi, y, Z2, Z23, Zeff, Zold, and zPathLength.

77  : G4VMscModel(nam)
78 {
79  masslimite = 0.6*MeV;
80  lambdalimit = 1.*mm;
81  fr = 0.02;
82  taubig = 8.0;
83  tausmall = 1.e-16;
84  taulim = 1.e-6;
86  tlimitminfix = 1.e-6*mm;
88  smallstep = 1.e10;
89  currentRange = 0. ;
90  rangeinit = 0.;
91  tlimit = 1.e10*mm;
92  tlimitmin = 10.*tlimitminfix;
93  tgeom = 1.e50*mm;
94  geombig = 1.e50*mm;
95  geommin = 1.e-3*mm;
97  presafety = 0.*mm;
98 
99  y = 0.;
100 
101  Zold = 0.;
102  Zeff = 1.;
103  Z2 = 1.;
104  Z23 = 1.;
105  lnZ = 0.;
106  coeffth1 = 0.;
107  coeffth2 = 0.;
108  coeffc1 = 0.;
109  coeffc2 = 0.;
110  scr1ini = fine_structure_const*fine_structure_const*
111  electron_mass_c2*electron_mass_c2/(0.885*0.885*4.*pi);
112  scr2ini = 3.76*fine_structure_const*fine_structure_const;
113  scr1 = 0.;
114  scr2 = 0.;
115 
116  theta0max = pi/6.;
117  rellossmax = 0.50;
118  third = 1./3.;
119  particle = 0;
120  theManager = G4LossTableManager::Instance();
121  firstStep = true;
122  inside = false;
123  insideskin = false;
124 
125  numlim = 0.01;
126  xsi = 3.;
127  ea = G4Exp(-xsi);
128  eaa = 1.-ea ;
129 
130  skindepth = skin*stepmin;
131 
132  mass = proton_mass_c2;
133  charge = ChargeSquare = 1.0;
135  = zPathLength = par1 = par2 = par3 = 0;
136 
138  fParticleChange = 0;
139  couple = 0;
140  SetSampleZ(false);
141 }
const G4ParticleDefinition * particle
const G4MaterialCutsCouple * couple
G4ParticleChangeForMSC * fParticleChange
G4double currentRadLength
const Double_t pi
const double MeV
G4LossTableManager * theManager
G4double currentKinEnergy
UrbanMscModel93::~UrbanMscModel93 ( )
virtual

Definition at line 145 of file UrbanMscModel93.cc.

146 {}
UrbanMscModel93::UrbanMscModel93 ( const UrbanMscModel93 )
private

Member Function Documentation

G4double UrbanMscModel93::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  particle,
G4double  KineticEnergy,
G4double  AtomicNumber,
G4double  AtomicWeight = 0.,
G4double  cut = 0.,
G4double  emax = DBL_MAX 
)

Definition at line 170 of file UrbanMscModel93.cc.

References EnergyCorrector::c, alignmentValidation::c1, charge, ChargeSquare, corr, hbarc, mass, MeV, funct::pow(), SetParticle(), mathSSE::sqrt(), pat::TAU, metsig::tau, w, and Z23.

175 {
176  static const G4double sigmafactor =
177  twopi*classic_electr_radius*classic_electr_radius;
178  static const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2*
179  Bohr_radius*Bohr_radius/(hbarc*hbarc);
180  static const G4double epsmin = 1.e-4 , epsmax = 1.e10;
181 
182  static const G4double Zdat[15] = { 4., 6., 13., 20., 26., 29., 32., 38., 47.,
183  50., 56., 64., 74., 79., 82. };
184 
185  static const G4double Tdat[22] = { 100*eV, 200*eV, 400*eV, 700*eV,
186  1*keV, 2*keV, 4*keV, 7*keV,
187  10*keV, 20*keV, 40*keV, 70*keV,
188  100*keV, 200*keV, 400*keV, 700*keV,
189  1*MeV, 2*MeV, 4*MeV, 7*MeV,
190  10*MeV, 20*MeV};
191 
192  // corr. factors for e-/e+ lambda for T <= Tlim
193  static const G4double celectron[15][22] =
194  {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
195  1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
196  1.112,1.108,1.100,1.093,1.089,1.087 },
197  {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
198  1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
199  1.109,1.105,1.097,1.090,1.086,1.082 },
200  {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
201  1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
202  1.131,1.124,1.113,1.104,1.099,1.098 },
203  {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
204  1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
205  1.112,1.105,1.096,1.089,1.085,1.098 },
206  {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
207  1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
208  1.073,1.070,1.064,1.059,1.056,1.056 },
209  {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
210  1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
211  1.074,1.070,1.063,1.059,1.056,1.052 },
212  {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
213  1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
214  1.068,1.064,1.059,1.054,1.051,1.050 },
215  {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
216  1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
217  1.039,1.037,1.034,1.031,1.030,1.036 },
218  {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
219  1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
220  1.031,1.028,1.024,1.022,1.021,1.024 },
221  {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
222  1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
223  1.020,1.017,1.015,1.013,1.013,1.020 },
224  {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
225  1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
226  0.995,0.993,0.993,0.993,0.993,1.011 },
227  {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
228  1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
229  0.974,0.972,0.973,0.974,0.975,0.987 },
230  {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
231  1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
232  0.950,0.947,0.949,0.952,0.954,0.963 },
233  {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
234  1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
235  0.941,0.938,0.940,0.944,0.946,0.954 },
236  {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
237  1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
238  0.933,0.930,0.933,0.936,0.939,0.949 }};
239 
240  static const G4double cpositron[15][22] = {
241  {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
242  1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
243  1.131,1.126,1.117,1.108,1.103,1.100 },
244  {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
245  1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
246  1.138,1.132,1.122,1.113,1.108,1.102 },
247  {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
248  1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
249  1.203,1.190,1.173,1.159,1.151,1.145 },
250  {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
251  1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
252  1.225,1.210,1.191,1.175,1.166,1.174 },
253  {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
254  1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
255  1.217,1.203,1.184,1.169,1.160,1.151 },
256  {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
257  1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
258  1.237,1.222,1.201,1.184,1.174,1.159 },
259  {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
260  1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
261  1.252,1.234,1.212,1.194,1.183,1.170 },
262  {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
263  2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
264  1.254,1.237,1.214,1.195,1.185,1.179 },
265  {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
266  2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
267  1.312,1.288,1.258,1.235,1.221,1.205 },
268  {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
269  2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
270  1.320,1.294,1.264,1.240,1.226,1.214 },
271  {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
272  2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
273  1.328,1.302,1.270,1.245,1.231,1.233 },
274  {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
275  2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
276  1.361,1.330,1.294,1.267,1.251,1.239 },
277  {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
278  3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
279  1.409,1.372,1.330,1.298,1.280,1.258 },
280  {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
281  3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
282  1.442,1.400,1.354,1.319,1.299,1.272 },
283  {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
284  3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
285  1.456,1.412,1.364,1.328,1.307,1.282 }};
286 
287  //data/corrections for T > Tlim
288  static const G4double Tlim = 10.*MeV;
289  static const G4double beta2lim = Tlim*(Tlim+2.*electron_mass_c2)/
290  ((Tlim+electron_mass_c2)*(Tlim+electron_mass_c2));
291  static const G4double bg2lim = Tlim*(Tlim+2.*electron_mass_c2)/
292  (electron_mass_c2*electron_mass_c2);
293 
294  static const G4double sig0[15] = {
295  0.2672*barn, 0.5922*barn, 2.653*barn, 6.235*barn,
296  11.69*barn , 13.24*barn , 16.12*barn, 23.00*barn ,
297  35.13*barn , 39.95*barn , 50.85*barn, 67.19*barn ,
298  91.15*barn , 104.4*barn , 113.1*barn};
299 
300  static const G4double hecorr[15] = {
301  120.70, 117.50, 105.00, 92.92, 79.23, 74.510, 68.29,
302  57.39, 41.97, 36.14, 24.53, 10.21, -7.855, -16.84,
303  -22.30};
304 
305  G4double sigma;
306  SetParticle(part);
307 
308  Z23 = pow(AtomicNumber,2./3.);
309 
310  // correction if particle .ne. e-/e+
311  // compute equivalent kinetic energy
312  // lambda depends on p*beta ....
313 
314  G4double eKineticEnergy = KineticEnergy;
315 
316  if(mass > electron_mass_c2)
317  {
318  G4double TAU = KineticEnergy/mass ;
319  G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
320  G4double w = c-2. ;
321  G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
322  eKineticEnergy = electron_mass_c2*tau ;
323  }
324 
325  G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
326  G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
327  /(eTotalEnergy*eTotalEnergy);
328  G4double bg2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
329  /(electron_mass_c2*electron_mass_c2);
330 
331  G4double eps = epsfactor*bg2/Z23;
332 
333  if (eps<epsmin) sigma = 2.*eps*eps;
334  else if(eps<epsmax) sigma = G4Log(1.+2.*eps)-2.*eps/(1.+2.*eps);
335  else sigma = G4Log(2.*eps)-1.+1./eps;
336 
337  sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
338 
339  // interpolate in AtomicNumber and beta2
340  G4double c1,c2,cc1,cc2,corr;
341 
342  // get bin number in Z
343  G4int iZ = 14;
344  while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
345  if (iZ==14) iZ = 13;
346  if (iZ==-1) iZ = 0 ;
347 
348  G4double ZZ1 = Zdat[iZ];
349  G4double ZZ2 = Zdat[iZ+1];
350  G4double ratZ = (AtomicNumber-ZZ1)*(AtomicNumber+ZZ1)/
351  ((ZZ2-ZZ1)*(ZZ2+ZZ1));
352 
353  if(eKineticEnergy <= Tlim)
354  {
355  // get bin number in T (beta2)
356  G4int iT = 21;
357  while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
358  if(iT==21) iT = 20;
359  if(iT==-1) iT = 0 ;
360 
361  // calculate betasquare values
362  G4double T = Tdat[iT], E = T + electron_mass_c2;
363  G4double b2small = T*(E+electron_mass_c2)/(E*E);
364 
365  T = Tdat[iT+1]; E = T + electron_mass_c2;
366  G4double b2big = T*(E+electron_mass_c2)/(E*E);
367  G4double ratb2 = (beta2-b2small)/(b2big-b2small);
368 
369  if (charge < 0.)
370  {
371  c1 = celectron[iZ][iT];
372  c2 = celectron[iZ+1][iT];
373  cc1 = c1+ratZ*(c2-c1);
374 
375  c1 = celectron[iZ][iT+1];
376  c2 = celectron[iZ+1][iT+1];
377  cc2 = c1+ratZ*(c2-c1);
378 
379  corr = cc1+ratb2*(cc2-cc1);
380 
381  sigma *= sigmafactor/corr;
382  }
383  else
384  {
385  c1 = cpositron[iZ][iT];
386  c2 = cpositron[iZ+1][iT];
387  cc1 = c1+ratZ*(c2-c1);
388 
389  c1 = cpositron[iZ][iT+1];
390  c2 = cpositron[iZ+1][iT+1];
391  cc2 = c1+ratZ*(c2-c1);
392 
393  corr = cc1+ratb2*(cc2-cc1);
394 
395  sigma *= sigmafactor/corr;
396  }
397  }
398  else
399  {
400  c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
401  c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
402  if((AtomicNumber >= ZZ1) && (AtomicNumber <= ZZ2))
403  sigma = c1+ratZ*(c2-c1) ;
404  else if(AtomicNumber < ZZ1)
405  sigma = AtomicNumber*AtomicNumber*c1/(ZZ1*ZZ1);
406  else if(AtomicNumber > ZZ2)
407  sigma = AtomicNumber*AtomicNumber*c2/(ZZ2*ZZ2);
408  }
409  return sigma;
410 
411 }
const double hbarc
Definition: MathUtil.h:18
const double w
Definition: UKUtility.cc:23
void SetParticle(const G4ParticleDefinition *)
const double MeV
T sqrt(T t)
Definition: SSEVec.h:48
JetCorrectorParameters corr
Definition: classes.h:5
part
Definition: HCALResponse.h:20
long double T
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
G4double UrbanMscModel93::ComputeGeomPathLength ( G4double  truePathLength)

Definition at line 645 of file UrbanMscModel93.cc.

References couple, currentKinEnergy, currentRange, firstStep, insideskin, lambda0, lambdaeff, mass, par1, par2, par3, particle, stepmin, metsig::tau, taulim, tausmall, third, tlimitminfix, tPathLength, and zPathLength.

646 {
647  firstStep = false;
648  lambdaeff = lambda0;
649  par1 = -1. ;
650  par2 = par3 = 0. ;
651 
652  // do the true -> geom transformation
654 
655  // z = t for very small tPathLength
656  if(tPathLength < tlimitminfix) return zPathLength;
657 
658  // this correction needed to run MSC with eIoni and eBrem inactivated
659  // and makes no harm for a normal run
662 
663  G4double tau = tPathLength/lambda0 ;
664 
665  if ((tau <= tausmall) || insideskin) {
668  return zPathLength;
669  }
670 
671  G4double zmean = tPathLength;
672  if (tPathLength < currentRange*dtrl) {
673  if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
674  else zmean = lambda0*(1.-G4Exp(-tau));
675  } else if(currentKinEnergy < mass || tPathLength == currentRange) {
676  par1 = 1./currentRange ;
677  par2 = 1./(par1*lambda0) ;
678  par3 = 1.+par2 ;
680  zmean = (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3) ;
681  else
682  zmean = 1./(par1*par3) ;
683  } else {
684  G4double T1 = GetEnergy(particle,currentRange-tPathLength,couple);
685  G4double lambda1 = GetTransportMeanFreePath(particle,T1);
686 
687  par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
688  par2 = 1./(par1*lambda0) ;
689  par3 = 1.+par2 ;
690  zmean = (1.-G4Exp(par3*G4Log(lambda1/lambda0)))/(par1*par3) ;
691  }
692 
693  zPathLength = zmean ;
694 
695  // sample z
696  if(samplez)
697  {
698  const G4double ztmax = 0.99 ;
699  G4double zt = zmean/tPathLength ;
700 
701  if (tPathLength > stepmin && zt < ztmax)
702  {
703  G4double u,cz1;
704  if(zt >= third)
705  {
706  G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
707  cz1 = 1.+cz ;
708  G4double u0 = cz/cz1 ;
709  G4double grej ;
710  do {
711  u = G4Exp(G4Log(G4UniformRand())/cz1) ;
712  grej = G4Exp(cz*G4Log(u/u0))*(1.-u)/(1.-u0) ;
713  } while (grej < G4UniformRand()) ;
714  }
715  else
716  {
717  cz1 = 1./zt-1.;
718  u = 1.-G4Exp(G4Log(G4UniformRand())/cz1) ;
719  }
720  zPathLength = tPathLength*u ;
721  }
722  }
723 
725  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda0 << G4endl;
726  return zPathLength;
727 }
const G4ParticleDefinition * particle
const G4MaterialCutsCouple * couple
G4double currentKinEnergy
G4double UrbanMscModel93::ComputeTheta0 ( G4double  truePathLength,
G4double  KineticEnergy 
)
inline

Definition at line 217 of file UrbanMscModel93.h.

References funct::abs(), charge, coeffth1, coeffth2, currentKinEnergy, currentRadLength, mass, MeV, mathSSE::sqrt(), and y.

Referenced by SampleCosineTheta().

219 {
220  // for all particles take the width of the central part
221  // from a parametrization similar to the Highland formula
222  // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
223  static const G4double c_highland = 13.6*CLHEP::MeV;
224  G4double invbetacp = sqrt((currentKinEnergy+mass)*(KineticEnergy+mass)/
226  KineticEnergy*(KineticEnergy+2.*mass)));
227  y = trueStepLength/currentRadLength;
228  G4double theta0 = c_highland*std::abs(charge)*sqrt(y)*invbetacp;
229  // correction factor from e- scattering data
230  theta0 *= (coeffth1+coeffth2*G4Log(y));
231 
232  return theta0;
233 }
G4double currentRadLength
const double MeV
T sqrt(T t)
Definition: SSEVec.h:48
G4double currentKinEnergy
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
G4double UrbanMscModel93::ComputeTruePathLengthLimit ( const G4Track &  track,
G4double &  currentMinimalStep 
)

Definition at line 428 of file UrbanMscModel93.cc.

References couple, currentKinEnergy, currentMaterialIndex, currentRange, firstStep, fr, geombig, geomlimit, geommin, inside, insideskin, lambda0, lambdalimit, mass, masslimite, MeV, particle, presafety, rangeinit, skindepth, smallstep, stepmin, tgeom, tlimit, tlimitmin, tlimitminfix, and tPathLength.

431 {
432  tPathLength = currentMinimalStep;
433  const G4DynamicParticle* dp = track.GetDynamicParticle();
434  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
435  G4StepStatus stepStatus = sp->GetStepStatus();
436  couple = track.GetMaterialCutsCouple();
437  SetCurrentCouple(couple);
438  currentMaterialIndex = couple->GetIndex();
439  currentKinEnergy = dp->GetKineticEnergy();
441  lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy);
442 
443  // stop here if small range particle
444  if(inside || tPathLength < tlimitminfix) {
445  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
446  }
447 
449 
450  presafety = sp->GetSafety();
451 
452  // G4cout << "Urban2::StepLimit tPathLength= "
453  // <<tPathLength<<" safety= " << presafety
454  // << " range= " <<currentRange<< " lambda= "<<lambda0
455  // << " Alg: " << steppingAlgorithm <<G4endl;
456 
457  // far from geometry boundary
458  if(currentRange < presafety)
459  {
460  inside = true;
461  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
462  }
463 
464  // standard version
465  //
466  if (steppingAlgorithm == fUseDistanceToBoundary)
467  {
468  //compute geomlimit and presafety
469  geomlimit = ComputeGeomLimit(track, presafety, currentRange);
470 
471  // is it far from boundary ?
472  if(currentRange < presafety)
473  {
474  inside = true;
475  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
476  }
477 
478  smallstep += 1.;
479  insideskin = false;
480 
481  if(firstStep || stepStatus == fGeomBoundary)
482  {
484  if(firstStep) smallstep = 1.e10;
485  else smallstep = 1.;
486 
487  //define stepmin here (it depends on lambda!)
488  //rough estimation of lambda_elastic/lambda_transport
489  G4double rat = currentKinEnergy/MeV ;
490  rat = 1.e-3/(rat*(10.+rat)) ;
491  //stepmin ~ lambda_elastic
492  stepmin = rat*lambda0;
493  skindepth = skin*stepmin;
494  //define tlimitmin
495  tlimitmin = 10.*stepmin;
497  //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
498  // << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
499  // constraint from the geometry
500  if((geomlimit < geombig) && (geomlimit > geommin))
501  {
502  // geomlimit is a geometrical step length
503  // transform it to true path length (estimation)
504  if((1.-geomlimit/lambda0) > 0.)
505  geomlimit = -lambda0*G4Log(1.-geomlimit/lambda0)+tlimitmin ;
506 
507  if(stepStatus == fGeomBoundary)
508  tgeom = geomlimit/facgeom;
509  else
510  tgeom = 2.*geomlimit/facgeom;
511  }
512  else
513  tgeom = geombig;
514  }
515 
516 
517  //step limit
518  tlimit = facrange*rangeinit;
519  if(tlimit < facsafety*presafety)
520  tlimit = facsafety*presafety;
521 
522  //lower limit for tlimit
524 
525  if(tlimit > tgeom) tlimit = tgeom;
526 
527  //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
528  // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
529 
530  // shortcut
531  if((tPathLength < tlimit) && (tPathLength < presafety) &&
532  (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
533  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
534 
535  // step reduction near to boundary
536  if(smallstep < skin)
537  {
538  tlimit = stepmin;
539  insideskin = true;
540  }
541  else if(geomlimit < geombig)
542  {
543  if(geomlimit > skindepth)
544  {
545  if(tlimit > geomlimit-0.999*skindepth)
546  tlimit = geomlimit-0.999*skindepth;
547  }
548  else
549  {
550  insideskin = true;
551  if(tlimit > stepmin) tlimit = stepmin;
552  }
553  }
554 
555  if(tlimit < stepmin) tlimit = stepmin;
556 
557  // randomize 1st step or 1st 'normal' step in volume
558  if(firstStep || ((smallstep == skin) && !insideskin))
559  {
560  G4double temptlimit = tlimit;
561  if(temptlimit > tlimitmin)
562  {
563  do {
564  temptlimit = G4RandGauss::shoot(tlimit,0.3*tlimit);
565  } while ((temptlimit < tlimitmin) ||
566  (temptlimit > 2.*tlimit-tlimitmin));
567  }
568  else
569  temptlimit = tlimitmin;
570  if(tPathLength > temptlimit) tPathLength = temptlimit;
571  }
572  else
573  {
575  }
576 
577  }
578  // for 'normal' simulation with or without magnetic field
579  // there no small step/single scattering at boundaries
580  else if(steppingAlgorithm == fUseSafety)
581  {
582  // compute presafety again if presafety <= 0 and no boundary
583  // i.e. when it is needed for optimization purposes
584  if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
585  presafety = ComputeSafety(sp->GetPosition(),tPathLength);
586 
587  // is far from boundary
588  if(currentRange < presafety)
589  {
590  inside = true;
591  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
592  }
593 
594  if(firstStep || stepStatus == fGeomBoundary)
595  {
596  rangeinit = currentRange;
597  fr = facrange;
598  // 9.1 like stepping for e+/e- only (not for muons,hadrons)
599  if(mass < masslimite)
600  {
601  if(lambda0 > currentRange)
602  rangeinit = lambda0;
603  if(lambda0 > lambdalimit)
604  fr *= 0.75+0.25*lambda0/lambdalimit;
605  }
606 
607  //lower limit for tlimit
608  G4double rat = currentKinEnergy/MeV ;
609  rat = 1.e-3/(rat*(10.+rat)) ;
610  tlimitmin = 10.*lambda0*rat;
612  }
613  //step limit
614  tlimit = fr*rangeinit;
615 
616  if(tlimit < facsafety*presafety)
617  tlimit = facsafety*presafety;
618 
619  //lower limit for tlimit
621 
623 
624  }
625 
626  // version similar to 7.1 (needed for some experiments)
627  else
628  {
629  if (stepStatus == fGeomBoundary)
630  {
631  if (currentRange > lambda0) tlimit = facrange*currentRange;
632  else tlimit = facrange*lambda0;
633 
636  }
637  }
638  //G4cout << "tPathLength= " << tPathLength
639  // << " currentMinimalStep= " << currentMinimalStep << G4endl;
640  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
641 }
const G4ParticleDefinition * particle
const G4MaterialCutsCouple * couple
const double MeV
G4double currentKinEnergy
G4double UrbanMscModel93::ComputeTrueStepLength ( G4double  geomStepLength)

Definition at line 731 of file UrbanMscModel93.cc.

References currentRange, insideskin, lambda0, par1, par3, tausmall, tlimitminfix, tPathLength, and zPathLength.

732 {
733  // step defined other than transportation
734  if(geomStepLength == zPathLength && tPathLength <= currentRange)
735  return tPathLength;
736 
737  // t = z for very small step
738  zPathLength = geomStepLength;
739  tPathLength = geomStepLength;
740  if(geomStepLength < tlimitminfix) return tPathLength;
741 
742  // recalculation
743  if((geomStepLength > lambda0*tausmall) && !insideskin)
744  {
745  if(par1 < 0.)
746  tPathLength = -lambda0*G4Log(1.-geomStepLength/lambda0) ;
747  else
748  {
749  if(par1*par3*geomStepLength < 1.)
750  tPathLength = (1.-G4Exp(G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
751  else
753  }
754  }
755  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
756 
757  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
758 
759  return tPathLength;
760 }
void UrbanMscModel93::Initialise ( const G4ParticleDefinition *  p,
const G4DataVector &   
)

Definition at line 150 of file UrbanMscModel93.cc.

References fParticleChange, MeV, SetParticle(), skindepth, and stepmin.

152 {
153  skindepth = skin*stepmin;
154 
155  // set values of some data members
156  SetParticle(p);
157 
158  if(p->GetPDGMass() > MeV) {
159  G4cout << "### WARNING: UrbanMscModel93 model is used for "
160  << p->GetParticleName() << " !!! " << G4endl;
161  G4cout << "### This model should be used only for e+-"
162  << G4endl;
163  }
164 
165  fParticleChange = GetParticleChangeForMSC(p);
166 }
G4ParticleChangeForMSC * fParticleChange
void SetParticle(const G4ParticleDefinition *)
const double MeV
G4double UrbanMscModel93::LatCorrelation ( )
private

Definition at line 1102 of file UrbanMscModel93.cc.

References currentTau, insideskin, kappa, kappami1, kappapl1, lambdaeff, taubig, taulim, tausmall, and third.

Referenced by SampleScattering().

1103 {
1104  G4double latcorr = 0.;
1105  if((currentTau >= tausmall) && !insideskin)
1106  {
1107  if(currentTau < taulim)
1110  else
1111  {
1112  G4double etau = 0.;
1113  if(currentTau < taubig) etau = G4Exp(-currentTau);
1114  latcorr = -kappa*currentTau;
1115  latcorr = G4Exp(latcorr)/kappami1;
1116  latcorr += 1.-kappa*etau/kappami1 ;
1117  latcorr *= 2.*lambdaeff*third ;
1118  }
1119  }
1120 
1121  return latcorr;
1122 }
static const G4double kappami1
static const G4double kappapl1
static const G4double kappa
UrbanMscModel93& UrbanMscModel93::operator= ( const UrbanMscModel93 right)
private
G4double UrbanMscModel93::SampleCosineTheta ( G4double  trueStepLength,
G4double  KineticEnergy 
)
private

Definition at line 856 of file UrbanMscModel93.cc.

References b, EnergyCorrector::c, alignmentValidation::c1, coeffc1, coeffc2, ComputeTheta0(), funct::cos(), couple, currentKinEnergy, currentRadLength, currentRange, currentTau, ztail::d, ea, eaa, GeV, i, insideskin, lambda0, lambdaeff, mass, timingPdfMaker::mean, gen::n, numlim, par1, par2, phi, funct::pow(), rellossmax, scr1, scr2, SimpleScattering(), funct::sin(), mathSSE::sqrt(), stepmin, metsig::tau, taubig, taulim, tausmall, theta0max, third, UpdateCache(), MetTreeProducer::var(), x, xsi, y, Zeff, and Zold.

Referenced by SampleScattering().

858 {
859  G4double cth = 1. ;
860  G4double tau = trueStepLength/lambda0 ;
861 
862  Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
863  couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ;
864 
865  if(Zold != Zeff)
866  UpdateCache();
867 
868  if(insideskin)
869  {
870  //no scattering, single or plural scattering
871  G4double mean = trueStepLength/stepmin ;
872 
873  G4int n = G4Poisson(mean);
874  if(n > 0)
875  {
876  //screening (Moliere-Bethe)
877  G4double mom2 = KineticEnergy*(2.*mass+KineticEnergy);
878  G4double beta2 = mom2/((KineticEnergy+mass)*(KineticEnergy+mass));
879  G4double ascr = scr1/mom2;
880  ascr *= 1.13+scr2/beta2;
881  G4double ascr1 = 1.+2.*ascr;
882  G4double bp1=ascr1+1.;
883  G4double bm1=ascr1-1.;
884 
885  // single scattering from screened Rutherford x-section
886  G4double ct,st,phi;
887  G4double sx=0.,sy=0.,sz=0.;
888  for(G4int i=1; i<=n; i++)
889  {
890  ct = ascr1-bp1*bm1/(2.*G4UniformRand()+bm1);
891  if(ct < -1.) ct = -1.;
892  if(ct > 1.) ct = 1.;
893  st = sqrt(1.-ct*ct);
894  phi = twopi*G4UniformRand();
895  sx += st*cos(phi);
896  sy += st*sin(phi);
897  sz += ct;
898  }
899  cth = sz/sqrt(sx*sx+sy*sy+sz*sz);
900  }
901  }
902  else
903  {
904  if(trueStepLength >= currentRange*dtrl)
905  {
906  if(par1*trueStepLength < 1.)
907  tau = -par2*G4Log(1.-par1*trueStepLength) ;
908  // for the case if ioni/brems are inactivated
909  // see the corresponding condition in ComputeGeomPathLength
910  else if(1.-KineticEnergy/currentKinEnergy > taulim)
911  tau = taubig ;
912  }
913  currentTau = tau ;
914  lambdaeff = trueStepLength/currentTau;
915  currentRadLength = couple->GetMaterial()->GetRadlen();
916 
917  if (tau >= taubig) cth = -1.+2.*G4UniformRand();
918  else if (tau >= tausmall)
919  {
920  G4double xmeanth, x2meanth;
921  if(tau < numlim) {
922  xmeanth = 1.0 - tau*(1.0 - 0.5*tau);
923  x2meanth= 1.0 - tau*(5.0 - 6.25*tau)*third;
924  } else {
925  xmeanth = G4Exp(-tau);
926  x2meanth = (1.+2.*G4Exp(-2.5*tau))*third;
927  }
928  G4double relloss = 1.-KineticEnergy/currentKinEnergy;
929 
930  if(relloss > rellossmax)
931  return SimpleScattering(xmeanth,x2meanth);
932 
933  G4double theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
934 
935  //G4cout << "Theta0= " << theta0 << " theta0max= " << theta0max
936  // << " sqrt(tausmall)= " << sqrt(tausmall) << G4endl;
937 
938  // protection for very small angles
939  G4double theta2 = theta0*theta0;
940 
941  if(theta2 < tausmall) { return cth; }
942 
943  if(theta0 > theta0max) {
944  return SimpleScattering(xmeanth,x2meanth);
945  }
946 
947  G4double x = theta2*(1.0 - theta2/12.);
948  if(theta2 > numlim) {
949  G4double sth = 2.*sin(0.5*theta0);
950  x = sth*sth;
951  }
952 
953  G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/eaa;
954  G4double x0 = 1. - xsi*x;
955 
956  // G4cout << " xmean1= " << xmean1 << " xmeanth= " << xmeanth << G4endl;
957 
958  if(xmean1 <= 0.999*xmeanth) {
959  return SimpleScattering(xmeanth,x2meanth);
960  }
961  // from e- and muon scattering data
962  G4double c = coeffc1+coeffc2*y;
963 
964  // tail should not be too big
965  if(c < 1.9) {
966  /*
967  if(KineticEnergy > 200*MeV && c < 1.6) {
968  G4cout << "UrbanMscModel93::SampleCosineTheta: E(GeV)= "
969  << KineticEnergy/GeV
970  << " !!** c= " << c
971  << " **!! length(mm)= " << trueStepLength << " Zeff= " << Zeff
972  << " " << couple->GetMaterial()->GetName()
973  << " tau= " << tau << G4endl;
974  }
975  */
976  c = 1.9;
977  }
978 
979  if(fabs(c-3.) < 0.001) { c = 3.001; }
980  else if(fabs(c-2.) < 0.001) { c = 2.001; }
981 
982  G4double c1 = c-1.;
983 
984  //from continuity of derivatives
985  G4double b = 1.+(c-xsi)*x;
986 
987  G4double b1 = b+1.;
988  G4double bx = c*x;
989 
990  G4double eb1 = pow(b1,c1);
991  G4double ebx = pow(bx,c1);
992  G4double d = ebx/eb1;
993 
994  // G4double xmean2 = (x0*eb1+ebx-(eb1*bx-b1*ebx)/(c-2.))/(eb1-ebx);
995  G4double xmean2 = (x0 + d - (bx - b1*d)/(c-2.))/(1. - d);
996 
997  G4double f1x0 = ea/eaa;
998  G4double f2x0 = c1/(c*(1. - d));
999  G4double prob = f2x0/(f1x0+f2x0);
1000 
1001  G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
1002 
1003  // sampling of costheta
1004  //G4cout << "c= " << c << " qprob= " << qprob << " eb1= " << eb1
1005  // << " c1= " << c1 << " b1= " << b1 << " bx= " << bx << " eb1= " << eb1
1006  // << G4endl;
1007  if(G4UniformRand() < qprob)
1008  {
1009  G4double var = 0;
1010  if(G4UniformRand() < prob) {
1011  cth = 1.+G4Log(ea+G4UniformRand()*eaa)*x;
1012  } else {
1013  var = (1.0 - d)*G4UniformRand();
1014  if(var < numlim*d) {
1015  var /= (d*c1);
1016  cth = -1.0 + var*(1.0 - 0.5*var*c)*(2. + (c - xsi)*x);
1017  } else {
1018  cth = 1. + x*(c - xsi - c*pow(var + d, -1.0/c1));
1019  //b-b1*bx/G4Exp(log(ebx+(eb1-ebx)*G4UniformRand())/c1) ;
1020  }
1021  }
1022  if(KineticEnergy > 5*GeV && cth < 0.9) {
1023  G4cout << "UrbanMscModel93::SampleCosineTheta: E(GeV)= "
1024  << KineticEnergy/GeV
1025  << " 1-cosT= " << 1 - cth
1026  << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
1027  << " tau= " << tau
1028  << " prob= " << prob << " var= " << var << G4endl;
1029  G4cout << " c= " << c << " qprob= " << qprob << " eb1= " << eb1
1030  << " ebx= " << ebx
1031  << " c1= " << c1 << " b= " << b << " b1= " << b1
1032  << " bx= " << bx << " d= " << d
1033  << " ea= " << ea << " eaa= " << eaa << G4endl;
1034  }
1035  }
1036  else {
1037  cth = -1.+2.*G4UniformRand();
1038  if(KineticEnergy > 5*GeV) {
1039  G4cout << "UrbanMscModel93::SampleCosineTheta: E(GeV)= "
1040  << KineticEnergy/GeV
1041  << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
1042  << " qprob= " << qprob << G4endl;
1043  }
1044  }
1045  }
1046  }
1047  return cth ;
1048 }
const G4MaterialCutsCouple * couple
int i
Definition: DBlmapReader.cc:9
const double GeV
Definition: MathUtil.h:16
G4double ComputeTheta0(G4double truePathLength, G4double KineticEnergy)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
G4double currentRadLength
tuple d
Definition: ztail.py:151
T sqrt(T t)
Definition: SSEVec.h:48
G4double currentKinEnergy
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
G4double SimpleScattering(G4double xmeanth, G4double x2meanth)
double b
Definition: hdecay.h:120
Definition: DDAxes.h:10
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
Definition: DDAxes.h:10
G4double UrbanMscModel93::SampleDisplacement ( )
private

Definition at line 1052 of file UrbanMscModel93.cc.

References currentTau, insideskin, kappa, kappami1, kappapl1, lambdaeff, mathSSE::sqrt(), taubig, taulim, tausmall, third, tPathLength, and zPathLength.

Referenced by SampleScattering().

1053 {
1054  // Compute rmean = sqrt(<r**2>) from theory
1055  G4double rmean = 0.0;
1056  if ((currentTau >= tausmall) && !insideskin) {
1057  if (currentTau < taulim) {
1059  (1.-kappapl1*currentTau*0.25)/6. ;
1060 
1061  } else {
1062  G4double etau = 0.0;
1063  if (currentTau<taubig) etau = G4Exp(-currentTau);
1064  rmean = -kappa*currentTau;
1065  rmean = -G4Exp(rmean)/(kappa*kappami1);
1066  rmean += currentTau-kappapl1/kappa+kappa*etau/kappami1;
1067  }
1068  if (rmean>0.) rmean = 2.*lambdaeff*sqrt(rmean*third);
1069  else rmean = 0.;
1070  }
1071 
1072  if(rmean == 0.) return rmean;
1073 
1074  // protection against z > t ...........................
1075  G4double rmax = (tPathLength-zPathLength)*(tPathLength+zPathLength);
1076  if(rmax <= 0.)
1077  rmax = 0.;
1078  else
1079  rmax = sqrt(rmax);
1080 
1081  if(rmean >= rmax) return rmax;
1082 
1083  return rmean;
1084  // VI comment out for the time being
1085  /*
1086  //sample r (Gaussian distribution with a mean of rmean )
1087  G4double r = 0.;
1088  G4double sigma = min(rmean,rmax-rmean);
1089  sigma /= 3.;
1090  G4double rlow = rmean-3.*sigma;
1091  G4double rhigh = rmean+3.*sigma;
1092  do {
1093  r = G4RandGauss::shoot(rmean,sigma);
1094  } while ((r < rlow) || (r > rhigh));
1095 
1096  return r;
1097  */
1098 }
static const G4double kappami1
T sqrt(T t)
Definition: SSEVec.h:48
static const G4double kappapl1
static const G4double kappa
G4ThreeVector & UrbanMscModel93::SampleScattering ( const G4ThreeVector &  oldDirection,
G4double  safety 
)

Definition at line 765 of file UrbanMscModel93.cc.

References funct::abs(), funct::cos(), couple, currentKinEnergy, currentRange, fParticleChange, lambda0, LatCorrelation(), particle, phi, colinearityKinematic::Phi, alignCSCRings::r, SampleCosineTheta(), SampleDisplacement(), funct::sin(), mathSSE::sqrt(), tausmall, tlimitminfix, and tPathLength.

767 {
768  fDisplacement.set(0.0,0.0,0.0);
769  G4double kineticEnergy = currentKinEnergy;
770  if (tPathLength > currentRange*dtrl) {
771  kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
772  } else {
773  kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple);
774  }
775  if((kineticEnergy <= eV) || (tPathLength <= tlimitminfix) ||
776  (tPathLength/tausmall < lambda0)) { return fDisplacement; }
777 
778  G4double cth = SampleCosineTheta(tPathLength,kineticEnergy);
779 
780  // protection against 'bad' cth values
781  if(std::fabs(cth) > 1.) { return fDisplacement; }
782 
783  // extra protection agaist high energy particles backscattered
784  // if(cth < 1.0 - 1000*tPathLength/lambda0 && kineticEnergy > 20*MeV) {
785  //G4cout << "Warning: large scattering E(MeV)= " << kineticEnergy
786  // << " s(mm)= " << tPathLength/mm
787  // << " 1-cosTheta= " << 1.0 - cth << G4endl;
788  // do Gaussian central scattering
789  // if(kineticEnergy > 0.5*GeV && cth < 0.9) {
790  /*
791  if(cth < 1.0 - 1000*tPathLength/lambda0 &&
792  cth < 0.9 && kineticEnergy > 500*MeV) {
793  G4ExceptionDescription ed;
794  ed << particle->GetParticleName()
795  << " E(MeV)= " << kineticEnergy/MeV
796  << " Step(mm)= " << tPathLength/mm
797  << " tau= " << tPathLength/lambda0
798  << " in " << CurrentCouple()->GetMaterial()->GetName()
799  << " CosTheta= " << cth
800  << " is too big";
801  G4Exception("UrbanMscModel93::SampleScattering","em0004",
802  JustWarning, ed,"");
803  }
804  */
805 
806  G4double sth = sqrt((1.0 - cth)*(1.0 + cth));
807  G4double phi = twopi*G4UniformRand();
808  G4double dirx = sth*cos(phi);
809  G4double diry = sth*sin(phi);
810 
811  G4ThreeVector newDirection(dirx,diry,cth);
812  newDirection.rotateUz(oldDirection);
813  fParticleChange->ProposeMomentumDirection(newDirection);
814 
815  if (latDisplasment && safety > tlimitminfix) {
816 
817  G4double r = SampleDisplacement();
818  /*
819  G4cout << "UrbanMscModel93::SampleSecondaries: e(MeV)= " << kineticEnergy
820  << " sinTheta= " << sth << " r(mm)= " << r
821  << " trueStep(mm)= " << tPathLength
822  << " geomStep(mm)= " << zPathLength
823  << G4endl;
824  */
825  if(r > 0.)
826  {
827  G4double latcorr = LatCorrelation();
828  if(latcorr > r) latcorr = r;
829 
830  // sample direction of lateral displacement
831  // compute it from the lateral correlation
832  G4double Phi = 0.;
833  if(std::abs(r*sth) < latcorr)
834  Phi = twopi*G4UniformRand();
835  else
836  {
837  G4double psi = std::acos(latcorr/(r*sth));
838  if(G4UniformRand() < 0.5)
839  Phi = phi+psi;
840  else
841  Phi = phi-psi;
842  }
843 
844  dirx = r*std::cos(Phi);
845  diry = r*std::sin(Phi);
846 
847  fDisplacement.set(dirx,diry,0.0);
848  fDisplacement.rotateUz(oldDirection);
849  }
850  }
851  return fDisplacement;
852 }
const G4ParticleDefinition * particle
const G4MaterialCutsCouple * couple
G4ParticleChangeForMSC * fParticleChange
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::map< std::string, int, std::less< std::string > > psi
G4double LatCorrelation()
T sqrt(T t)
Definition: SSEVec.h:48
G4double currentKinEnergy
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
G4double SampleDisplacement()
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy)
Definition: DDAxes.h:10
void UrbanMscModel93::SetParticle ( const G4ParticleDefinition *  p)
inlineprivate

Definition at line 182 of file UrbanMscModel93.h.

References charge, ChargeSquare, mass, AlCaHLTBitMon_ParallelJobs::p, and particle.

Referenced by ComputeCrossSectionPerAtom(), Initialise(), and StartTracking().

183 {
184  if (p != particle) {
185  particle = p;
186  mass = p->GetPDGMass();
187  charge = p->GetPDGCharge()/CLHEP::eplus;
189  }
190 }
const G4ParticleDefinition * particle
G4double UrbanMscModel93::SimpleScattering ( G4double  xmeanth,
G4double  x2meanth 
)
inlineprivate

Definition at line 238 of file UrbanMscModel93.h.

References a.

Referenced by SampleCosineTheta().

239 {
240  // 'large angle scattering'
241  // 2 model functions with correct xmean and x2mean
242  G4double a = (2.*xmeanth+9.*x2meanth-3.)/(2.*xmeanth-3.*x2meanth+1.);
243  G4double prob = (a+2.)*xmeanth/a;
244 
245  // sampling
246  G4double cth = 1.;
247  if(G4UniformRand() < prob) {
248  cth = -1.+2.*G4Exp(G4Log(G4UniformRand())/(a+1.));
249  } else {
250  cth = -1.+2.*G4UniformRand();
251  }
252  return cth;
253 }
double a
Definition: hdecay.h:121
void UrbanMscModel93::StartTracking ( G4Track *  track)

Definition at line 415 of file UrbanMscModel93.cc.

References firstStep, geombig, inside, insideskin, SetParticle(), stepmin, tlimit, tlimitmin, and tlimitminfix.

416 {
417  SetParticle(track->GetDynamicParticle()->GetDefinition());
418  firstStep = true;
419  inside = false;
420  insideskin = false;
421  tlimit = geombig;
423  tlimitmin = 10.*stepmin ;
424 }
void SetParticle(const G4ParticleDefinition *)
void UrbanMscModel93::UpdateCache ( )
inlineprivate

Definition at line 195 of file UrbanMscModel93.h.

References ChargeSquare, coeffc1, coeffc2, coeffth1, coeffth2, lnZ, scr1, scr1ini, scr2, scr2ini, Z2, Z23, Zeff, and Zold.

Referenced by SampleCosineTheta().

196 {
197  lnZ = G4Log(Zeff);
198  //new correction in theta0 formula
199  coeffth1 = (1.-8.7780e-2/Zeff)*(0.87+0.03*lnZ);
200  coeffth2 = (4.0780e-2+1.7315e-4*Zeff)*(0.87+0.03*lnZ);
201  // tail parameters
202  G4double lnZ1 = G4Log(Zeff+1.);
203  coeffc1 = 2.943-0.197*lnZ1;
204  coeffc2 = 0.0987-0.0143*lnZ1;
205  // for single scattering
206  Z2 = Zeff*Zeff;
207  Z23 = G4Exp(2.*lnZ/3.);
208  scr1 = scr1ini*Z23;
210 
211  Zold = Zeff;
212 }

Member Data Documentation

G4double UrbanMscModel93::charge
private
G4double UrbanMscModel93::ChargeSquare
private
G4double UrbanMscModel93::coeffc1
private

Definition at line 170 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::coeffc2
private

Definition at line 170 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::coeffth1
private

Definition at line 169 of file UrbanMscModel93.h.

Referenced by ComputeTheta0(), UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::coeffth2
private

Definition at line 169 of file UrbanMscModel93.h.

Referenced by ComputeTheta0(), UpdateCache(), and UrbanMscModel93().

const G4MaterialCutsCouple* UrbanMscModel93::couple
private
G4double UrbanMscModel93::currentKinEnergy
private
G4int UrbanMscModel93::currentMaterialIndex
private

Definition at line 164 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::currentRadLength
private

Definition at line 157 of file UrbanMscModel93.h.

Referenced by ComputeTheta0(), SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::currentRange
private
G4double UrbanMscModel93::currentTau
private
G4double UrbanMscModel93::ea
private

Definition at line 159 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::eaa
private

Definition at line 159 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4bool UrbanMscModel93::firstStep
private
G4ParticleChangeForMSC* UrbanMscModel93::fParticleChange
private

Definition at line 120 of file UrbanMscModel93.h.

Referenced by Initialise(), SampleScattering(), and UrbanMscModel93().

G4double UrbanMscModel93::fr
private

Definition at line 127 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::geombig
private

Definition at line 138 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), StartTracking(), and UrbanMscModel93().

G4double UrbanMscModel93::geomlimit
private

Definition at line 140 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::geommin
private

Definition at line 139 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4bool UrbanMscModel93::inside
private

Definition at line 174 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), StartTracking(), and UrbanMscModel93().

G4bool UrbanMscModel93::insideskin
private
G4double UrbanMscModel93::lambda0
private
G4double UrbanMscModel93::lambdaeff
private
G4double UrbanMscModel93::lambdalimit
private

Definition at line 127 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::lnZ
private

Definition at line 168 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::mass
private
G4double UrbanMscModel93::masslimite
private

Definition at line 127 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::numlim
private

Definition at line 159 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::par1
private
G4double UrbanMscModel93::par2
private

Definition at line 150 of file UrbanMscModel93.h.

Referenced by ComputeGeomPathLength(), SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::par3
private
const G4ParticleDefinition* UrbanMscModel93::particle
private
G4double UrbanMscModel93::presafety
private

Definition at line 144 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::rangeinit
private

Definition at line 156 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::rellossmax
private

Definition at line 161 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::scr1
private

Definition at line 171 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::scr1ini
private

Definition at line 171 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::scr2
private

Definition at line 171 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::scr2ini
private

Definition at line 171 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::skindepth
private

Definition at line 141 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), Initialise(), and UrbanMscModel93().

G4double UrbanMscModel93::smallstep
private

Definition at line 142 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::stepmin
private
G4double UrbanMscModel93::taubig
private
G4double UrbanMscModel93::taulim
private
G4double UrbanMscModel93::tausmall
private
G4double UrbanMscModel93::tgeom
private

Definition at line 136 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4LossTableManager* UrbanMscModel93::theManager
private

Definition at line 123 of file UrbanMscModel93.h.

Referenced by UrbanMscModel93().

G4double UrbanMscModel93::theta0max
private

Definition at line 161 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::third
private
G4double UrbanMscModel93::tlimit
private

Definition at line 133 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), StartTracking(), and UrbanMscModel93().

G4double UrbanMscModel93::tlimitmin
private

Definition at line 134 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), StartTracking(), and UrbanMscModel93().

G4double UrbanMscModel93::tlimitminfix
private
G4double UrbanMscModel93::tPathLength
private
G4double UrbanMscModel93::xsi
private

Definition at line 159 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::y
private
G4double UrbanMscModel93::Z2
private

Definition at line 168 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::Z23
private

Definition at line 168 of file UrbanMscModel93.h.

Referenced by ComputeCrossSectionPerAtom(), UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::Zeff
private

Definition at line 168 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::Zold
private

Definition at line 167 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::zPathLength
private