source: trunk/source/processes/electromagnetic/polarisation/src/G4PolarizedMollerCrossSection.cc@ 1317

Last change on this file since 1317 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 10.1 KB
Line 
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: G4PolarizedMollerCrossSection.cc,v 1.5 2007/11/01 17:32:34 schaelic Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28// -------------------------------------------------------------------
29//
30// GEANT4 Class file
31//
32//
33// File name: G4PolarizedMollerCrossSection
34//
35// Author: Andreas Schaelicke
36//
37// Creation date: 12.01.2006
38//
39// Modifications:
40// 16-01-06 included cross section as calculated by P.Starovoitov
41//
42// Class Description:
43// * calculates the differential cross section
44// incomming electron K1(along positive z direction) scatters at an electron K2 at rest
45// * phi denotes the angle between the scattering plane (defined by the
46// outgoing electron) and X-axis
47// * all stokes vectors refer to spins in the Global System (X,Y,Z)
48//
49
50#include "G4PolarizedMollerCrossSection.hh"
51
52G4PolarizedMollerCrossSection::G4PolarizedMollerCrossSection()
53{
54 SetXmax(.5);
55}
56G4PolarizedMollerCrossSection::~G4PolarizedMollerCrossSection() {}
57void 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
214G4double 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
227G4double 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
263G4StokesVector G4PolarizedMollerCrossSection::GetPol2()
264{
265 // Note, mean polarization can not contain correlation
266 // effects.
267 return 1./phi0 * phi2;
268}
269G4StokesVector G4PolarizedMollerCrossSection::GetPol3()
270{
271 // Note, mean polarization can not contain correlation
272 // effects.
273 return 1./phi0 * phi3;
274}
Note: See TracBrowser for help on using the repository browser.