ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolarizedMollerCrossSection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PolarizedMollerCrossSection.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: G4PolarizedMollerCrossSection
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 //
40 // Class Description:
41 // * calculates the differential cross section
42 // incomming electron K1(along positive z direction) scatters at an electron K2 at rest
43 // * phi denotes the angle between the scattering plane (defined by the
44 // outgoing electron) and X-axis
45 // * all stokes vectors refer to spins in the Global System (X,Y,Z)
46 //
47 
48 #include "G4PolarizedMollerCrossSection.hh"
49 #include "G4PhysicalConstants.hh"
50 
51 G4PolarizedMollerCrossSection::G4PolarizedMollerCrossSection() :
52  phi0(0.)
53 {
54  SetXmax(.5);
55 }
56 G4PolarizedMollerCrossSection::~G4PolarizedMollerCrossSection() {}
57 void G4PolarizedMollerCrossSection::Initialize(
58  G4double e,
59  G4double gamma,
60  G4double /*phi*/,
61  const G4StokesVector & pol0,
62  const G4StokesVector & pol1,
63  G4int flag)
64 {
65  G4double re2 = classic_electr_radius * classic_electr_radius;
66  G4double gamma2=gamma*gamma;
67  G4double gmo = (gamma - 1.);
68  G4double gmo2 = (gamma - 1.)*(gamma - 1.);
69  G4double gpo = (gamma + 1.);
70  G4double pref = gamma2*re2/(gmo2*(gamma + 1.0));
71  G4double sqrttwo=std::sqrt(2.);
72  G4double f = (-1. + e);
73  G4double e2 = e*e;
74  G4double f2 = f*f;
75  // G4double w = e*(1. - e);
76 
77  G4bool polarized=(!pol0.IsZero())||(!pol1.IsZero());
78 
79  if (flag==0) polarized=false;
80  // Unpolarised part of XS
81  phi0 = 0.;
82  phi0+= gmo2/gamma2;
83  phi0+= ((1. - 2.*gamma)/gamma2)*(1./e + 1./(1.-e));
84  phi0+= 1./(e*e) + 1./((1. - e)*(1. - e));
85  phi0*=0.25;
86  // Initial state polarisarion dependence
87  if (polarized) {
88  G4double usephi=1.;
89  if (flag<=1) usephi=0.;
90  // G4cout<<"Polarized differential moller cross section"<<G4endl;
91  // G4cout<<"Initial state polarisation contributions"<<G4endl;
92  // G4cout<<"Diagonal Matrix Elements"<<G4endl;
93  G4double xx = (gamma - f*e*gmo*(3. + gamma))/(4.*f*e*gamma2);
94  G4double yy = (-1. + f*e*gmo2 + 2.*gamma)/(4.*f*e*gamma2);
95  G4double zz = (-(e*gmo*(3. + gamma)) + e2*gmo*(3. + gamma) +
96  gamma*(-1. + 2.*gamma))/(4.*f*e*gamma2);
97 
98  phi0 += xx*pol0.x()*pol1.x() + yy*pol0.y()*pol1.y() + zz*pol0.z()*pol1.z();
99 
100  if (usephi==1.) {
101  // G4cout<<"Non-diagonal Matrix Elements"<<G4endl;
102  G4double xy = 0.;
103  G4double xz = -((-1. + 2.*e)*gmo)/(2.*sqrttwo*gamma2*
104  std::sqrt(-((f*e)/gpo)));
105  G4double yx = 0.;
106  G4double yz = 0.;
107  G4double zx = -((-1. + 2.*e)*gmo)/(2.*sqrttwo*gamma2*
108  std::sqrt(-((f*e)/gpo)));
109  G4double zy = 0.;
110  phi0+=yx*pol0.y()*pol1.x() + xy*pol0.x()*pol1.y();
111  phi0+=zx*pol0.z()*pol1.x() + xz*pol0.x()*pol1.z();
112  phi0+=zy*pol0.z()*pol1.y() + yz*pol0.y()*pol1.z();
113  }
114  }
115  // Final state polarisarion dependence
116  phi2=G4ThreeVector();
117  phi3=G4ThreeVector();
118 
119  if (flag>=1) {
120  //
121  // Final Electron P1
122  //
123 
124  // initial electron K1
125  if (!pol0.IsZero()) {
126  G4double xxP1K1 = (std::sqrt(gpo/(1. + e2*gmo + gamma - 2.*e*gamma))*
127  (gamma - e*gpo))/(4.*e2*gamma);
128  G4double xyP1K1 = 0.;
129  G4double xzP1K1 = (-1. + 2.*e*gamma)/(2.*sqrttwo*f*gamma*
130  std::sqrt(e*e2*(1. + e + gamma - e*gamma)));
131  G4double yxP1K1 = 0.;
132  G4double yyP1K1 = (-gamma2 + e*(-1. + gamma*(2. + gamma)))/(4.*f*e2*gamma2);
133  G4double yzP1K1 = 0.;
134  G4double zxP1K1 = (1. + 2.*e2*gmo - 2.*e*gamma)/(2.*sqrttwo*f*e*gamma*
135  std::sqrt(e*(1. + e + gamma - e*gamma)));
136  G4double zyP1K1 = 0.;
137  G4double zzP1K1 = (-gamma + e*(1. - 2.*e*gmo + gamma))/(4.*f*e2*gamma*
138  std::sqrt(1. - (2.*e)/(f*gpo)));
139  phi2[0] += xxP1K1*pol0.x() + xyP1K1*pol0.y() + xzP1K1*pol0.z();
140  phi2[1] += yxP1K1*pol0.x() + yyP1K1*pol0.y() + yzP1K1*pol0.z();
141  phi2[2] += zxP1K1*pol0.x() + zyP1K1*pol0.y() + zzP1K1*pol0.z();
142  }
143  // initial electron K2
144  if (!pol1.IsZero()) {
145  G4double xxP1K2 = ((1. + e*(-3. + gamma))*std::sqrt(gpo/(1. + e2*gmo + gamma -
146  2.*e*gamma)))/(4.*f*e*gamma);
147  G4double xyP1K2 = 0.;
148  G4double xzP1K2 = (-2. + 2.*e + gamma)/(2.*sqrttwo*f2*gamma*
149  std::sqrt(e*(1. + e + gamma - e*gamma)));
150  G4double yxP1K2 = 0.;
151  G4double yyP1K2 = (1. - 2.*gamma + e*(-1. + gamma*(2. + gamma)))/(4.*f2*e*gamma2);
152  G4double yzP1K2 = 0.;
153  G4double zxP1K2 = (2.*e*(1. + e*gmo - 2.*gamma) + gamma)/(2.*sqrttwo*f2*gamma*
154  std::sqrt(e*(1. + e + gamma - e*gamma)));
155  G4double zyP1K2 = 0.;
156  G4double zzP1K2 = (1. - 2.*gamma + e*(-1. - 2.*e*gmo + 3.*gamma))/
157  (4.*f2*e*gamma*std::sqrt(1. - (2.*e)/(f*gpo)));
158  phi2[0] += xxP1K2*pol1.x() + xyP1K2*pol1.y() + xzP1K2*pol1.z();
159  phi2[1] += yxP1K2*pol1.x() + yyP1K2*pol1.y() + yzP1K2*pol1.z();
160  phi2[2] += zxP1K2*pol1.x() + zyP1K2*pol1.y() + zzP1K2*pol1.z();
161  }
162  //
163  // Final Electron P2
164  //
165 
166  // initial electron K1
167  if (!pol0.IsZero()) {
168 
169 
170  G4double xxP2K1 = (-1. + e + e*gamma)/(4.*f2*gamma*
171  std::sqrt((e*(2. + e*gmo))/gpo));
172  G4double xyP2K1 = 0.;
173  G4double xzP2K1 = -((1. + 2.*f*gamma)*std::sqrt(f/(-2. + e - e*gamma)))/
174  (2.*sqrttwo*f2*e*gamma);
175  G4double yxP2K1 = 0.;
176  G4double yyP2K1 = (1. - 2.*gamma + e*(-1. + gamma*(2. + gamma)))/(4.*f2*e*gamma2);
177  G4double yzP2K1 = 0.;
178  G4double zxP2K1 = (1. + 2.*e*(-2. + e + gamma - e*gamma))/(2.*sqrttwo*f*e*
179  std::sqrt(-(f*(2. + e*gmo)))*gamma);
180  G4double zyP2K1 = 0.;
181  G4double zzP2K1 = (std::sqrt((e*gpo)/(2. + e*gmo))*
182  (-3. + e*(5. + 2.*e*gmo - 3.*gamma) + 2.*gamma))/(4.*f2*e*gamma);
183 
184  phi3[0] += xxP2K1*pol0.x() + xyP2K1*pol0.y() + xzP2K1*pol0.z();
185  phi3[1] += yxP2K1*pol0.x() + yyP2K1*pol0.y() + yzP2K1*pol0.z();
186  phi3[2] += zxP2K1*pol0.x() + zyP2K1*pol0.y() + zzP2K1*pol0.z();
187  }
188  // initial electron K2
189  if (!pol1.IsZero()) {
190 
191  G4double xxP2K2 = (-2. - e*(-3. + gamma) + gamma)/
192  (4.*f*e*gamma* std::sqrt((e*(2. + e*gmo))/gpo));
193  G4double xyP2K2 = 0.;
194  G4double xzP2K2 = ((-2.*e + gamma)*std::sqrt(f/(-2. + e - e*gamma)))/
195  (2.*sqrttwo*f*e2*gamma);
196  G4double yxP2K2 = 0.;
197  G4double yyP2K2 = (-gamma2 + e*(-1. + gamma*(2. + gamma)))/(4.*f*e2*gamma2);
198  G4double yzP2K2 = 0.;
199  G4double zxP2K2 = (gamma + 2.*e*(-1. + e - e*gamma))/
200  (2.*sqrttwo*e2* std::sqrt(-(f*(2. + e*gmo)))*gamma);
201  G4double zyP2K2 = 0.;
202  G4double zzP2K2 = (std::sqrt((e*gpo)/(2. + e*gmo))*
203  (-2. + e*(3. + 2.*e*gmo - gamma) + gamma))/(4.*f*e2*gamma);
204  phi3[0] += xxP2K2*pol1.x() + xyP2K2*pol1.y() + xzP2K2*pol1.z();
205  phi3[1] += yxP2K2*pol1.x() + yyP2K2*pol1.y() + yzP2K2*pol1.z();
206  phi3[2] += zxP2K2*pol1.x() + zyP2K2*pol1.y() + zzP2K2*pol1.z();
207  }
208  }
209  phi0 *= pref;
210  phi2 *= pref;
211  phi3 *= pref;
212 }
213 
214 G4double G4PolarizedMollerCrossSection::XSection(const G4StokesVector & pol2,
215  const G4StokesVector & pol3)
216 {
217  G4double xs=0.;
218  xs+=phi0;
219 
220  G4bool polarized=(!pol2.IsZero())||(!pol3.IsZero());
221  if (polarized) {
222  xs+=phi2*pol2 + phi3*pol3;
223  }
224  return xs;
225 }
226 
227 G4double G4PolarizedMollerCrossSection::TotalXSection(
228  G4double xmin, G4double xmax, G4double gamma,
229  const G4StokesVector & pol0,const G4StokesVector & pol1)
230 {
231  G4double xs=0.;
232 
233  G4double x=xmin;
234 
235  if (xmax != 1./2.) G4cout<<" warning xmax expected to be 1/2 but is "<<xmax<< G4endl;
236 
237  // re -> electron radius^2;
238  G4double re2 = classic_electr_radius * classic_electr_radius;
239  G4double gamma2=gamma*gamma;
240  G4double gmo2 = (gamma - 1.)*(gamma - 1.);
241  G4double logMEM = std::log(1./x - 1.);
242  G4double pref = twopi*gamma2*re2/(gmo2*(gamma + 1.0));
243  // unpolarise XS
244  G4double sigma0 = 0.;
245  sigma0 += (gmo2/gamma2)*(0.5 - x);
246  sigma0 += ((1. - 2.*gamma)/gamma2)*logMEM;
247  sigma0 += 1./x - 1./(1. - x);
248  // longitudinal part
249  G4double sigma2=0.;
250  sigma2 += ((gamma2 + 2.*gamma - 3.)/gamma2)*(0.5 - x);
251  sigma2 += (1./gamma - 2.)*logMEM;
252  // transverse part
253  G4double sigma3=0.;
254  sigma3 += (2.*(1. - gamma)/gamma2)*(0.5 - x);
255  sigma3 += (1. - 3.*gamma)/(2.*gamma2)*logMEM;
256  // total cross section
257  xs+=pref*(sigma0 + sigma2*pol0.z()*pol1.z() + sigma3*(pol0.x()*pol1.x()+pol0.y()*pol1.y()));
258 
259  return xs;
260 }
261 
262 
263 G4StokesVector G4PolarizedMollerCrossSection::GetPol2()
264 {
265  // Note, mean polarization can not contain correlation
266  // effects.
267  return 1./phi0 * phi2;
268 }
269 G4StokesVector G4PolarizedMollerCrossSection::GetPol3()
270 {
271  // Note, mean polarization can not contain correlation
272  // effects.
273  return 1./phi0 * phi3;
274 }