21 #include "CLHEP/Units/PhysicalConstants.h"
22 #include "Randomize.hh"
23 #include "G4Electron.hh"
24 #include "G4LossTableManager.hh"
25 #include "G4ParticleChangeForMSC.hh"
27 #include "G4Poisson.hh"
33 using namespace CLHEP;
35 static const G4double
kappa = 2.5;
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;
93 mass = proton_mass_c2;
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;
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;
132 static const G4double Zdat[15] = {4., 6., 13., 20., 26., 29., 32., 38., 47., 50., 56., 64., 74., 79., 82.};
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,
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}};
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}};
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);
209 static const G4double sig0[15] = {0.2672 * barn,
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};
231 Z23 =
pow(AtomicNumber, 2. / 3.);
237 G4double eKineticEnergy = KineticEnergy;
239 if (
mass > electron_mass_c2) {
240 G4double
TAU = KineticEnergy /
mass;
241 G4double
c =
mass *
TAU * (
TAU + 2.) / (electron_mass_c2 * (
TAU + 1.));
244 eKineticEnergy = electron_mass_c2 *
tau;
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);
251 G4double
eps = epsfactor * bg2 /
Z23;
255 else if (
eps < epsmax)
256 sigma = G4Log(1. + 2. *
eps) - 2. *
eps / (1. + 2. *
eps);
258 sigma = G4Log(2. *
eps) - 1. + 1. /
eps;
260 sigma *=
ChargeSquare * AtomicNumber * AtomicNumber / (beta2 * bg2);
263 G4double
c1, c2, cc1, cc2,
corr;
267 while ((iZ >= 0) && (Zdat[iZ] >= AtomicNumber))
274 G4double ZZ1 = Zdat[iZ];
275 G4double ZZ2 = Zdat[iZ + 1];
276 G4double ratZ = (AtomicNumber - ZZ1) * (AtomicNumber + ZZ1) / ((ZZ2 - ZZ1) * (ZZ2 + ZZ1));
278 if (eKineticEnergy <= Tlim) {
281 while ((iT >= 0) && (Tdat[iT] >= eKineticEnergy))
289 G4double
T = Tdat[iT], E =
T + electron_mass_c2;
290 G4double b2small =
T * (E + electron_mass_c2) / (E * E);
293 E =
T + electron_mass_c2;
294 G4double b2big =
T * (E + electron_mass_c2) / (E * E);
295 G4double ratb2 = (beta2 - b2small) / (b2big - b2small);
298 c1 = celectron[iZ][iT];
299 c2 = celectron[iZ + 1][iT];
300 cc1 =
c1 + ratZ * (c2 -
c1);
302 c1 = celectron[iZ][iT + 1];
303 c2 = celectron[iZ + 1][iT + 1];
304 cc2 =
c1 + ratZ * (c2 -
c1);
306 corr = cc1 + ratb2 * (cc2 - cc1);
308 sigma *= sigmafactor /
corr;
310 c1 = cpositron[iZ][iT];
311 c2 = cpositron[iZ + 1][iT];
312 cc1 =
c1 + ratZ * (c2 -
c1);
314 c1 = cpositron[iZ][iT + 1];
315 c2 = cpositron[iZ + 1][iT + 1];
316 cc2 =
c1 + ratZ * (c2 -
c1);
318 corr = cc1 + ratb2 * (cc2 - cc1);
320 sigma *= sigmafactor /
corr;
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);
351 const G4DynamicParticle*
dp =
track.GetDynamicParticle();
352 G4StepPoint* sp =
track.GetStep()->GetPreStepPoint();
353 G4StepStatus stepStatus = sp->GetStepStatus();
363 return ConvertTrueToGeom(
tPathLength, currentMinimalStep);
380 return ConvertTrueToGeom(
tPathLength, currentMinimalStep);
385 if (steppingAlgorithm == fUseDistanceToBoundary) {
392 return ConvertTrueToGeom(
tPathLength, currentMinimalStep);
398 if (
firstStep || stepStatus == fGeomBoundary) {
408 rat = 1.e-3 / (rat * (10. + rat));
425 if (stepStatus == fGeomBoundary)
451 return ConvertTrueToGeom(
tPathLength, currentMinimalStep);
473 G4double temptlimit =
tlimit;
490 else if (steppingAlgorithm == fUseSafety) {
499 return ConvertTrueToGeom(
tPathLength, currentMinimalStep);
502 if (
firstStep || stepStatus == fGeomBoundary) {
515 rat = 1.e-3 / (rat * (10. + rat));
537 if (stepStatus == fGeomBoundary) {
551 return ConvertTrueToGeom(
tPathLength, currentMinimalStep);
599 G4double lambda1 = GetTransportMeanFreePath(
particle, T1);
611 const G4double ztmax = 0.99;
617 G4double cz = 0.5 * (3. *
zt - 1.) / (1. -
zt);
619 G4double u0 = cz / cz1;
622 u = G4Exp(G4Log(G4UniformRand()) / cz1);
623 grej = G4Exp(cz * G4Log(u / u0)) * (1. - u) / (1. - u0);
624 }
while (grej < G4UniformRand());
627 u = 1. - G4Exp(G4Log(G4UniformRand()) / cz1);
657 if (
par1 *
par3 * geomStepLength < 1.)
674 fDisplacement.set(0.0, 0.0, 0.0);
682 return fDisplacement;
688 if (std::fabs(cth) > 1.) {
689 return fDisplacement;
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);
720 G4ThreeVector newDirection(dirx, diry, cth);
721 newDirection.rotateUz(oldDirection);
742 Phi = twopi * G4UniformRand();
744 G4double
psi = std::acos(latcorr / (
r * sth));
745 if (G4UniformRand() < 0.5)
754 fDisplacement.set(dirx, diry, 0.0);
755 fDisplacement.rotateUz(oldDirection);
758 return fDisplacement;
767 Zeff =
couple->GetMaterial()->GetTotNbOfElectPerVolume() /
couple->GetMaterial()->GetTotNbOfAtomsPerVolume();
776 G4int
n = G4Poisson(
mean);
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.;
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);
796 st =
sqrt(1. - ct * ct);
797 phi = twopi * G4UniformRand();
806 if (
par1 * trueStepLength < 1.)
818 cth = -1. + 2. * G4UniformRand();
820 G4double xmeanth, x2meanth;
822 xmeanth = 1.0 -
tau * (1.0 - 0.5 *
tau);
825 xmeanth = G4Exp(-
tau);
826 x2meanth = (1. + 2. * G4Exp(-2.5 *
tau)) *
third;
833 G4double theta0 =
ComputeTheta0(trueStepLength, KineticEnergy);
839 G4double theta2 = theta0 * theta0;
849 G4double
x = theta2 * (1.0 - theta2 / 12.);
851 G4double sth = 2. *
sin(0.5 * theta0);
855 G4double xmean1 = 1. - (1. - (1. +
xsi) *
ea) *
x /
eaa;
856 G4double x0 = 1. -
xsi *
x;
860 if (xmean1 <= 0.999 * xmeanth) {
881 if (fabs(
c - 3.) < 0.001) {
883 }
else if (fabs(
c - 2.) < 0.001) {
887 G4double
c1 =
c - 1.;
890 G4double
b = 1. + (
c -
xsi) *
x;
892 G4double
b1 =
b + 1.;
897 G4double
d = ebx / eb1;
900 G4double xmean2 = (x0 +
d - (
bx -
b1 *
d) / (
c - 2.)) / (1. -
d);
902 G4double f1x0 =
ea /
eaa;
903 G4double f2x0 =
c1 / (
c * (1. -
d));
904 G4double
prob = f2x0 / (f1x0 + f2x0);
906 G4double qprob = xmeanth / (
prob * xmean1 + (1. -
prob) * xmean2);
912 if (G4UniformRand() < qprob) {
914 if (G4UniformRand() <
prob) {
915 cth = 1. + G4Log(
ea + G4UniformRand() *
eaa) *
x;
917 var = (1.0 -
d) * G4UniformRand();
920 cth = -1.0 +
var * (1.0 - 0.5 *
var *
c) * (2. + (
c -
xsi) *
x);
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
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;
950 G4double rmean = 0.0;
1002 G4double latcorr = 0.;
1011 latcorr = G4Exp(latcorr) /
kappami1;