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