Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4PolarizedIonisationBhabhaXS.cc
Go to the documentation of this file.
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
52
53G4PolarizedIonisationBhabhaXS::~G4PolarizedIonisationBhabhaXS() {}
54
55void 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
272G4double 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
286G4double 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
340G4StokesVector G4PolarizedIonisationBhabhaXS::GetPol2()
341{
342 // Note, mean polarization can not contain correlation effects.
343 return G4StokesVector(1. / fPhi0 * fPhi2);
344}
345G4StokesVector G4PolarizedIonisationBhabhaXS::GetPol3()
346{
347 // Note, mean polarization can not contain correlation effects.
348 return G4StokesVector(1. / fPhi0 * fPhi3);
349}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
G4PolarizedIonisationBhabhaXS::TotalXSection
G4double TotalXSection(G4double xmin, G4double xmax, G4double y, const G4StokesVector &pol0, const G4StokesVector &pol1) override
Definition G4PolarizedIonisationBhabhaXS.cc:286
G4PolarizedIonisationBhabhaXS::Initialize
void Initialize(G4double x, G4double y, G4double phi, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0) override
Definition G4PolarizedIonisationBhabhaXS.cc:55
G4PolarizedIonisationBhabhaXS::XSection
G4double XSection(const G4StokesVector &pol2, const G4StokesVector &pol3) override
Definition G4PolarizedIonisationBhabhaXS.cc:272
G4StokesVector::IsZero
G4bool IsZero() const
Definition G4StokesVector.cc:65
G4VPolarizedXS::SetXmax
void SetXmax(G4double xmax)
Definition G4VPolarizedXS.hh:84