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) override
 
G4double ComputeGeomPathLength (G4double truePathLength) override
 
G4double ComputeTheta0 (G4double truePathLength, G4double KineticEnergy)
 
G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep) override
 
G4double ComputeTrueStepLength (G4double geomStepLength) override
 
void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
G4ThreeVector & SampleScattering (const G4ThreeVector &, G4double safety) override
 
void StartTracking (G4Track *) override
 
 UrbanMscModel93 (const G4String &nam="UrbanMsc93")
 
 ~UrbanMscModel93 () override
 

Private Member Functions

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

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 39 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  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;
47  tlimitminfix = 1.e-6 * mm;
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;
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;
96 
98  fParticleChange = nullptr;
99  couple = nullptr;
100  SetSampleZ(false);
101 }
const G4ParticleDefinition * particle
const G4MaterialCutsCouple * couple
G4ParticleChangeForMSC * fParticleChange
G4double currentRadLength
const Double_t pi
const double MeV
G4LossTableManager * theManager
G4double currentKinEnergy
UrbanMscModel93::~UrbanMscModel93 ( )
override

Definition at line 105 of file UrbanMscModel93.cc.

105 {}
UrbanMscModel93::UrbanMscModel93 ( const UrbanMscModel93 )
privatedelete

Member Function Documentation

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

Definition at line 125 of file UrbanMscModel93.cc.

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

126  {
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 }
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:19
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:30
G4double UrbanMscModel93::ComputeGeomPathLength ( G4double  truePathLength)
override

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

556  {
557  firstStep = false;
558  lambdaeff = lambda0;
559  par1 = -1.;
560  par2 = par3 = 0.;
561 
562  // do the true -> geom transformation
564 
565  // z = t for very small tPathLength
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
573 
574  G4double tau = tPathLength / lambda0;
575 
576  if ((tau <= tausmall) || insideskin) {
578  if (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;
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)
635  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda0 << G4endl;
636  return zPathLength;
637 }
const G4ParticleDefinition * particle
const G4MaterialCutsCouple * couple
G4double currentKinEnergy
G4double UrbanMscModel93::ComputeTheta0 ( G4double  truePathLength,
G4double  KineticEnergy 
)
inline

Definition at line 196 of file UrbanMscModel93.h.

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

Referenced by SampleCosineTheta().

196  {
197  // for all particles take the width of the central part
198  // from a parametrization similar to the Highland formula
199  // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
200  static const G4double c_highland = 13.6 * CLHEP::MeV;
201  G4double invbetacp =
202  sqrt((currentKinEnergy + mass) * (KineticEnergy + mass) /
203  (currentKinEnergy * (currentKinEnergy + 2. * mass) * KineticEnergy * (KineticEnergy + 2. * mass)));
204  y = trueStepLength / currentRadLength;
205  G4double theta0 = c_highland * std::abs(charge) * sqrt(y) * invbetacp;
206  // correction factor from e- scattering data
207  theta0 *= (coeffth1 + coeffth2 * G4Log(y));
208 
209  return theta0;
210 }
G4double currentRadLength
const double MeV
T sqrt(T t)
Definition: SSEVec.h:19
G4double currentKinEnergy
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
G4double UrbanMscModel93::ComputeTruePathLengthLimit ( const G4Track &  track,
G4double &  currentMinimalStep 
)
override

Definition at line 349 of file UrbanMscModel93.cc.

References couple, currentKinEnergy, currentMaterialIndex, currentRange, Calorimetry_cff::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.

349  {
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();
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) {
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) {
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)
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)
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)
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)
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)
547  }
548  }
549  //G4cout << "tPathLength= " << tPathLength
550  // << " currentMinimalStep= " << currentMinimalStep << G4endl;
551  return ConvertTrueToGeom(tPathLength, currentMinimalStep);
552 }
const G4ParticleDefinition * particle
const G4MaterialCutsCouple * couple
const double MeV
G4double currentKinEnergy
G4double UrbanMscModel93::ComputeTrueStepLength ( G4double  geomStepLength)
override

Definition at line 641 of file UrbanMscModel93.cc.

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

641  {
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
661  }
662  }
663  if (tPathLength < geomStepLength)
664  tPathLength = geomStepLength;
665 
666  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
667 
668  return tPathLength;
669 }
void UrbanMscModel93::Initialise ( const G4ParticleDefinition *  p,
const G4DataVector &   
)
override

Definition at line 109 of file UrbanMscModel93.cc.

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

109  {
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 }
G4ParticleChangeForMSC * fParticleChange
void SetParticle(const G4ParticleDefinition *)
const double MeV
G4double UrbanMscModel93::LatCorrelation ( )
private

Definition at line 1001 of file UrbanMscModel93.cc.

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

Referenced by SampleScattering().

1001  {
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 }
static const G4double kappami1
static const G4double kappapl1
static const G4double kappa
UrbanMscModel93& UrbanMscModel93::operator= ( const UrbanMscModel93 right)
privatedelete
G4double UrbanMscModel93::SampleCosineTheta ( G4double  trueStepLength,
G4double  KineticEnergy 
)
private

Definition at line 763 of file UrbanMscModel93.cc.

References b, l1GtPatternGenerator_cfi::bx, HltBtagPostValidation_cff::c, alignmentValidation::c1, coeffc1, coeffc2, ComputeTheta0(), funct::cos(), couple, currentKinEnergy, currentRadLength, currentRange, currentTau, ztail::d, ea, eaa, ecalTB2006H4_GenSimDigiReco_cfg::G4cout, GeV, mps_fire::i, insideskin, lambda0, lambdaeff, mass, SiStripPI::mean, dqmiodumpmetadata::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(), trigObjTnPSource_cfi::var, x, xsi, y, Zeff, and Zold.

