Geant4  10.01.p03
G4PolarizedComptonCrossSection.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 // $Id: G4PolarizedComptonCrossSection.cc 68046 2013-03-13 14:31:38Z gcosmo $
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4PolarizedComptonCrossSection
32 //
33 // Author: Andreas Schaelicke
34 //
35 // Creation date: 15.05.2005
36 //
37 // Modifications:
38 //
39 // Class Description:
40 // determine the polarization of the final state
41 // in a Compton scattering process employing the differential
42 // cross section by F.W.Lipps & H.A.Tolhoek
43 // ( Physica 20 (1954) 395 )
44 // recalculated by P.Starovoitov
45 //
47 #include "G4PhysicalConstants.hh"
48 #include "Randomize.hh"
49 
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52  : gammaPol2(false), electronPol3(false)
53 {
54  SetYmin(0.);
55 
56  // G4cout<<"G4PolarizedComptonCrossSection() init\n";
57 
58  re2 = classic_electr_radius * classic_electr_radius * sqr(4*pi/hbarc);
59  // G4double unit_conversion = hbarc_squared ;
60  // G4cout<<" (keV)^2* m^2 ="<<unit_conversion<<"\n";
61  phi0 = 0.; polXS = 0.; unpXS = 0.;
62  phi2 = G4ThreeVector(0., 0., 0.);
63  phi3 = G4ThreeVector(0., 0., 0.);
64  polxx = polyy = polzz = polxz = polzx = polyz = polzy = polxy = polyx = 0.;
65  diffXSFactor = 1.;
66  totalXSFactor = 1.;
67 }
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 {}
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75  const G4StokesVector & pol0,
76  const G4StokesVector & pol1,
77  G4int flag)
78 {
79  G4double cosT = 1. - (1./eps - 1.)/X;
80  if(cosT > 1.+1.e-8) cosT = 1.;
81  if(cosT < -1.-1.e-8) cosT = -1.;
82  G4double cosT2 = cosT*cosT;
83  G4double cosT3 = cosT2*cosT;
84  G4double sinT2 = 1. - cosT2;
85  if(sinT2 > 1. + 1.e-8) sinT2 = 1.;
86  if(sinT2 < 0.) sinT2 = 0.;
87  G4double sinT = std::sqrt(sinT2);
88  G4double cos2T = 2.*cosT2 - 1.;
89  G4double sin2T = 2.*sinT*cosT;
90  G4double eps2 = sqr(eps);
91  DefineCoefficients(pol0,pol1);
92  diffXSFactor = re2/(4.*X);
93 
94  // unpolarized Cross Section
95  unpXS = (eps2 + 1. - eps*sinT2)/(2.*eps);
96  // initial polarization dependence
97  polXS = -sinT2*pol0.x() + (1. - eps)*sinT*polzx + ((eps2 - 1.)/eps)*cosT*polzz;
98  polXS *= 0.5;
99 
100  phi0 = unpXS + polXS;
101 
102  if (flag == 2 ){
103  // polarization of outgoing photon
104  G4double PHI21 = -sinT2 + 0.5*(cos2T + 3.)*pol0.x() - ((1. - eps)/eps)*sinT*polzx;
105  PHI21 *= 0.5;
106  G4double PHI22 = cosT*pol0.y() + ((1. - eps)/(2.*eps))*sinT*polzy;
107  G4double PHI23 = ((eps2 + 1.)/eps)*cosT*pol0.z() - ((1. - eps)/eps)*(eps*cosT2 + 1.)*pol1.z();
108  PHI23 += 0.5*(1. - eps)*sin2T*pol1.x();
109  PHI23 += (eps - 1.)*(-sinT2*polxz + sinT*polyy - 0.5*sin2T*polxx);
110  PHI23 *= 0.5;
111  phi2 = G4ThreeVector(PHI21, PHI22, PHI23);
112 
113  // polarization of outgoing electron
114  G4double PHI32 = -sinT2*polxy + ((1. - eps)/eps)*sinT*polyz + 0.5*(cos2T + 3.)*pol1.y();
115  PHI32 *= 0.5;
116 
117  G4double PHI31 = 0., PHI31add = 0., PHI33 = 0., PHI33add = 0.;
118 
119  if ((1. - eps) > 1.e-12){
120  G4double helpVar = std::sqrt(eps2 - 2.*cosT*eps + 1.);
121 
122  PHI31 = (1. - eps)*(1. + cosT)*sinT*pol0.z();
123  PHI31 += (-eps*cosT3 + eps*cosT2 + (eps - 2.)*cosT + eps)*pol1.x();
124  PHI31 += -(eps*cosT2 - eps*cosT + cosT + 1.)*sinT*pol1.z();
125  PHI31 /= 2.*helpVar;
126 
127  PHI31add = -eps*sqr(1. - cosT)*(1. + cosT)*polxx;
128  PHI31add += (1. - eps)*sinT2*polyy;
129  PHI31add += -(-eps2 + cosT*(cosT*eps - eps + 1.)*eps + eps - 1.)*sinT*polxz/eps;
130  PHI31add /= 2.*helpVar;
131 
132  PHI33 = ((1. - eps)/eps)*(-eps*cosT2 + eps*(eps + 1.)*cosT - 1.)*pol0.z();
133  PHI33 += -(eps*cosT2 + (1. - eps)*eps*cosT + 1.)*sinT*pol1.x();
134  PHI33 += -(-eps2*cosT3 + eps*(eps2 - eps + 1.)*cosT2 - cosT + eps2)*pol1.z()/eps;
135  PHI33 /= -2.*helpVar;
136 
137  PHI33add = (eps*(eps - cosT - 1.)*cosT + 1.)*sinT*polxx;
138  PHI33add += -(-eps2 + cosT*eps + eps - 1.)*sinT2*polxz;
139  PHI33add += (eps - 1.)*(cosT - eps)*sinT*polyy;
140  PHI33add /= -2.*helpVar;
141  }else{
142  PHI31 = -pol1.z() - (X - 1.)*std::sqrt(1. - eps)*pol1.x()/std::sqrt(2.*X);
143  PHI31add = -(-X*X*pol1.z() - 2.*X*(2.*pol0.z() - pol1.z()) - (4.*pol0.x() + 5.)*pol1.z())*(1. - eps)/(4.*X);
144 
145  PHI33 = pol1.x() - (X - 1.)*std::sqrt(1. - eps)*pol1.z()/std::sqrt(2.*X);
146  PHI33add = -(X*X - 2.*X + 4.*pol0.x() + 5.)*(1. - eps)*pol1.x()/(4.*X);
147  }
148  phi3 = G4ThreeVector(PHI31 + PHI31add, PHI32, PHI33 + PHI33add);
149 
150  }
151  unpXS *= diffXSFactor;
152  polXS *= diffXSFactor;
153  phi0 *= diffXSFactor;
154  phi2 *= diffXSFactor;
155  phi3 *= diffXSFactor;
156 
157 }
158 
160 {
161  gammaPol2 = !(pol2==G4StokesVector::ZERO);
163 
164  G4double phi = 0.;
165  // polarization independent part
166  phi += phi0;
167 
168 
169  if (gammaPol2) {
170  // part depending on the polarization of the final photon
171  phi += phi2*pol2;
172  }
173 
174  if (electronPol3) {
175  // part depending on the polarization of the final electron
176  phi += phi3*pol3;
177  }
178 
179  // return cross section.
180  return phi;
181 }
182 
183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
184 
185 
187  const G4StokesVector & pol0,
188  const G4StokesVector & pol1)
189 {
190 
191  // G4double k0 = gammaEnergy / electron_mass_c2 ;
192  G4double k1 = 1 + 2*k0 ;
193 
194 // // pi*re^2
195 // G4double re=2.81794e-15; //m
196 // G4double barn=1.e-28; //m^2
197  G4double Z=theZ;
198 
199  G4double unit = Z*pi*classic_electr_radius
200  * classic_electr_radius ; // *1./barn;
201 
202  G4double pre = unit/(sqr(k0)*sqr(1.+2.*k0));
203 
204  G4double xs_0 = ((k0 - 2.)*k0 -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.);
205  G4double xs_pol = (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.);
206 
207  return pre*(xs_0/k0 + pol0.p3()*pol1.z()*xs_pol);
208 }
209 
210 
211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
212 
213 
215 {
216  // Note, mean polarization can not contain correlation
217  // effects.
218  return 1./phi0 * phi2;
219 }
220 
221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
222 
224 {
225  // Note, mean polarization can not contain correlation
226  // effects.
227  return 1./phi0 * phi3;
228 }
229 
230 
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
234  const G4StokesVector & pol1)
235 {
236  polxx=pol0.x()*pol1.x();
237  polyy=pol0.y()*pol1.y();
238  polzz=pol0.z()*pol1.z();
239 
240  polxz=pol0.x()*pol1.z();
241  polzx=pol0.z()*pol1.x();
242 
243  polyz=pol0.y()*pol1.z();
244  polzy=pol0.z()*pol1.y();
245 
246  polxy=pol0.x()*pol1.y();
247  polyx=pol0.y()*pol1.x();
248 }
249 
250 
251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void DefineCoefficients(const G4StokesVector &pol0, const G4StokesVector &pol1)
G4double XSection(const G4StokesVector &pol2, const G4StokesVector &pol3)
CLHEP::Hep3Vector G4ThreeVector
const G4double pi
static const G4double eps
int G4int
Definition: G4Types.hh:78
G4double p3() const
G4double TotalXSection(G4double xmin, G4double xmax, G4double y, const G4StokesVector &pol0, const G4StokesVector &pol1)
virtual void Initialize(G4double eps, G4double X, G4double phi, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0)
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
static const G4StokesVector ZERO