Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolarizedAnnihilationCrossSection.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 // -------------------------------------------------------------------
27 // $Id: G4PolarizedAnnihilationCrossSection.cc 68046 2013-03-13 14:31:38Z gcosmo $
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4PolarizedAnnihilationCrossSection
34 //
35 // Author: Andreas Schaelicke
36 //
37 // Creation date: 22.03.2006
38 //
39 // Modifications:
40 // 24-03-06 included cross section as given in W.McMaster, Nuovo Cimento 7, 1960, 395
41 // 27-07-06 new calculation by P.Starovoitov
42 // 15.10.07 introduced a more general framework for cross sections (AS)
43 // 16.10.07 some minor corrections in formula longPart
44 //
45 // Class Description:
46 // * calculates the differential cross section in e+e- -> gamma gamma
47 //
48 
50 #include "G4PhysicalConstants.hh"
51 
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 {
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 
70 {
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 
75 void 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 
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 
230  const G4StokesVector & pol3)
231 {
232  G4double xs=phi0+pol2*phi2+pol3*phi3;
233  return xs;
234 }
235 //
236 // calculate total cross section
237 //
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......
273 {
274  // Note, mean polarization can not contain correlation
275  // effects.
276  return 1./phi0 * phi2;
277 }
278 
279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
280 
282 {
283  // Note, mean polarization can not contain correlation
284  // effects.
285  return 1./phi0 * phi3;
286 }
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 void 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 
308 {
309  return 0.5*(1.-std::sqrt((y-1.)/(y+1.)));
310 }
312 {
313  return 0.5*(1.+std::sqrt((y-1.)/(y+1.)));
314 }
315 
316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
318 {
319  return dice;
320 }
321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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......
CLHEP::Hep3Vector G4ThreeVector
double x() const
static const G4double eps
int G4int
Definition: G4Types.hh:78
void setY(double)
double z() const
void setZ(double)
void setX(double)
virtual G4double XSection(const G4StokesVector &pol2, const G4StokesVector &pol3) override
virtual G4double TotalXSection(G4double xmin, G4double xmax, G4double y, const G4StokesVector &pol0, const G4StokesVector &pol1) override
double y() const
static constexpr double pi
Definition: G4SIunits.hh:75
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
virtual void Initialize(G4double eps, G4double gamma, G4double phi, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0) override