Referenced by SampleScattering().

763  {
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 }
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:19
G4double currentKinEnergy
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
G4double SimpleScattering(G4double xmeanth, G4double x2meanth)
d
Definition: ztail.py:151
double b
Definition: hdecay.h:118
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:30
G4double UrbanMscModel93::SampleDisplacement ( )
private

Definition at line 948 of file UrbanMscModel93.cc.

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

Referenced by SampleScattering().

948  {
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 }
static const G4double kappami1
T sqrt(T t)
Definition: SSEVec.h:19
static const G4double kappapl1
static const G4double kappa
G4ThreeVector & UrbanMscModel93::SampleScattering ( const G4ThreeVector &  oldDirection,
G4double  safety 
)
override

Definition at line 673 of file UrbanMscModel93.cc.

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

673  {
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 }
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:19
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 165 of file UrbanMscModel93.h.

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

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

165  {
166  if (p != particle) {
167  particle = p;
168  mass = p->GetPDGMass();
169  charge = p->GetPDGCharge() / CLHEP::eplus;
171  }
172 }
const G4ParticleDefinition * particle
G4double UrbanMscModel93::SimpleScattering ( G4double  xmeanth,
G4double  x2meanth 
)
inlineprivate

Definition at line 214 of file UrbanMscModel93.h.

References a, and TtFullHadEvtBuilder_cfi::prob.

Referenced by SampleCosineTheta().

214  {
215  // 'large angle scattering'
216  // 2 model functions with correct xmean and x2mean
217  G4double a = (2. * xmeanth + 9. * x2meanth - 3.) / (2. * xmeanth - 3. * x2meanth + 1.);
218  G4double prob = (a + 2.) * xmeanth / a;
219 
220  // sampling
221  G4double cth = 1.;
222  if (G4UniformRand() < prob) {
223  cth = -1. + 2. * G4Exp(G4Log(G4UniformRand()) / (a + 1.));
224  } else {
225  cth = -1. + 2. * G4UniformRand();
226  }
227  return cth;
228 }
double a
Definition: hdecay.h:119
void UrbanMscModel93::StartTracking ( G4Track *  track)
override

Definition at line 337 of file UrbanMscModel93.cc.

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

337  {
338  SetParticle(track->GetDynamicParticle()->GetDefinition());
339  firstStep = true;
340  inside = false;
341  insideskin = false;
342  tlimit = geombig;
344  tlimitmin = 10. * stepmin;
345 }
void SetParticle(const G4ParticleDefinition *)
void UrbanMscModel93::UpdateCache ( )
inlineprivate

Definition at line 176 of file UrbanMscModel93.h.

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

Referenced by SampleCosineTheta().

176  {
177  lnZ = G4Log(Zeff);
178  //new correction in theta0 formula
179  coeffth1 = (1. - 8.7780e-2 / Zeff) * (0.87 + 0.03 * lnZ);
180  coeffth2 = (4.0780e-2 + 1.7315e-4 * Zeff) * (0.87 + 0.03 * lnZ);
181  // tail parameters
182  G4double lnZ1 = G4Log(Zeff + 1.);
183  coeffc1 = 2.943 - 0.197 * lnZ1;
184  coeffc2 = 0.0987 - 0.0143 * lnZ1;
185  // for single scattering
186  Z2 = Zeff * Zeff;
187  Z23 = G4Exp(2. * lnZ / 3.);
188  scr1 = scr1ini * Z23;
189  scr2 = scr2ini * Z2 * ChargeSquare;
190 
191  Zold = Zeff;
192 }

Member Data Documentation

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

Definition at line 154 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::coeffc2
private

Definition at line 154 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::coeffth1
private

Definition at line 153 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::coeffth2
private

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

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::currentRadLength
private

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

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::eaa
private

Definition at line 143 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

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

Definition at line 104 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::fr
private

Definition at line 111 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::geombig
private

Definition at line 122 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::geomlimit
private

Definition at line 124 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::geommin
private

Definition at line 123 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4bool UrbanMscModel93::inside
private

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

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::lnZ
private

Definition at line 152 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

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

Definition at line 111 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::numlim
private

Definition at line 143 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

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

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

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::rangeinit
private

Definition at line 140 of file UrbanMscModel93.h.

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4double UrbanMscModel93::rellossmax
private

Definition at line 145 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

G4double UrbanMscModel93::scr1
private

Definition at line 155 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::scr1ini
private

Definition at line 155 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::scr2
private

Definition at line 155 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::scr2ini
private

Definition at line 155 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::skindepth
private

Definition at line 125 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::smallstep
private

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

Referenced by ComputeTruePathLengthLimit(), and UrbanMscModel93().

G4LossTableManager* UrbanMscModel93::theManager
private

Definition at line 107 of file UrbanMscModel93.h.

Referenced by UrbanMscModel93().

G4double UrbanMscModel93::theta0max
private

Definition at line 145 of file UrbanMscModel93.h.

Referenced by SampleCosineTheta(), and UrbanMscModel93().

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

Definition at line 117 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::tlimitmin
private

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

Referenced by SampleCosineTheta(), and UrbanMscModel93().

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

Definition at line 152 of file UrbanMscModel93.h.

Referenced by UpdateCache(), and UrbanMscModel93().

G4double UrbanMscModel93::Z23
private

Definition at line 152 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::Zeff
private

Definition at line 152 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::Zold
private

Definition at line 151 of file UrbanMscModel93.h.

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

G4double UrbanMscModel93::zPathLength
private