CMS 3D CMS Logo

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 59 of file UrbanMscModel93.h.

Constructor & Destructor Documentation

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

Definition at line 38 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.

39  : G4VMscModel(nam)
40 {
41  masslimite = 0.6*MeV;
42  lambdalimit = 1.*mm;
43  fr = 0.02;
44  taubig = 8.0;
45  tausmall = 1.e-16;
46  taulim = 1.e-6;
48  tlimitminfix = 1.e-6*mm;
50  smallstep = 1.e10;
51  currentRange = 0. ;
52  rangeinit = 0.;
53  tlimit = 1.e10*mm;
54  tlimitmin = 10.*tlimitminfix;
55  tgeom = 1.e50*mm;
56  geombig = 1.e50*mm;
57  geommin = 1.e-3*mm;
59  presafety = 0.*mm;
60 
61  y = 0.;
62 
63  Zold = 0.;
64  Zeff = 1.;
65  Z2 = 1.;
66  Z23 = 1.;
67  lnZ = 0.;
68  coeffth1 = 0.;
69  coeffth2 = 0.;
70  coeffc1 = 0.;
71  coeffc2 = 0.;
72  scr1ini = fine_structure_const*fine_structure_const*
73  electron_mass_c2*electron_mass_c2/(0.885*0.885*4.*pi);
74  scr2ini = 3.76*fine_structure_const*fine_structure_const;
75  scr1 = 0.;
76  scr2 = 0.;
77 
78  theta0max = pi/6.;
79  rellossmax = 0.50;
80  third = 1./3.;
81  particle = 0;
82  theManager = G4LossTableManager::Instance();
83  firstStep = true;
84  inside = false;
85  insideskin = false;
86 
87  numlim = 0.01;
88  xsi = 3.;
89  ea = G4Exp(-xsi);
90  eaa = 1.-ea ;
91 
92  skindepth = skin*stepmin;
93 
94  mass = proton_mass_c2;
95  charge = ChargeSquare = 1.0;
97  = zPathLength = par1 = par2 = par3 = 0;
98 
100  fParticleChange = 0;
101  couple = 0;
102  SetSampleZ(false);
103 }
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 107 of file UrbanMscModel93.cc.

108 {}
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 132 of file UrbanMscModel93.cc.

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

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

608 {
609  firstStep = false;
610  lambdaeff = lambda0;
611  par1 = -1. ;
612  par2 = par3 = 0. ;
613 
614  // do the true -> geom transformation
616 
617  // z = t for very small tPathLength
618  if(tPathLength < tlimitminfix) return zPathLength;
619 
620  // this correction needed to run MSC with eIoni and eBrem inactivated
621  // and makes no harm for a normal run
624 
625  G4double tau = tPathLength/lambda0 ;
626 
627  if ((tau <= tausmall) || insideskin) {
630  return zPathLength;
631  }
632 
633  G4double zmean;
634  if (tPathLength < currentRange*dtrl) {
635  if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
636  else zmean = lambda0*(1.-G4Exp(-tau));
637  } else if(currentKinEnergy < mass || tPathLength == currentRange) {
638  par1 = 1./currentRange ;
639  par2 = 1./(par1*lambda0) ;
640  par3 = 1.+par2 ;
642  zmean = (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3) ;
643  else
644  zmean = 1./(par1*par3) ;
645  } else {
646  G4double T1 = GetEnergy(particle,currentRange-tPathLength,couple);
647  G4double lambda1 = GetTransportMeanFreePath(particle,T1);
648 
649  par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
650  par2 = 1./(par1*lambda0) ;
651  par3 = 1.+par2 ;
652  zmean = (1.-G4Exp(par3*G4Log(lambda1/lambda0)))/(par1*par3) ;
653  }
654 
655  zPathLength = zmean ;
656 
657  // sample z
658  if(samplez)
659  {
660  const G4double ztmax = 0.99 ;
661  G4double zt = zmean/tPathLength ;
662 
663  if (tPathLength > stepmin && zt < ztmax)
664  {
665  G4double u,cz1;
666  if(zt >= third)
667  {
668  G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
669  cz1 = 1.+cz ;
670  G4double u0 = cz/cz1 ;
671  G4double grej ;
672  do {
673  u = G4Exp(G4Log(G4UniformRand())/cz1) ;
674  grej = G4Exp(cz*G4Log(u/u0))*(1.-u)/(1.-u0) ;
675  } while (grej < G4UniformRand()) ;
676  }
677  else
678  {
679  cz1 = 1./zt-1.;
680  u = 1.-G4Exp(G4Log(G4UniformRand())/cz1) ;
681  }
682  zPathLength = tPathLength*u ;
683  }
684  }
685 
687  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda0 << G4endl;
688  return zPathLength;
689 }
const G4ParticleDefinition * particle
const G4MaterialCutsCouple * couple
G4double currentKinEnergy
G4double UrbanMscModel93::ComputeTheta0 ( G4double  truePathLength,
G4double  KineticEnergy 
)
inline

