Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4PolarizedBhabhaCrossSection.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// $Id$
27// -------------------------------------------------------------------
28//
29// GEANT4 Class file
30//
31//
32// File name: G4PolarizedBhabhaCrossSection
33//
34// Author: Andreas Schaelicke
35//
36// Creation date: 12.01.2006
37//
38// Modifications:
39// 16-01-06 included cross section as calculated by P.Starovoitov
40// 24-08-06 bugfix in total cross section (A. Schaelicke)
41// 07-11-06 modify reference system for polarisation vectors
42// (A. Schaelicke & P.Starovoitov)
43//
44// Class Description:
45// * calculates the differential cross section
46// incomming positron Kpl(along positive z direction) scatters at
47// an electron Kmn at rest
48// * phi denotes the angle between the scattering plane (defined by the
49// outgoing electron) and X-axis
50// * all stokes vectors refer to spins in the Global System (X,Y,Z)
51//
52
53#include "G4PolarizedBhabhaCrossSection.hh"
54#include "G4PhysicalConstants.hh"
55
56G4PolarizedBhabhaCrossSection::G4PolarizedBhabhaCrossSection() : phi0(1.)
57{
58}
59G4PolarizedBhabhaCrossSection::~G4PolarizedBhabhaCrossSection()
60{
61}
62void G4PolarizedBhabhaCrossSection::Initialize(
63 G4double e,
64 G4double gamma,
65 G4double /*phi*/,
66 const G4StokesVector & pol0,
67 const G4StokesVector & pol1,
68 G4int flag)
69{
70 SetXmax(1.);
71
72 G4double re2 = classic_electr_radius * classic_electr_radius;
73 G4double gamma2 = gamma*gamma;
74 G4double gamma3 = gamma2*gamma;
75 G4double gmo = (gamma - 1.);
76 G4double gmo2 = (gamma - 1.)*(gamma - 1.);
77 G4double gmo3 = gmo2*(gamma - 1.);
78 G4double gpo = (gamma + 1.);
79 G4double gpo2 = (gamma + 1.)*(gamma + 1.);
80 G4double gpo3 = gpo2*(gamma + 1.);
81 G4double gpo12 = std::sqrt(gpo);
82 G4double gpo32 = gpo*gpo12;
83 G4double gpo52 = gpo2*gpo12;
84
85 G4double pref = re2/(gamma - 1.0);
86 G4double sqrttwo=std::sqrt(2.);
87 G4double d = std::sqrt(1./e - 1.);
88 G4double e2 = e*e;
89 G4double e3 = e2*e;
90
91 // *** new ***
92 G4double gmo12 = std::sqrt(gmo);
93 G4double gmo32 = gmo*gmo12;
94 G4double egmp32 = std::pow(e*(2 + e*gmo)*gpo,(3./2.));
95 G4double e32 = e*std::sqrt(e);
96
97 G4bool polarized=(!pol0.IsZero())||(!pol1.IsZero());
98
99 if (flag==0) polarized=false;
100 // Unpolarised part of XS
101 // *AS* UnpME . OK
102 phi0 = 0.;
103 phi0+= e2*gmo3/gpo3;
104 phi0+= -2.*e*gamma*gmo2/gpo3;
105 phi0+= (3.*gamma2 + 6.*gamma + 4.)*gmo/gpo3;
106 phi0+= -(2.*gamma2 + 4.*gamma + 1.)/(e*gpo2);
107 phi0+= gamma2/(e2*(gamma2 - 1.));
108 phi0*=0.25;
109 // Initial state polarisarion dependence
110 if (polarized) {
111 // G4cout<<"Polarized differential Bhabha cross section"<<G4endl;
112 // G4cout<<"Initial state polarisation contributions"<<G4endl;
113 // G4cout<<"Diagonal Matrix Elements"<<G4endl;
114 // *** new ***
115 G4double xx = -((e*gmo - gamma)*(-1 - gamma + e*(e*gmo - gamma)*(3 + gamma)))/(4*e*gpo3);
116 G4double yy = (e3*gmo3 - 2*e2*gmo2*gamma - gpo*(1 + 2*gamma) + e*(-2 + gamma2 + gamma3))/(4*e*gpo3);
117 G4double zz = ((e*gmo - gamma)*(e2*gmo*(3 + gamma) - e*gamma*(3 + gamma) + gpo*(1 + 2*gamma)))/(4*e*gpo3);
118 // ***
119
120 phi0 += xx*pol0.x()*pol1.x() + yy*pol0.y()*pol1.y() + zz*pol0.z()*pol1.z();
121
122 {
123 G4double xy = 0;
124 G4double xz = (d*(e*gmo - gamma)*(-1 + 2*e*gmo - gamma))/(2*sqrttwo*gpo52);
125 G4double yx = 0;
126 G4double yz = 0;
127 G4double zx = xz;
128 G4double zy = 0;
129 // G4cout<<"Non-diagonal Matrix Elements"<<G4endl;
130 phi0+=yx*pol0.y()*pol1.x() + xy*pol0.x()*pol1.y();
131 phi0+=zx*pol0.z()*pol1.x() + xz*pol0.x()*pol1.z();
132 phi0+=zy*pol0.z()*pol1.y() + yz*pol0.y()*pol1.z();
133 }
134 }
135 // Final state polarisarion dependence
136 phi2=G4ThreeVector();
137 phi3=G4ThreeVector();
138
139 if (flag>=1) {
140 //
141 // Final Positron Ppl
142 //
143 // initial positron Kpl
144 if (!pol0.IsZero()) {
145
146 G4double xxPplKpl = -((-1 + e)*(e*gmo - gamma)*(-(gamma*gpo) + e*(-2 + gamma + gamma2)))/
147 (4*e2*gpo*std::sqrt(gmo*gpo*(-1 + e + gamma - e*gamma)* (1 + e + gamma - e*gamma)));
148 G4double xyPplKpl = 0;
149 G4double xzPplKpl = ((e*gmo - gamma)*(-1 - gamma + e*gmo*(1 + 2*gamma)))/
150 (2*sqrttwo*e32*gmo*gpo2*std::sqrt(1 + e + gamma - e*gamma));
151 G4double yxPplKpl = 0;
152 G4double yyPplKpl = (gamma2*gpo + e2*gmo2*(3 + gamma) -
153 e*gmo*(1 + 2*gamma*(2 + gamma)))/(4*e2*gmo*gpo2);
154 G4double yzPplKpl = 0;
155 G4double zxPplKpl = ((e*gmo - gamma)*(1 + e*(-1 + 2*e*gmo - 2*gamma)*gmo + gamma))/
156 (2*sqrttwo*e*gmo*gpo2*std::sqrt(e*(1 + e + gamma - e*gamma)));
157 G4double zyPplKpl = 0;
158 G4double zzPplKpl = -((e*gmo - gamma)*std::sqrt((1 - e)/(e - e*gamma2 + gpo2))*
159 (2*e2*gmo2 + gamma + gamma2 - e*(-2 + gamma + gamma2)))/
160 (4*e2*(-1 + gamma2));
161
162 phi2[0] += xxPplKpl*pol0.x() + xyPplKpl*pol0.y() + xzPplKpl*pol0.z();
163 phi2[1] += yxPplKpl*pol0.x() + yyPplKpl*pol0.y() + yzPplKpl*pol0.z();
164 phi2[2] += zxPplKpl*pol0.x() + zyPplKpl*pol0.y() + zzPplKpl*pol0.z();
165 }
166 // initial electron Kmn
167 if (!pol1.IsZero()) {
168 G4double xxPplKmn = ((-1 + e)*(e*(-2 + gamma)*gmo + gamma))/(4*e*gpo32*std::sqrt(1 + e2*gmo + gamma - 2*e*gamma));
169 G4double xyPplKmn = 0;
170 G4double xzPplKmn = (-1 + e*gmo + gmo*gamma)/(2*sqrttwo*gpo2* std::sqrt(e*(1 + e + gamma - e*gamma)));
171 G4double yxPplKmn = 0;
172 G4double yyPplKmn = (-1 - 2*gamma + e*gmo*(3 + gamma))/(4*e*gpo2);
173 G4double yzPplKmn = 0;
174 G4double zxPplKmn = (1 + 2*e2*gmo2 + gamma + gamma2 + e*(1 + (3 - 4*gamma)*gamma))/
175 (2*sqrttwo*gpo2*std::sqrt(e*(1 + e + gamma - e*gamma)));
176 G4double zyPplKmn = 0;
177 G4double zzPplKmn = -(std::sqrt((1 - e)/(e - e*gamma2 + gpo2))*
178 (2*e2*gmo2 + gamma + 2*gamma2 + e*(2 + gamma - 3*gamma2)))/(4*e*gpo);
179
180 phi2[0] += xxPplKmn*pol1.x() + xyPplKmn*pol1.y() + xzPplKmn*pol1.z();
181 phi2[1] += yxPplKmn*pol1.x() + yyPplKmn*pol1.y() + yzPplKmn*pol1.z();
182 phi2[2] += zxPplKmn*pol1.x() + zyPplKmn*pol1.y() + zzPplKmn*pol1.z();
183 }
184//
185// Final Electron Pmn
186//
187 // initial positron Kpl
188 if (!pol0.IsZero()) {
189 G4double xxPmnKpl = ((-1 + e*gmo)*(2 + gamma))/(4*gpo* std::sqrt(e*(2 + e*gmo)*gpo));
190 G4double xyPmnKpl = 0;
191 G4double xzPmnKpl = (std::sqrt((-1 + e)/(-2 + e - e*gamma))*
192 (e + gamma + e*gamma - 2*(-1 + e)*gamma2))/(2*sqrttwo*e*gpo2);
193 G4double yxPmnKpl = 0;
194 G4double yyPmnKpl = (-1 - 2*gamma + e*gmo*(3 + gamma))/(4*e*gpo2);
195 G4double yzPmnKpl = 0;
196 G4double zxPmnKpl = -((-1 + e)*(1 + 2*e*gmo)*(e*gmo - gamma))/
197 (2*sqrttwo*e*std::sqrt(-((-1 + e)*(2 + e*gmo)))*gpo2);
198 G4double zyPmnKpl = 0;
199 G4double zzPmnKpl = (-2 + 2*e2*gmo2 + gamma*(-1 + 2*gamma) +
200 e*(-2 + (5 - 3*gamma)*gamma))/(4*std::sqrt(e*(2 + e*gmo))* gpo32);
201
202 phi3[0] += xxPmnKpl*pol0.x() + xyPmnKpl*pol0.y() + xzPmnKpl*pol0.z();
203 phi3[1] += yxPmnKpl*pol0.x() + yyPmnKpl*pol0.y() + yzPmnKpl*pol0.z();
204 phi3[2] += zxPmnKpl*pol0.x() + zyPmnKpl*pol0.y() + zzPmnKpl*pol0.z();
205 }
206 // initial electron Kmn
207 if (!pol1.IsZero()) {
208 G4double xxPmnKmn = -((2 + e*gmo)*(-1 + e*gmo - gamma)*(e*gmo - gamma)*
209 (-2 + gamma))/(4*gmo*egmp32);
210 G4double xyPmnKmn = 0;
211 G4double xzPmnKmn = ((e*gmo - gamma)*
212 std::sqrt((-1 + e + gamma - e*gamma)/(2 + e*gmo))*
213 (e + gamma - e*gamma + gamma2))/
214 (2*sqrttwo*e2*gmo32*gpo2);
215 G4double yxPmnKmn = 0;
216 G4double yyPmnKmn = (gamma2*gpo + e2*gmo2*(3 + gamma) -
217 e*gmo*(1 + 2*gamma*(2 + gamma)))/(4*e2*gmo*gpo2);
218 G4double yzPmnKmn = 0;
219 G4double zxPmnKmn = -((-1 + e)*(e*gmo - gamma)*(e*gmo + 2*e2*gmo2 - gamma*gpo))/
220 (2*sqrttwo*e2*std::sqrt(-((-1 + e)*(2 + e*gmo)))* gmo*gpo2);
221 G4double zyPmnKmn = 0;
222 G4double zzPmnKmn = ((e*gmo - gamma)*std::sqrt(e/((2 + e*gmo)*gpo))*
223 (-(e*(-2 + gamma)*gmo) + 2*e2*gmo2 + (-2 + gamma)*gpo))/(4*e2*(-1 + gamma2));
224
225 phi3[0] += xxPmnKmn*pol1.x() + xyPmnKmn*pol1.y() + xzPmnKmn*pol1.z();
226 phi3[1] += yxPmnKmn*pol1.x() + yyPmnKmn*pol1.y() + yzPmnKmn*pol1.z();
227 phi3[2] += zxPmnKmn*pol1.x() + zyPmnKmn*pol1.y() + zzPmnKmn*pol1.z();
228 }
229 }
230 phi0 *= pref;
231 phi2 *= pref;
232 phi3 *= pref;
233
234}
235
236G4double G4PolarizedBhabhaCrossSection::XSection(const G4StokesVector & pol2,
237 const G4StokesVector & pol3)
238{
239 G4double xs=0.;
240 xs+=phi0;
241
242 G4bool polarized=(!pol2.IsZero())||(!pol3.IsZero());
243 if (polarized) {
244 xs+=phi2*pol2 + phi3*pol3;
245 }
246 return xs;
247}
248
249G4double G4PolarizedBhabhaCrossSection::TotalXSection(
250 G4double xmin, G4double xmax, G4double gamma,
251 const G4StokesVector & pol0,const G4StokesVector & pol1)
252{
253 G4double xs=0.;
254
255 G4double x=xmin;
256
257 if (xmax != 1.) G4cout<<" warning xmax expected to be 1 but is "<<xmax<< G4endl;
258
259 // re -> electron radius^2;
260 G4double re2 = classic_electr_radius * classic_electr_radius;
261 G4double gamma2=gamma*gamma;
262 G4double gmo2 = (gamma - 1.)*(gamma - 1.);
263 G4double gpo2 = (gamma + 1.)*(gamma + 1.);
264 G4double gpo3 = gpo2*(gamma + 1.);
265 G4double logMEM = std::log(x);
266 G4double pref = twopi*re2/(gamma - 1.0);
267 // unpolarise XS
268 G4double sigma0 = 0.;
269 sigma0 += -gmo2*(gamma - 1.)*x*x*x/3. + gmo2*gamma*x*x;
270 sigma0 += -(gamma - 1.)*(3.*gamma*(gamma + 2.) +4.)*x;
271 sigma0 += (gamma*(gamma*(gamma*(4.*gamma - 1.) - 21.) - 7.)+13.)/(3.*(gamma - 1.));
272 sigma0 /= gpo3;
273 sigma0 += logMEM*(2. - 1./gpo2);
274 sigma0 += gamma2/((gamma2 - 1.)*x);
275 // longitudinal part
276 G4double sigma2=0.;
277 sigma2 += logMEM*gamma*(gamma + 1.)*(2.*gamma + 1.);
278 sigma2 += gamma*(7.*gamma*(gamma + 1.) - 2.)/3.;
279 sigma2 += -(3.*gamma + 1.)*(gamma2 + gamma - 1.)*x;
280 sigma2 += (gamma - 1.)*gamma*(gamma + 3.)*x*x;
281 sigma2 += -gmo2*(gamma + 3.)*x*x*x/3.;
282 sigma2 /= gpo3;
283 // transverse part
284 G4double sigma3=0.;
285 sigma3 += 0.5*(gamma + 1.)*(3.*gamma + 1.)*logMEM;
286 sigma3 += (gamma*(5.*gamma - 4.) - 13.)/6.;
287 sigma3 += 0.5*(gamma2 + 3.)*x;
288 sigma3 += - 2.*(gamma - 1.)*gamma*x*x; // *AS* changed sign
289 sigma3 += 2.*gmo2*x*x*x/3.;
290 sigma3 /= gpo3;
291 // total cross section
292 xs+=pref*(sigma0 + sigma2*pol0.z()*pol1.z() + sigma3*(pol0.x()*pol1.x()+pol0.y()*pol1.y()));
293
294 return xs;
295}
296
297
298G4StokesVector G4PolarizedBhabhaCrossSection::GetPol2()
299{
300 // Note, mean polarization can not contain correlation
301 // effects.
302 return 1./phi0 * phi2;
303}
304G4StokesVector G4PolarizedBhabhaCrossSection::GetPol3()
305{
306 // Note, mean polarization can not contain correlation
307 // effects.
308 return 1./phi0 * phi3;
309}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
G4PolarizedBhabhaCrossSection::TotalXSection
G4double TotalXSection(G4double xmin, G4double xmax, G4double y, const G4StokesVector &pol0, const G4StokesVector &pol1)
Definition: G4PolarizedBhabhaCrossSection.cc:249
G4PolarizedBhabhaCrossSection::XSection
G4double XSection(const G4StokesVector &pol2, const G4StokesVector &pol3)
Definition: G4PolarizedBhabhaCrossSection.cc:236
G4PolarizedBhabhaCrossSection::Initialize
void Initialize(G4double x, G4double y, G4double phi, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0)
Definition: G4PolarizedBhabhaCrossSection.cc:62
G4StokesVector::IsZero
G4bool IsZero() const
Definition: G4StokesVector.hh:78