source: trunk/source/processes/electromagnetic/polarisation/src/G4PolarizedAnnihilationCrossSection.cc@ 1199

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 11.4 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// -------------------------------------------------------------------
27// $Id: G4PolarizedAnnihilationCrossSection.cc,v 1.6 2007/11/01 17:32:34 schaelic Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4PolarizedAnnihilationCrossSection
35//
36// Author: Andreas Schaelicke
37//
38// Creation date: 22.03.2006
39//
40// Modifications:
41// 24-03-06 included cross section as given in W.McMaster, Nuovo Cimento 7, 1960, 395
42// 27-07-06 new calculation by P.Starovoitov
43// 15.10.07 introduced a more general framework for cross sections (AS)
44// 16.10.07 some minor corrections in formula longPart
45//
46// Class Description:
47// * calculates the differential cross section in e+e- -> gamma gamma
48//
49
50#include "G4PolarizedAnnihilationCrossSection.hh"
51
52G4PolarizedAnnihilationCrossSection::G4PolarizedAnnihilationCrossSection() :
53 polxx(0.), polyy(0.), polzz(0.), polxz(0.), polzx(0.), polxy(0.),
54 polyx(0.), polyz(0.), polzy(0.),
55 re2(1.), diffXSFactor(1.), totalXSFactor(1.),
56 phi0(0.)
57{
58 re2 = classic_electr_radius * classic_electr_radius;
59 phi2 = G4ThreeVector(0., 0., 0.);
60 phi3 = G4ThreeVector(0., 0., 0.);
61 dice = 0.;
62 polXS= 0.;
63 unpXS = 0.;
64 ISPxx=ISPyy=ISPzz=ISPnd=0.;
65}
66
67//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68
69G4PolarizedAnnihilationCrossSection::~G4PolarizedAnnihilationCrossSection()
70{
71}
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74
75void G4PolarizedAnnihilationCrossSection::TotalXS()
76{
77 // total cross section is sum of
78 // + unpol xsec "sigma0"
79 // + longitudinal polarised cross section "sigma_zz" times pol_3(e-)*pol_3(e+)
80 // + transverse contribution "(sigma_xx+sigma_yy)/2" times pol_T(e-)*pol_T(e+)
81 // (Note: if both beams are transversely polarised, i.e. pol_T(e-)!=0 and
82 // pol_T(e+)!=0, and sigma_xx!=sigma_yy, then the diff. cross section will
83 // exhibit a azimuthal asymmetry even if pol_T(e-)*pol_T(e+)==0)
84
85
86}
87
88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89
90void G4PolarizedAnnihilationCrossSection::Initialize(
91 G4double eps,
92 G4double gam,
93 G4double , // phi
94 const G4StokesVector & pol0, // positron polarization
95 const G4StokesVector & pol1, // electron polarization
96 G4int flag)
97{
98
99 diffXSFactor=re2/(gam - 1.);
100 DefineCoefficients(pol0,pol1);
101 //
102 // prepare dicing
103 //
104 dice = 0.;
105 G4double symmXS = 0.125*((-1./sqr(gam + 1.))/sqr(eps) +
106 ((sqr(gam) + 4.*gam - 1.)/sqr(gam + 1.))/eps - 1.);
107 //
108 //
109 //
110 G4ThreeVector epsVector(1./sqr(eps), 1./eps, 1.);
111 G4ThreeVector oneEpsVector(1./sqr(1. - eps), 1./(1.-eps), 1.);
112 G4ThreeVector sumEpsVector(epsVector + oneEpsVector);
113 G4ThreeVector difEpsVector(epsVector - oneEpsVector);
114 G4ThreeVector calcVector(0., 0., 0.);
115 //
116 // temporary variables
117 //
118 G4double helpVar2 = 0., helpVar1 = 0.;
119 //
120 // unpolarised contribution
121 //
122 helpVar1 = (gam*gam + 4.*gam + 1.)/sqr(gam + 1.);
123 helpVar2 = -1./sqr(gam + 1.);
124 calcVector = G4ThreeVector(helpVar2, helpVar1, -1.);
125 unpXS = 0.125 * calcVector * sumEpsVector;
126
127 // initial particles polarised contribution
128 helpVar2 = 1./sqr(gam + 1.);
129 helpVar1 = -(gam*gam + 4.*gam + 1.)/sqr(gam + 1.);
130 calcVector = G4ThreeVector(helpVar2, helpVar1, 0.5*(gam + 3.));
131 ISPxx = 0.25*(calcVector * sumEpsVector)/(gam - 1.);
132
133 helpVar1 = 1./sqr(gam + 1.);
134 calcVector = G4ThreeVector(-helpVar1, 2.*gam*helpVar1, -1.);
135 ISPyy = 0.125 * calcVector * sumEpsVector;
136
137 helpVar1 = 1./(gam - 1.);
138 helpVar2 = 1./sqr(gam + 1.);
139 calcVector = G4ThreeVector(-(gam*gam + 1.)*helpVar2,(gam*gam*(gam + 1.) + 7.*gam + 3.)*helpVar2, -(gam + 3.));
140 ISPzz = 0.125*helpVar1*(calcVector * sumEpsVector);
141
142 helpVar1 = std::sqrt(std::fabs(eps*(1. - eps)*2.*(gam + 1.) - 1.));
143 calcVector = G4ThreeVector(-1./(gam*gam - 1.), 2./(gam - 1.), 0.);
144 ISPnd = 0.125*(calcVector * difEpsVector) * helpVar1;
145
146 polXS = 0.;
147 polXS += ISPxx*polxx;
148 polXS += ISPyy*polyy;
149 polXS += ISPzz*polzz;
150 polXS += ISPnd*(polzx + polxz);
151 phi0 = unpXS + polXS;
152 dice = symmXS;
153 // if(polzz != 0.) dice *= (1. + std::fabs(polzz*ISPzz/unpXS));
154 if(polzz != 0.) {
155 dice *= (1. + (polzz*ISPzz/unpXS));
156 if (dice<0.) dice=0.;
157 }
158 // prepare final state coefficients
159 if (flag==2) {
160 //
161 // circular polarisation
162 //
163 G4double circ1 = 0., circ2 = 0., circ3 = 0.;
164 helpVar1 = 8.*sqr(1. - eps)*sqr(eps)*(gam - 1.)*sqr(gam + 1.)/std::sqrt(gam*gam - 1.);
165 helpVar2 = sqr(gam + 1.)*sqr(eps)*(-2.*eps + 3.) - (gam*gam + 3.*gam + 2.)*eps;
166 circ1 = helpVar2 + gam;
167 circ1 /= helpVar1;
168 circ2 = helpVar2 + 1.;
169 circ2 /= helpVar1;
170 helpVar1 = std::sqrt(std::fabs(eps*(1. - eps)*2.*(gam + 1.) - 1.));
171 helpVar1 /= std::sqrt(gam*gam - 1.);
172 calcVector = G4ThreeVector(1., -2.*gam, 0.);
173 circ3 = 0.125*(calcVector * sumEpsVector)/(gam + 1.);
174 circ3 *= helpVar1;
175
176 phi2.setZ( circ2*pol1.z() + circ1*pol0.z() + circ3*(pol1.x() + pol0.x()));
177 phi3.setZ(-circ1*pol1.z() - circ2*pol0.z() - circ3*(pol1.x() + pol0.x()));
178 //
179 // common to both linear polarisation
180 //
181 calcVector = G4ThreeVector(-1., 2.*gam, 0.);
182 G4double linearZero = 0.125*(calcVector * sumEpsVector)/sqr(gam + 1.);
183 //
184 // Linear Polarisation #1
185 //
186 helpVar1 = std::sqrt(std::fabs(2.*(gam + 1.)*(1. - eps)*eps - 1.))/((gam + 1.)*eps*(1. - eps));
187 helpVar2 = helpVar1*helpVar1;
188 //
189 // photon 1
190 //
191 G4double diagContrib = 0.125*helpVar2*(polxx + polyy - polzz);
192 G4double nonDiagContrib = 0.125*helpVar1*(-polxz/(1. - eps) + polzx/eps);
193
194 phi2.setX(linearZero + diagContrib + nonDiagContrib);
195 //
196 // photon 2
197 //
198 nonDiagContrib = 0.125*helpVar1*(polxz/eps - polzx/(1. - eps));
199
200
201 phi3.setX(linearZero + diagContrib + nonDiagContrib);
202 //
203 // Linear Polarisation #2
204 //
205 helpVar1 = std::sqrt(gam*gam - 1.)*(2.*(gam + 1.)*eps*(1. - eps) - 1.);
206 helpVar1 /= 8.*sqr(1. - eps)*sqr(eps)*sqr(gam + 1.)*(gam - 1.);
207 helpVar2 = std::sqrt((gam*gam - 1.)*std::fabs(2.*(gam + 1.)*eps*(1. - eps) - 1.));
208 helpVar2 /= 8.*sqr(1. - eps)*sqr(eps)*sqr(gam + 1.)*(gam - 1.);
209
210 G4double contrib21 = (-polxy + polyx)*helpVar1;
211 G4double contrib32 = -(eps*(gam + 1.) - 1.)*polyz + (eps*(gam + 1.) - gam)*polzy;
212
213 contrib32 *=helpVar2;
214 phi2.setY(contrib21 + contrib32);
215
216 contrib32 = -(eps*(gam + 1.) - gam)*polyz + (eps*(gam + 1.) - 1.)*polzy;
217 contrib32 *=helpVar2;
218 phi3.setY(contrib21 + contrib32);
219
220 }
221 phi0 *= diffXSFactor;
222 phi2 *= diffXSFactor;
223 phi3 *= diffXSFactor;
224}
225
226
227//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
228
229G4double G4PolarizedAnnihilationCrossSection::XSection(const G4StokesVector & pol2,
230 const G4StokesVector & pol3)
231{
232 G4double xs=phi0+pol2*phi2+pol3*phi3;
233 return xs;
234}
235//
236// calculate total cross section
237//
238G4double G4PolarizedAnnihilationCrossSection::TotalXSection(
239 G4double ,G4double ,G4double gam,
240 const G4StokesVector & pol0,const G4StokesVector & pol1)
241{
242 totalXSFactor =pi*re2/(gam + 1.); // atomic number ignored
243 DefineCoefficients(pol0,pol1);
244
245 G4double xs = 0.;
246
247
248 G4double gam2 = gam*gam;
249 G4double sqrtgam1 = std::sqrt(gam2 - 1.);
250 G4double logMEM = std::log(gam+sqrtgam1);
251 G4double unpME = (gam*(gam + 4.) + 1.)*logMEM;
252 unpME += -(gam + 3.)*sqrtgam1;
253 unpME /= 4.*(gam2 - 1.);
254// G4double longPart = - 2.*(gam*(gam + 4.) + 1.)*logMEM;
255// longPart += (gam*(gam + 4.) + 7.)*sqrtgam1;
256// longPart /= 4.*sqr(gam - 1.)*(gam + 1.);
257 G4double longPart = (3+gam*(gam*(gam + 1.) + 7.))*logMEM;
258 longPart += - (5.+ gam*(3*gam + 4.))*sqrtgam1;
259 longPart /= 4.*sqr(gam - 1.)*(gam + 1.);
260 G4double tranPart = -(5*gam + 1.)*logMEM;
261 tranPart += (gam + 5.)*sqrtgam1;
262 tranPart /= 4.*sqr(gam - 1.)*(gam + 1.);
263
264 xs += unpME;
265 xs += polzz*longPart;
266 xs += (polxx + polyy)*tranPart;
267
268 return xs*totalXSFactor;
269}
270
271//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
272G4StokesVector G4PolarizedAnnihilationCrossSection::GetPol2()
273{
274 // Note, mean polarization can not contain correlation
275 // effects.
276 return 1./phi0 * phi2;
277}
278
279//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
280
281G4StokesVector G4PolarizedAnnihilationCrossSection::GetPol3()
282{
283 // Note, mean polarization can not contain correlation
284 // effects.
285 return 1./phi0 * phi3;
286}
287//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288void G4PolarizedAnnihilationCrossSection::DefineCoefficients(const G4StokesVector & pol0,
289 const G4StokesVector & pol1)
290{
291 polxx=pol0.x()*pol1.x();
292 polyy=pol0.y()*pol1.y();
293 polzz=pol0.z()*pol1.z();
294
295 polxz=pol0.x()*pol1.z();
296 polzx=pol0.z()*pol1.x();
297
298 polyz=pol0.y()*pol1.z();
299 polzy=pol0.z()*pol1.y();
300
301 polxy=pol0.x()*pol1.y();
302 polyx=pol0.y()*pol1.x();
303}
304
305//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
306
307G4double G4PolarizedAnnihilationCrossSection::GetXmin(G4double y)
308{
309 return 0.5*(1.-std::sqrt((y-1.)/(y+1.)));
310}
311G4double G4PolarizedAnnihilationCrossSection::GetXmax(G4double y)
312{
313 return 0.5*(1.+std::sqrt((y-1.)/(y+1.)));
314}
315
316//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
317G4double G4PolarizedAnnihilationCrossSection::DiceEpsilon()
318{
319 return dice;
320}
321//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
322G4double G4PolarizedAnnihilationCrossSection::getVar(G4int choice)
323{
324 if (choice == -1) return polXS/unpXS;
325 if (choice == 0) return unpXS;
326 if (choice == 1) return ISPxx;
327 if (choice == 2) return ISPyy;
328 if (choice == 3) return ISPzz;
329 if (choice == 4) return ISPnd;
330 return 0;
331}
332//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.