Definition at line 207 of file UrbanMscModel93.h.

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

Referenced by SampleCosineTheta().

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

Definition at line 390 of file UrbanMscModel93.cc.

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

393 {
394  tPathLength = currentMinimalStep;
395  const G4DynamicParticle* dp = track.GetDynamicParticle();
396  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
397  G4StepStatus stepStatus = sp->GetStepStatus();
398  couple = track.GetMaterialCutsCouple();
399  SetCurrentCouple(couple);
400  currentMaterialIndex = couple->GetIndex();
401  currentKinEnergy = dp->GetKineticEnergy();
403  lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy);
404 
405  // stop here if small range particle
406  if(inside || tPathLength < tlimitminfix) {
407  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
408  }
409 
411 
412  presafety = sp->GetSafety();
413 
414  // G4cout << "Urban2::StepLimit tPathLength= "
415  // <<tPathLength<<" safety= " << presafety
416  // << " range= " <<currentRange<< " lambda= "<<lambda0
417  // << " Alg: " << steppingAlgorithm <<G4endl;
418 
419  // far from geometry boundary
420  if(currentRange < presafety)
421  {
422  inside = true;
423  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
424  }
425 
426  // standard version
427  //
428  if (steppingAlgorithm == fUseDistanceToBoundary)
429  {
430  //compute geomlimit and presafety
431  geomlimit = ComputeGeomLimit(track, presafety, currentRange);
432 
433  // is it far from boundary ?
434  if(currentRange < presafety)
435  {
436  inside = true;
437  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
438  }
439 
440  smallstep += 1.;
441  insideskin = false;
442 
443  if(firstStep || stepStatus == fGeomBoundary)
444  {
446  if(firstStep) smallstep = 1.e10;
447  else smallstep = 1.;
448 
449  //define stepmin here (it depends on lambda!)
450  //rough estimation of lambda_elastic/lambda_transport
451  G4double rat = currentKinEnergy/MeV ;
452  rat = 1.e-3/(rat*(10.+rat)) ;
453  //stepmin ~ lambda_elastic
454  stepmin = rat*lambda0;
455  skindepth = skin*stepmin;
456  //define tlimitmin
457  tlimitmin = 10.*stepmin;
459  //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
460  // << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
461  // constraint from the geometry
462  if((geomlimit < geombig) && (geomlimit > geommin))
463  {
464  // geomlimit is a geometrical step length
465  // transform it to true path length (estimation)
466  if((1.-geomlimit/lambda0) > 0.)
467  geomlimit = -lambda0*G4Log(1.-geomlimit/lambda0)+tlimitmin ;
468 
469  if(stepStatus == fGeomBoundary)
470  tgeom = geomlimit/facgeom;
471  else
472  tgeom = 2.*geomlimit/facgeom;
473  }
474  else
475  tgeom = geombig;
476  }
477 
478 
479  //step limit
480  tlimit = facrange*rangeinit;
481  if(tlimit < facsafety*presafety)
482  tlimit = facsafety*presafety;
483 
484  //lower limit for tlimit
486 
487  if(tlimit > tgeom) tlimit = tgeom;
488 
489  //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
490  // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
491 
492  // shortcut
493  if((tPathLength < tlimit) && (tPathLength < presafety) &&
494  (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
495  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
496 
497  // step reduction near to boundary
498  if(smallstep < skin)
499  {
500  tlimit = stepmin;
501  insideskin = true;
502  }
503  else if(geomlimit < geombig)
504  {
505  if(geomlimit > skindepth)
506  {
507  if(tlimit > geomlimit-0.999*skindepth)
508  tlimit = geomlimit-0.999*skindepth;
509  }
510  else
511  {
512  insideskin = true;
513  if(tlimit > stepmin) tlimit = stepmin;
514  }
515  }
516 
517  if(tlimit < stepmin) tlimit = stepmin;
518 
519  // randomize 1st step or 1st 'normal' step in volume
520  if(firstStep || ((smallstep == skin) && !insideskin))
521  {
522  G4double temptlimit = tlimit;
523  if(temptlimit > tlimitmin)
524  {
525  do {
526  temptlimit = G4RandGauss::shoot(tlimit,0.3*tlimit);
527  } while ((temptlimit < tlimitmin) ||
528  (temptlimit > 2.*tlimit-tlimitmin));
529  }
530  else
531  temptlimit = tlimitmin;
532  if(tPathLength > temptlimit) tPathLength = temptlimit;
533  }
534  else
535  {
537  }
538 
539  }
540  // for 'normal' simulation with or without magnetic field
541  // there no small step/single scattering at boundaries
542  else if(steppingAlgorithm == fUseSafety)
543  {
544  // compute presafety again if presafety <= 0 and no boundary
545  // i.e. when it is needed for optimization purposes
546  if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
547  presafety = ComputeSafety(sp->GetPosition(),tPathLength);
548 
549  // is far from boundary
550  if(currentRange < presafety)
551  {
552  inside = true;
553  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
554  }
555 
556  if(firstStep || stepStatus == fGeomBoundary)
557  {
558  rangeinit = currentRange;
559  fr = facrange;
560  // 9.1 like stepping for e+/e- only (not for muons,hadrons)
561  if(mass < masslimite)
562  {
563  if(lambda0 > currentRange)
564  rangeinit = lambda0;
565  if(lambda0 > lambdalimit)
566  fr *= 0.75+0.25*lambda0/lambdalimit;
567  }
568 
569  //lower limit for tlimit
570  G4double rat = currentKinEnergy/MeV ;
571  rat = 1.e-3/(rat*(10.+rat)) ;
572  tlimitmin = 10.*lambda0*rat;
574  }
575  //step limit
576  tlimit = fr*rangeinit;
577 
578  if(tlimit < facsafety*presafety)
579  tlimit = facsafety*presafety;
580 
581  //lower limit for tlimit
583 
585 
586  }
587 
588  // version similar to 7.1 (needed for some experiments)
589  else
590  {
591  if (stepStatus == fGeomBoundary)
592  {
593  if (currentRange > lambda0) tlimit = facrange*currentRange;
594  else tlimit = facrange*lambda0;
595 
598  }
599  }
600  //G4cout << "tPathLength= " << tPathLength
601  // << " currentMinimalStep= " << currentMinimalStep << G4endl;
602  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
603 }
const G4ParticleDefinition * particle
const G4MaterialCutsCouple * couple
const double MeV
G4double currentKinEnergy
auto dp
Definition: deltaR.h:22
Definition: sp.h:21
G4double UrbanMscModel93::ComputeTrueStepLength ( G4double  geomStepLength)

Definition at line 693 of file UrbanMscModel93.cc.

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

694 {
695  // step defined other than transportation
696  if(geomStepLength == zPathLength && tPathLength <= currentRange)
697  return tPathLength;
698 
699  // t = z for very small step
700  zPathLength = geomStepLength;
701  tPathLength = geomStepLength;
702  if(geomStepLength < tlimitminfix) return tPathLength;
703 
704  // recalculation
705  if((geomStepLength > lambda0*tausmall) && !insideskin)
706  {
707  if(par1 < 0.)
708  tPathLength = -lambda0*G4Log(1.-geomStepLength/lambda0) ;
709  else
710  {
711  if(par1*par3*geomStepLength < 1.)
712  tPathLength = (1.-G4Exp(G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
713  else
715  }
716  }
717  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
718 
719  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
720 
721  return tPathLength;
722 }
void UrbanMscModel93::Initialise ( const G4ParticleDefinition *  p,
const G4DataVector &   
)

Definition at line 112 of file UrbanMscModel93.cc.

References fParticleChange, ecalTB2006H4_GenSimDigiReco_cfg::G4cout, MeV, SetParticle(), skindepth, and stepmin.

114 {
115  skindepth = skin*stepmin;
116 
117  // set values of some data members
118  SetParticle(p);
119 
120  if(p->GetPDGMass() > MeV) {
121  G4cout << "### WARNING: UrbanMscModel93 model is used for "
122  << p->GetParticleName() << " !!! " << G4endl;
123  G4cout << "### This model should be used only for e+-"
124  << G4endl;
125  }
126 
127  fParticleChange = GetParticleChangeForMSC(p);
128 }
G4ParticleChangeForMSC * fParticleChange
void SetParticle(const G4ParticleDefinition *)
const double MeV
G4double UrbanMscModel93::LatCorrelation ( )
private

Definition at line 1064 of file UrbanMscModel93.cc.

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

Referenced by SampleScattering().

1065 {
1066  G4double latcorr = 0.;
1067  if((currentTau >= tausmall) && !insideskin)
1068  {
1069  if(currentTau < taulim)
1072  else
1073  {
1074  G4double etau = 0.;
1075  if(currentTau < taubig) etau = G4Exp(-currentTau);
1076  latcorr = -kappa*currentTau;
1077  latcorr = G4Exp(latcorr)/kappami1;
1078  latcorr += 1.-kappa*etau/kappami1 ;
1079  latcorr *= 2.*lambdaeff*third ;
1080  }
1081  }
1082 
1083  return latcorr;
1084 }
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 818 of file UrbanMscModel93.cc.

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

Referenced by SampleScattering().

820 {
821  G4double cth = 1. ;
822  G4double tau = trueStepLength/lambda0 ;
823 
824  Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
825  couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ;
826 
827  if(Zold != Zeff)
828  UpdateCache();
829 
830  if(insideskin)
831  {
832  //no scattering, single or plural scattering
833  G4double mean = trueStepLength/stepmin ;
834 
835  G4int n = G4Poisson(mean);
836  if(n > 0)
837  {
838  //screening (Moliere-Bethe)
839  G4double mom2 = KineticEnergy*(2.*mass+KineticEnergy);
840  G4double beta2 = mom2/((KineticEnergy+mass)*(KineticEnergy+mass));
841  G4double ascr = scr1/mom2;
842  ascr *= 1.13+scr2/beta2;
843  G4double ascr1 = 1.+2.*ascr;
844  G4double bp1=ascr1+1.;
845  G4double bm1=ascr1-1.;
846 
847  // single scattering from screened Rutherford x-section
848  G4double ct,st,phi;
849  G4double sx=0.,sy=0.,sz=0.;
850  for(G4int i=1; i<=n; i++)
851  {
852  ct = ascr1-bp1*bm1/(2.*G4UniformRand()+bm1);
853  if(ct < -1.) ct = -1.;
854  if(ct > 1.) ct = 1.;
855  st = sqrt(1.-ct*ct);
856  phi = twopi*G4UniformRand();
857  sx += st*cos(phi);
858  sy += st*sin(phi);
859  sz += ct;
860  }
861  cth = sz/sqrt(sx*sx+sy*sy+sz*sz);
862  }
863  }
864  else
865  {
866  if(trueStepLength >= currentRange*dtrl)
867  {
868  if(par1*trueStepLength < 1.)
869  tau = -par2*G4Log(1.-par1*trueStepLength) ;
870  // for the case if ioni/brems are inactivated
871  // see the corresponding condition in ComputeGeomPathLength
872  else if(1.-KineticEnergy/currentKinEnergy > taulim)
873  tau = taubig ;
874  }
875  currentTau = tau ;
876  lambdaeff = trueStepLength/currentTau;
877  currentRadLength = couple->GetMaterial()->GetRadlen();
878 
879  if (tau >= taubig) cth = -1.+2.*G4UniformRand();
880  else if (tau >= tausmall)
881  {
882  G4double xmeanth, x2meanth;
883  if(tau < numlim) {
884  xmeanth = 1.0 - tau*(1.0 - 0.5*tau);
885  x2meanth= 1.0 - tau*(5.0 - 6.25*tau)*third;
886  } else {
887  xmeanth = G4Exp(-tau);
888  x2meanth = (1.+2.*G4Exp(-2.5*tau))*third;
889  }
890  G4double relloss = 1.-KineticEnergy/currentKinEnergy;
891 
892  if(relloss > rellossmax)
893  return SimpleScattering(xmeanth,x2meanth);
894 
895  G4double theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
896 
897  //G4cout << "Theta0= " << theta0 << " theta0max= " << theta0max
898  // << " sqrt(tausmall)= " << sqrt(tausmall) << G4endl;
899 
900  // protection for very small angles
901  G4double theta2 = theta0*theta0;
902 
903  if(theta2 < tausmall) { return cth; }
904 
905  if(theta0 > theta0max) {
906  return SimpleScattering(xmeanth,x2meanth);
907  }
908 
909  G4double x = theta2*(1.0 - theta2/12.);
910  if(theta2 > numlim) {
911  G4double sth = 2.*sin(0.5*theta0);
912  x = sth*sth;
913  }
914 
915  G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/eaa;
916  G4double x0 = 1. - xsi*x;
917 
918  // G4cout << " xmean1= " << xmean1 << " xmeanth= " << xmeanth << G4endl;
919 
920  if(xmean1 <= 0.999*xmeanth) {
921  return SimpleScattering(xmeanth,x2meanth);
922  }
923  // from e- and muon scattering data
924  G4double c = coeffc1+coeffc2*y;
925 
926  // tail should not be too big
927  if(c < 1.9) {
928  /*
929  if(KineticEnergy > 200*MeV && c < 1.6) {
930  G4cout << "UrbanMscModel93::SampleCosineTheta: E(GeV)= "
931  << KineticEnergy/GeV
932  << " !!** c= " << c
933  << " **!! length(mm)= " << trueStepLength << " Zeff= " << Zeff
934  << " " << couple->GetMaterial()->GetName()
935  << " tau= " << tau << G4endl;
936  }
937  */
938  c = 1.9;
939  }
940 
941  if(fabs(c-3.) < 0.001) { c = 3.001; }
942  else if(fabs(c-2.) < 0.001) { c = 2.001; }
943 
944  G4double c1 = c-1.;
945 
946  //from continuity of derivatives
947  G4double b = 1.+(c-xsi)*x;
948 
949  G4double b1 = b+1.;
950  G4double bx = c*x;
951 
952  G4double eb1 = pow(b1,c1);
953  G4double ebx = pow(bx,c1);
954  G4double d = ebx/eb1;
955 
956  // G4double xmean2 = (x0*eb1+ebx-(eb1*bx-b1*ebx)/(c-2.))/(eb1-ebx);
957  G4double xmean2 = (x0 + d - (bx - b1*d)/(c-2.))/(1. - d);
958 
959  G4double f1x0 = ea/eaa;
960  G4double f2x0 = c1/(c*(1. - d));
961  G4double prob = f2x0/(f1x0+f2x0);
962 
963  G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
964 
965  // sampling of costheta
966  //G4cout << "c= " << c << " qprob= " << qprob << " eb1= " << eb1
967  // << " c1= " << c1 << " b1= " << b1 << " bx= " << bx << " eb1= " << eb1
968  // << G4endl;
969  if(G4UniformRand() < qprob)
970  {
971  G4double var = 0;
972  if(G4UniformRand() < prob) {
973  cth = 1.+G4Log(ea+G4UniformRand()*eaa)*x;
974  } else {
975  var = (1.0 - d)*G4UniformRand();
976  if(var < numlim*d) {
977  var /= (d*c1);
978  cth = -1.0 + var*(1.0 - 0.5*var*c)*(2. + (c - xsi)*x);
979  } else {
980  cth = 1. + x*(c - xsi - c*pow(var + d, -1.0/c1));
981  //b-b1*bx/G4Exp(log(ebx+(eb1-ebx)*G4UniformRand())/c1) ;
982  }
983  }
984  if(KineticEnergy > 5*GeV && cth < 0.9) {
985  G4cout << "UrbanMscModel93::SampleCosineTheta: E(GeV)= "
986  << KineticEnergy/GeV
987  << " 1-cosT= " << 1 - cth
988  << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
989  << " tau= " << tau
990  << " prob= " << prob << " var= " << var << G4endl;
991  G4cout << " c= " << c << " qprob= " << qprob << " eb1= " << eb1
992  << " ebx= " << ebx
993  << " c1= " << c1 << " b= " << b << " b1= " << b1
994  << " bx= " << bx << " d= " << d
995  << " ea= " << ea << " eaa= " << eaa << G4endl;
996  }
997  }
998  else {
999  cth = -1.+2.*G4UniformRand();
1000  if(KineticEnergy > 5*GeV) {
1001  G4cout << "UrbanMscModel93::SampleCosineTheta: E(GeV)= "
1002  << KineticEnergy/GeV
1003  << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
1004  << " qprob= " << qprob << G4endl;
1005  }
1006  }
1007  }
1008  }
1009  return cth ;
1010 }
const G4MaterialCutsCouple * couple
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
T sqrt(T t)
Definition: SSEVec.h:18
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
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
G4double UrbanMscModel93::SampleDisplacement ( )
private

Definition at line 1014 of file UrbanMscModel93.cc.

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

Referenced by SampleScattering().

1015 {
1016  // Compute rmean = sqrt(<r**2>) from theory
1017  G4double rmean = 0.0;
1018  if ((currentTau >= tausmall) && !insideskin) {
1019  if (currentTau < taulim) {
1021  (1.-kappapl1*currentTau*0.25)/6. ;
1022 
1023  } else {
1024  G4double etau = 0.0;
1025  if (currentTau<taubig) etau = G4Exp(-currentTau);
1026  rmean = -kappa*currentTau;
1027  rmean = -G4Exp(rmean)/(kappa*kappami1);
1028  rmean += currentTau-kappapl1/kappa+kappa*etau/kappami1;
1029  }
1030  if (rmean>0.) rmean = 2.*lambdaeff*sqrt(rmean*third);
1031  else rmean = 0.;
1032  }
1033 
1034  if(rmean == 0.) return rmean;
1035 
1036  // protection against z > t ...........................
1037  G4double rmax = (tPathLength-zPathLength)*(tPathLength+zPathLength);
1038  if(rmax <= 0.)
1039  rmax = 0.;
1040  else
1041  rmax = sqrt(rmax);
1042 
1043  if(rmean >= rmax) return rmax;
1044 
1045  return rmean;
1046  // VI comment out for the time being
1047  /*
1048  //sample r (Gaussian distribution with a mean of rmean )
1049  G4double r = 0.;
1050  G4double sigma = min(rmean,rmax-rmean);
1051  sigma /= 3.;
1052  G4double rlow = rmean-3.*sigma;
1053  G4double rhigh = rmean+3.*sigma;
1054  do {
1055  r = G4RandGauss::shoot(rmean,sigma);
1056  } while ((r < rlow) || (r > rhigh));
1057 
1058  return r;
1059  */
1060 }
static const G4double kappami1
T sqrt(T t)
Definition: SSEVec.h:18
static const G4double kappapl1
static const G4double kappa
G4ThreeVector & UrbanMscModel93::SampleScattering ( const G4ThreeVector &  oldDirection,
G4double  safety 
)

Definition at line 727 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.

729 {
730  fDisplacement.set(0.0,0.0,0.0);
731  G4double kineticEnergy = currentKinEnergy;
732  if (tPathLength > currentRange*dtrl) {
733  kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
734  } else {
735  kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple);
736  }
737  if((kineticEnergy <= eV) || (tPathLength <= tlimitminfix) ||
738  (tPathLength/tausmall < lambda0)) { return fDisplacement; }
739 
740  G4double cth = SampleCosineTheta(tPathLength,kineticEnergy);
741 
742  // protection against 'bad' cth values
743  if(std::fabs(cth) > 1.) { return fDisplacement; }
744 
745  // extra protection agaist high energy particles backscattered
746  // if(cth < 1.0 - 1000*tPathLength/lambda0 && kineticEnergy > 20*MeV) {
747  //G4cout << "Warning: large scattering E(MeV)= " << kineticEnergy
748  // << " s(mm)= " << tPathLength/mm
749  // << " 1-cosTheta= " << 1.0 - cth << G4endl;
750  // do Gaussian central scattering
751  // if(kineticEnergy > 0.5*GeV && cth < 0.9) {
752  /*
753  if(cth < 1.0 - 1000*tPathLength/lambda0 &&
754  cth < 0.9 && kineticEnergy > 500*MeV) {
755  G4ExceptionDescription ed;
756  ed << particle->GetParticleName()
757  << " E(MeV)= " << kineticEnergy/MeV
758  << " Step(mm)= " << tPathLength/mm
759  << " tau= " << tPathLength/lambda0
760  << " in " << CurrentCouple()->GetMaterial()->GetName()
761  << " CosTheta= " << cth
762  << " is too big";
763  G4Exception("UrbanMscModel93::SampleScattering","em0004",
764  JustWarning, ed,"");
765  }
766  */
767 
768  G4double sth = sqrt((1.0 - cth)*(1.0 + cth));
769  G4double phi = twopi*G4UniformRand();
770  G4double dirx = sth*cos(phi);
771  G4double diry = sth*sin(phi);
772 
773  G4ThreeVector newDirection(dirx,diry,cth);
774  newDirection.rotateUz(oldDirection);
775  fParticleChange->ProposeMomentumDirection(newDirection);
776 
777  if (latDisplasment && safety > tlimitminfix) {
778 
779  G4double r = SampleDisplacement();
780  /*
781  G4cout << "UrbanMscModel93::SampleSecondaries: e(MeV)= " << kineticEnergy
782  << " sinTheta= " << sth << " r(mm)= " << r
783  << " trueStep(mm)= " << tPathLength
784  << " geomStep(mm)= " << zPathLength
785  << G4endl;
786  */
787  if(r > 0.)
788  {
789  G4double latcorr = LatCorrelation();
790  if(latcorr > r) latcorr = r;
791 
792  // sample direction of lateral displacement
793  // compute it from the lateral correlation
794  G4double Phi = 0.;
795  if(std::abs(r*sth) < latcorr)
796  Phi = twopi*G4UniformRand();
797  else
798  {
799  G4double psi = std::acos(latcorr/(r*sth));
800  if(G4UniformRand() < 0.5)
801  Phi = phi+psi;
802  else
803  Phi = phi-psi;
804  }
805 
806  dirx = r*std::cos(Phi);
807  diry = r*std::sin(Phi);
808 
809  fDisplacement.set(dirx,diry,0.0);
810  fDisplacement.rotateUz(oldDirection);
811  }
812  }
813  return fDisplacement;
814 }
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:18
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)
void UrbanMscModel93::SetParticle ( const G4ParticleDefinition *  p)
inlineprivate

Definition at line 172 of file UrbanMscModel93.h.

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

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

173 {
174  if (p != particle) {
175  particle = p;
176  mass = p->GetPDGMass();
177  charge = p->GetPDGCharge()/CLHEP::eplus;
179  }
180 }
const G4ParticleDefinition * particle
G4double UrbanMscModel93::SimpleScattering ( G4double  xmeanth,
G4double  x2meanth 
)
inlineprivate

Definition at line 228 of file UrbanMscModel93.h.

References a, and TtFullHadEvtBuilder_cfi::prob.

Referenced by SampleCosineTheta().

229 {
230  // 'large angle scattering'
231  // 2 model functions with correct xmean and x2mean
232  G4double a = (2.*xmeanth+9.*x2meanth-3.)/(2.*xmeanth-3.*x2meanth+1.);
233  G4double prob = (a+2.)*xmeanth/a;
234 
235  // sampling
236  G4double cth = 1.;
237  if(G4UniformRand() < prob) {
238  cth = -1.+2.*G4Exp(G4Log(G4UniformRand())/(a+1.));
239  } else {
240  cth = -1.+2.*G4UniformRand();
241  }
242  return cth;
243 }
double a
Definition: hdecay.h:121
void UrbanMscModel93::StartTracking ( G4Track *  track)

Definition at line 377 of file UrbanMscModel93.cc.

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

378 {
379  SetParticle(track->GetDynamicParticle()->GetDefinition());
380  firstStep = true;
381  inside = false;
382  insideskin = false;
383  tlimit = geombig;
385  tlimitmin = 10.*stepmin ;
386 }
void SetParticle(const G4ParticleDefinition *)
void UrbanMscModel93::UpdateCache ( )
inlineprivate

Definition at line 185 of file UrbanMscModel93.h.

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

Referenced by SampleCosineTheta().

186 {
187  lnZ = G4Log(Zeff);
188  //new correction in theta0 formula
189  coeffth1 = (1.-8.7780e-2/Zeff)*(0.87+0.03*lnZ);
190  coeffth2 = (4.0780e-2+1.7315e-4*Zeff)*(0.87+0.03*lnZ);
191  // tail parameters
192  G4double lnZ1 = G4Log(Zeff+1.);
193  coeffc1 = 2.943-0.197*lnZ1;
194  coeffc2 = 0.0987-0.0143*lnZ1;
195  // for single scattering
196  Z2 = Zeff*Zeff;
197  Z23 = G4Exp(2.*lnZ/3.);
198  scr1 = scr1ini*Z23;
200 
201  Zold = Zeff;
202 }

Member Data Documentation

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

Definition at line 160 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::coeffc2
private

Definition at line 160 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::coeffth1
private

Definition at line 159 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::coeffth2
private

Definition at line 159 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 154 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::currentRadLength
private

Definition at line 147 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 149 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::eaa
private

Definition at line 149 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

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

Definition at line 110 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::fr
private

Definition at line 117 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::geombig
private

Definition at line 128 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::geomlimit
private

Definition at line 130 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::geommin
private

Definition at line 129 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4bool UrbanMscModel93::inside
private

Definition at line 164 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 117 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::lnZ
private

Definition at line 158 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

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

Definition at line 117 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::numlim
private

Definition at line 149 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

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

Definition at line 140 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 134 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::rangeinit
private

Definition at line 146 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::rellossmax
private

Definition at line 151 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::scr1
private

Definition at line 161 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::scr1ini
private

Definition at line 161 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::scr2
private

Definition at line 161 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::scr2ini
private

Definition at line 161 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::skindepth
private

Definition at line 131 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::smallstep
private

Definition at line 132 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 126 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4LossTableManager* UrbanMscModel93::theManager
private

Definition at line 113 of file UrbanMscModel93.h.

Referenced by UrbanMscModel93().

G4double UrbanMscModel93::theta0max
private

Definition at line 151 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

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

Definition at line 123 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::tlimitmin
private

Definition at line 124 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 149 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

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

Definition at line 158 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::Z23
private

Definition at line 158 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::Zeff
private

Definition at line 158 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::Zold
private

Definition at line 157 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::zPathLength
private