Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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
46G4PolarizedIonisationBhabhaXS::G4PolarizedIonisationBhabhaXS()
47 : fPhi0(1.)
48{
49 fPhi2 = G4ThreeVector();
50 fPhi3 = G4ThreeVector();
51}
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