| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // -------------------------------------------------------------------
27 //
28 // Geant4 Class file
29 //
30 // File name: G4PolarizedIonisationBhabhaXS
31 //
32 // Author: Andreas Schaelicke
33 //
34 // Class Description:
35 // * calculates the differential cross section
36 // incoming positron Kpl(along positive z direction) scatters at
37 // an electron Kmn at rest
38 // * phi denotes the angle between the scattering plane (defined by the
39 // outgoing electron) and X-axis
40 // * all stokes vectors refer to spins in the Global System (X,Y,Z)
41
42 #include "G4PolarizedIonisationBhabhaXS.hh"
43
44 #include "G4PhysicalConstants.hh"
45
46 G4PolarizedIonisationBhabhaXS::G4PolarizedIonisationBhabhaXS()
47 : fPhi0(1.)
48 {
49 fPhi2 = G4ThreeVector();
50 fPhi3 = G4ThreeVector();
51 }
52
53 G4PolarizedIonisationBhabhaXS::~G4PolarizedIonisationBhabhaXS() {}
54
55 void G4PolarizedIonisationBhabhaXS::Initialize(G4double e, G4double gamma,
56 G4double /*phi*/,
57 const G4StokesVector& pol0,
58 const G4StokesVector& pol1,
59 G4int flag)
60 {
61 SetXmax(1.);
62
63 constexpr G4double re2 = classic_electr_radius * classic_electr_radius;
64 G4double gamma2 = gamma * gamma;
65 G4double gamma3 = gamma2 * gamma;
66 G4double gmo = (gamma - 1.);
67 G4double gmo2 = (gamma - 1.) * (gamma - 1.);
68 G4double gmo3 = gmo2 * (gamma - 1.);
69 G4double gpo = (gamma + 1.);
70 G4double gpo2 = (gamma + 1.) * (gamma + 1.);
71 G4double gpo3 = gpo2 * (gamma + 1.);
72 G4double gpo12 = std::sqrt(gpo);
73 G4double gpo32 = gpo * gpo12;
74 G4double gpo52 = gpo2 * gpo12;
75
76 G4double pref = re2 / (gamma - 1.0);
77 constexpr G4double sqrttwo = 1.41421356237309504880; // sqrt(2.)
78 G4double d = std::sqrt(1. / e - 1.);
79 G4double e2 = e * e;
80 G4double e3 = e2 * e;
81
82 G4double gmo12 = std::sqrt(gmo);
83 G4double gmo32 = gmo * gmo12;
84 G4double egmp32 = std::pow(e * (2 + e * gmo) * gpo, (3. / 2.));
85 G4double e32 = e * std::sqrt(e);
86
87 G4bool polarized = (!pol0.IsZero()) || (!pol1.IsZero());
88
89 if(flag == 0)
90 polarized = false;
91 // Unpolarised part of XS
92 fPhi0 = e2 * gmo3 / gpo3;
93 fPhi0 -= 2. * e * gamma * gmo2 / gpo3;
94 fPhi0 += (3. * gamma2 + 6. * gamma + 4.) * gmo / gpo3;
95 fPhi0 -= (2. * gamma2 + 4. * gamma + 1.) / (e * gpo2);
96 fPhi0 += gamma2 / (e2 * (gamma2 - 1.));
97 fPhi0 *= 0.25;
98 // Initial state polarisation dependence
99 if(polarized)
100 {
101 G4double xx = -((e * gmo - gamma) *
102 (-1. - gamma + e * (e * gmo - gamma) * (3. + gamma))) /
103 (4. * e * gpo3);
104 G4double yy = (e3 * gmo3 - 2. * e2 * gmo2 * gamma -
105 gpo * (1. + 2. * gamma) + e * (-2. + gamma2 + gamma3)) /
106 (4. * e * gpo3);
107 G4double zz =
108 ((e * gmo - gamma) * (e2 * gmo * (3. + gamma) - e * gamma * (3. + gamma) +
109 gpo * (1. + 2. * gamma))) /
110 (4. * e * gpo3);
111
112 fPhi0 += xx * pol0.x() * pol1.x() + yy * pol0.y() * pol1.y() +
113 zz * pol0.z() * pol1.z();
114
115 {
116 G4double xy = 0.;
117 G4double xz = (d * (e * gmo - gamma) * (-1. + 2. * e * gmo - gamma)) /
118 (2. * sqrttwo * gpo52);
119 G4double yx = 0.;
120 G4double yz = 0.;
121 G4double zx = xz;
122 G4double zy = 0.;
123 fPhi0 += yx * pol0.y() * pol1.x() + xy * pol0.x() * pol1.y();
124 fPhi0 += zx * pol0.z() * pol1.x() + xz * pol0.x() * pol1.z();
125 fPhi0 += zy * pol0.z() * pol1.y() + yz * pol0.y() * pol1.z();
126 }
127 }
128 // Final state polarisarion dependence
129 fPhi2 = G4ThreeVector();
130 fPhi3 = G4ThreeVector();
131
132 if(flag >= 1)
133 {
134 // Final Positron Ppl
135 // initial positron Kpl
136 if(!pol0.IsZero())
137 {
138 G4double xxPplKpl = -((-1. + e) * (e * gmo - gamma) *
139 (-(gamma * gpo) + e * (-2. + gamma + gamma2))) /
140 (4. * e2 * gpo *
141 std::sqrt(gmo * gpo * (-1. + e + gamma - e * gamma) *
142 (1. + e + gamma - e * gamma)));
143 G4double xyPplKpl = 0.;
144 G4double xzPplKpl =
145 ((e * gmo - gamma) * (-1. - gamma + e * gmo * (1. + 2. * gamma))) /
146 (2. * sqrttwo * e32 * gmo * gpo2 *
147 std::sqrt(1. + e + gamma - e * gamma));
148 G4double yxPplKpl = 0.;
149 G4double yyPplKpl = (gamma2 * gpo + e2 * gmo2 * (3. + gamma) -
150 e * gmo * (1. + 2. * gamma * (2. + gamma))) /
151 (4. * e2 * gmo * gpo2);
152 G4double yzPplKpl = 0.;
153 G4double zxPplKpl =
154 ((e * gmo - gamma) *
155 (1. + e * (-1. + 2. * e * gmo - 2. * gamma) * gmo + gamma)) /
156 (2. * sqrttwo * e * gmo * gpo2 *
157 std::sqrt(e * (1. + e + gamma - e * gamma)));
158 G4double zyPplKpl = 0.;
159 G4double zzPplKpl =
160 -((e * gmo - gamma) * std::sqrt((1. - e) / (e - e * gamma2 + gpo2)) *
161 (2. * e2 * gmo2 + gamma + gamma2 - e * (-2. + gamma + gamma2))) /
162 (4. * e2 * (-1. + gamma2));
163
164 fPhi2[0] +=
165 xxPplKpl * pol0.x() + xyPplKpl * pol0.y() + xzPplKpl * pol0.z();
166 fPhi2[1] +=
167 yxPplKpl * pol0.x() + yyPplKpl * pol0.y() + yzPplKpl * pol0.z();
168 fPhi2[2] +=
169 zxPplKpl * pol0.x() + zyPplKpl * pol0.y() + zzPplKpl * pol0.z();
170 }
171 // initial electron Kmn
172 if(!pol1.IsZero())
173 {
174 G4double xxPplKmn =
175 ((-1. + e) * (e * (-2. + gamma) * gmo + gamma)) /
176 (4. * e * gpo32 * std::sqrt(1. + e2 * gmo + gamma - 2. * e * gamma));
177 G4double xyPplKmn = 0.;
178 G4double xzPplKmn =
179 (-1. + e * gmo + gmo * gamma) /
180 (2. * sqrttwo * gpo2 * std::sqrt(e * (1. + e + gamma - e * gamma)));
181 G4double yxPplKmn = 0.;
182 G4double yyPplKmn =
183 (-1. - 2. * gamma + e * gmo * (3. + gamma)) / (4. * e * gpo2);
184 G4double yzPplKmn = 0.;
185 G4double zxPplKmn =
186 (1. + 2. * e2 * gmo2 + gamma + gamma2 +
187 e * (1. + (3. - 4. * gamma) * gamma)) /
188 (2. * sqrttwo * gpo2 * std::sqrt(e * (1. + e + gamma - e * gamma)));
189 G4double zyPplKmn = 0.;
190 G4double zzPplKmn = -(std::sqrt((1. - e) / (e - e * gamma2 + gpo2)) *
191 (2. * e2 * gmo2 + gamma + 2. * gamma2 +
192 e * (2. + gamma - 3. * gamma2))) /
193 (4. * e * gpo);
194
195 fPhi2[0] +=
196 xxPplKmn * pol1.x() + xyPplKmn * pol1.y() + xzPplKmn * pol1.z();
197 fPhi2[1] +=
198 yxPplKmn * pol1.x() + yyPplKmn * pol1.y() + yzPplKmn * pol1.z();
199 fPhi2[2] +=
200 zxPplKmn * pol1.x() + zyPplKmn * pol1.y() + zzPplKmn * pol1.z();
201 }
202 // Final Electron Pmn
203 // initial positron Kpl
204 if(!pol0.IsZero())
205 {
206 G4double xxPmnKpl = ((-1. + e * gmo) * (2. + gamma)) /
207 (4. * gpo * std::sqrt(e * (2. + e * gmo) * gpo));
208 G4double xyPmnKpl = 0.;
209 G4double xzPmnKpl = (std::sqrt((-1. + e) / (-2. + e - e * gamma)) *
210 (e + gamma + e * gamma - 2. * (-1. + e) * gamma2)) /
211 (2. * sqrttwo * e * gpo2);
212 G4double yxPmnKpl = 0.;
213 G4double yyPmnKpl =
214 (-1. - 2. * gamma + e * gmo * (3. + gamma)) / (4. * e * gpo2);
215 G4double yzPmnKpl = 0.;
216 G4double zxPmnKpl =
217 -((-1. + e) * (1. + 2. * e * gmo) * (e * gmo - gamma)) /
218 (2. * sqrttwo * e * std::sqrt(-((-1. + e) * (2. + e * gmo))) * gpo2);
219 G4double zyPmnKpl = 0;
220 G4double zzPmnKpl = (-2. + 2. * e2 * gmo2 + gamma * (-1. + 2. * gamma) +
221 e * (-2. + (5. - 3. * gamma) * gamma)) /
222 (4. * std::sqrt(e * (2. + e * gmo)) * gpo32);
223
224 fPhi3[0] +=
225 xxPmnKpl * pol0.x() + xyPmnKpl * pol0.y() + xzPmnKpl * pol0.z();
226 fPhi3[1] +=
227 yxPmnKpl * pol0.x() + yyPmnKpl * pol0.y() + yzPmnKpl * pol0.z();
228 fPhi3[2] +=
229 zxPmnKpl * pol0.x() + zyPmnKpl * pol0.y() + zzPmnKpl * pol0.z();
230 }
231 // initial electron Kmn
232 if(!pol1.IsZero())
233 {
234 G4double xxPmnKmn = -((2. + e * gmo) * (-1. + e * gmo - gamma) *
235 (e * gmo - gamma) * (-2. + gamma)) /
236 (4. * gmo * egmp32);
237 G4double xyPmnKmn = 0.;
238 G4double xzPmnKmn =
239 ((e * gmo - gamma) *
240 std::sqrt((-1. + e + gamma - e * gamma) / (2. + e * gmo)) *
241 (e + gamma - e * gamma + gamma2)) /
242 (2. * sqrttwo * e2 * gmo32 * gpo2);
243 G4double yxPmnKmn = 0.;
244 G4double yyPmnKmn = (gamma2 * gpo + e2 * gmo2 * (3. + gamma) -
245 e * gmo * (1. + 2. * gamma * (2. + gamma))) /
246 (4. * e2 * gmo * gpo2);
247 G4double yzPmnKmn = 0.;
248 G4double zxPmnKmn =
249 -((-1. + e) * (e * gmo - gamma) *
250 (e * gmo + 2. * e2 * gmo2 - gamma * gpo)) /
251 (2. * sqrttwo * e2 * std::sqrt(-((-1. + e) * (2. + e * gmo))) * gmo *
252 gpo2);
253 G4double zyPmnKmn = 0.;
254 G4double zzPmnKmn =
255 ((e * gmo - gamma) * std::sqrt(e / ((2. + e * gmo) * gpo)) *
256 (-(e * (-2. + gamma) * gmo) + 2. * e2 * gmo2 + (-2. + gamma) * gpo)) /
257 (4. * e2 * (-1. + gamma2));
258
259 fPhi3[0] +=
260 xxPmnKmn * pol1.x() + xyPmnKmn * pol1.y() + xzPmnKmn * pol1.z();
261 fPhi3[1] +=
262 yxPmnKmn * pol1.x() + yyPmnKmn * pol1.y() + yzPmnKmn * pol1.z();
263 fPhi3[2] +=
264 zxPmnKmn * pol1.x() + zyPmnKmn * pol1.y() + zzPmnKmn * pol1.z();
265 }
266 }
267 fPhi0 *= pref;
268 fPhi2 *= pref;
269 fPhi3 *= pref;
270 }
271
272 G4double G4PolarizedIonisationBhabhaXS::XSection(const G4StokesVector& pol2,
273 const G4StokesVector& pol3)
274 {
275 G4double xs = 0.;
276 xs += fPhi0;
277
278 G4bool polarized = (!pol2.IsZero()) || (!pol3.IsZero());
279 if(polarized)
280 {
281 xs += fPhi2 * pol2 + fPhi3 * pol3;
282 }
283 return xs;
284 }
285
286 G4double G4PolarizedIonisationBhabhaXS::TotalXSection(
287 G4double xmin, G4double /*xmax*/, G4double gamma,
288 const G4StokesVector& pol0,
289 const G4StokesVector& pol1)
290 {
291 // VI: In this model an incorrect G4Exception was used in which
292 // xmax was compared with 1. In the case of electron this
293 // value is 0.5, for positrons 1.0. The computation of the
294 // cross section is left unchanged, so part of integral
295 // from 0.5 to 1.0 is computed for electrons, which is not
296 // correct, however, this part is significantly smaller then
297 // from the interval xmin - 0.5.
298 G4double xs = 0.;
299 G4double x = xmin;
300
301 constexpr G4double re2 = classic_electr_radius * classic_electr_radius;
302 G4double gamma2 = gamma * gamma;
303 G4double gmo2 = (gamma - 1.) * (gamma - 1.);
304 G4double gpo2 = (gamma + 1.) * (gamma + 1.);
305 G4double gpo3 = gpo2 * (gamma + 1.);
306 G4double logMEM = std::log(x);
307 G4double pref = twopi * re2 / (gamma - 1.0);
308 // unpolarised XS
309 G4double sigma0 = 0.;
310 sigma0 += -gmo2 * (gamma - 1.) * x * x * x / 3. + gmo2 * gamma * x * x;
311 sigma0 += -(gamma - 1.) * (3. * gamma * (gamma + 2.) + 4.) * x;
312 sigma0 += (gamma * (gamma * (gamma * (4. * gamma - 1.) - 21.) - 7.) + 13.) /
313 (3. * (gamma - 1.));
314 sigma0 /= gpo3;
315 sigma0 += logMEM * (2. - 1. / gpo2);
316 sigma0 += gamma2 / ((gamma2 - 1.) * x);
317 // longitudinal part
318 G4double sigma2 = 0.;
319 sigma2 += logMEM * gamma * (gamma + 1.) * (2. * gamma + 1.);
320 sigma2 += gamma * (7. * gamma * (gamma + 1.) - 2.) / 3.;
321 sigma2 += -(3. * gamma + 1.) * (gamma2 + gamma - 1.) * x;
322 sigma2 += (gamma - 1.) * gamma * (gamma + 3.) * x * x;
323 sigma2 += -gmo2 * (gamma + 3.) * x * x * x / 3.;
324 sigma2 /= gpo3;
325 // transverse part
326 G4double sigma3 = 0.;
327 sigma3 += 0.5 * (gamma + 1.) * (3. * gamma + 1.) * logMEM;
328 sigma3 += (gamma * (5. * gamma - 4.) - 13.) / 6.;
329 sigma3 += 0.5 * (gamma2 + 3.) * x;
330 sigma3 += -2. * (gamma - 1.) * gamma * x * x;
331 sigma3 += 2. * gmo2 * x * x * x / 3.;
332 sigma3 /= gpo3;
333 // total cross section
334 xs += pref * (sigma0 + sigma2 * pol0.z() * pol1.z() +
335 sigma3 * (pol0.x() * pol1.x() + pol0.y() * pol1.y()));
336
337 return xs;
338 }
339
340 G4StokesVector G4PolarizedIonisationBhabhaXS::GetPol2()
341 {
342 // Note, mean polarization can not contain correlation effects.
343 return G4StokesVector(1. / fPhi0 * fPhi2);
344 }
345 G4StokesVector G4PolarizedIonisationBhabhaXS::GetPol3()
346 {
347 // Note, mean polarization can not contain correlation effects.
348 return G4StokesVector(1. / fPhi0 * fPhi3);
349 }
350