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