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