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