| Geant4 Cross Reference |
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 // Geant4 Class file
27 //
28 // File name: G4PolarizedComptonXS
29 //
30 // Author: Andreas Schaelicke
31 //
32 // Class Description:
33 // determine the polarization of the final state in a Compton scattering
34 // process employing the differential cross section by F.W.Lipps & H.A.Tolhoek
35 // ( Physica 20 (1954) 395 )
36 // recalculated by P.Starovoitov
37
38 #include "G4PolarizedComptonXS.hh"
39
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41 G4PolarizedComptonXS::G4PolarizedComptonXS()
42 {
43 SetYmin(0.);
44
45 fPhi0 = 0.;
46 fPolXS = 0.;
47 fUnpXS = 0.;
48 fPhi2 = G4ThreeVector(0., 0., 0.);
49 fPhi3 = G4ThreeVector(0., 0., 0.);
50 polxx = polyy = polzz = polxz = polzx = polyz = polzy = polxy = polyx = 0.;
51 }
52
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 G4PolarizedComptonXS::~G4PolarizedComptonXS() {}
55
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 void G4PolarizedComptonXS::Initialize(G4double eps, G4double X,
58 G4double, // phi
59 const G4StokesVector& pol0,
60 const G4StokesVector& pol1, G4int flag)
61 {
62 G4double cosT = 1. - (1. / eps - 1.) / X;
63 if(cosT > 1. + 1.e-8)
64 cosT = 1.;
65 else if(cosT < -1. - 1.e-8)
66 cosT = -1.;
67 G4double cosT2 = cosT * cosT;
68 G4double cosT3 = cosT2 * cosT;
69 G4double sinT2 = 1. - cosT2;
70 if(sinT2 > 1. + 1.e-8)
71 sinT2 = 1.;
72 else if(sinT2 < 0.)
73 sinT2 = 0.;
74 G4double sinT = std::sqrt(sinT2);
75 G4double cos2T = 2. * cosT2 - 1.;
76 G4double sin2T = 2. * sinT * cosT;
77 G4double eps2 = sqr(eps);
78 DefineCoefficients(pol0, pol1);
79 G4double diffXSFactor = re2 / (4. * X);
80
81 // unpolarized Cross Section
82 fUnpXS = (eps2 + 1. - eps * sinT2) / (2. * eps);
83 // initial polarization dependence
84 fPolXS = -sinT2 * pol0.x() + (1. - eps) * sinT * polzx +
85 ((eps2 - 1.) / eps) * cosT * polzz;
86 fPolXS *= 0.5;
87
88 fPhi0 = fUnpXS + fPolXS;
89
90 if(flag == 2)
91 {
92 // polarization of outgoing photon
93 G4double phi21 = -sinT2 + 0.5 * (cos2T + 3.) * pol0.x() -
94 ((1. - eps) / eps) * sinT * polzx;
95 phi21 *= 0.5;
96 G4double phi22 = cosT * pol0.y() + ((1. - eps) / (2. * eps)) * sinT * polzy;
97 G4double phi23 = ((eps2 + 1.) / eps) * cosT * pol0.z() -
98 ((1. - eps) / eps) * (eps * cosT2 + 1.) * pol1.z();
99 phi23 += 0.5 * (1. - eps) * sin2T * pol1.x();
100 phi23 += (eps - 1.) * (-sinT2 * polxz + sinT * polyy - 0.5 * sin2T * polxx);
101 phi23 *= 0.5;
102
103 fPhi2 = G4ThreeVector(phi21, phi22, phi23);
104
105 // polarization of outgoing electron
106 G4double phi32 = -sinT2 * polxy + ((1. - eps) / eps) * sinT * polyz +
107 0.5 * (cos2T + 3.) * pol1.y();
108 phi32 *= 0.5;
109
110 G4double phi31 = 0.;
111 G4double phi31add = 0.;
112 G4double phi33 = 0.;
113 G4double phi33add = 0.;
114
115 if((1. - eps) > 1.e-12)
116 {
117 G4double helpVar = std::sqrt(eps2 - 2. * cosT * eps + 1.);
118
119 phi31 = (1. - eps) * (1. + cosT) * sinT * pol0.z();
120 phi31 +=
121 (-eps * cosT3 + eps * cosT2 + (eps - 2.) * cosT + eps) * pol1.x();
122 phi31 += -(eps * cosT2 - eps * cosT + cosT + 1.) * sinT * pol1.z();
123 phi31 /= 2. * helpVar;
124
125 phi31add = -eps * sqr(1. - cosT) * (1. + cosT) * polxx;
126 phi31add += (1. - eps) * sinT2 * polyy;
127 phi31add += -(-eps2 + cosT * (cosT * eps - eps + 1.) * eps + eps - 1.) *
128 sinT * polxz / eps;
129 phi31add /= 2. * helpVar;
130
131 phi33 = ((1. - eps) / eps) *
132 (-eps * cosT2 + eps * (eps + 1.) * cosT - 1.) * pol0.z();
133 phi33 += -(eps * cosT2 + (1. - eps) * eps * cosT + 1.) * sinT * pol1.x();
134 phi33 +=
135 -(-eps2 * cosT3 + eps * (eps2 - eps + 1.) * cosT2 - cosT + eps2) *
136 pol1.z() / eps;
137 phi33 /= -2. * helpVar;
138
139 phi33add = (eps * (eps - cosT - 1.) * cosT + 1.) * sinT * polxx;
140 phi33add += -(-eps2 + cosT * eps + eps - 1.) * sinT2 * polxz;
141 phi33add += (eps - 1.) * (cosT - eps) * sinT * polyy;
142 phi33add /= -2. * helpVar;
143 }
144 else
145 {
146 phi31 = -pol1.z() -
147 (X - 1.) * std::sqrt(1. - eps) * pol1.x() / std::sqrt(2. * X);
148 phi31add = -(-X * X * pol1.z() - 2. * X * (2. * pol0.z() - pol1.z()) -
149 (4. * pol0.x() + 5.) * pol1.z()) *
150 (1. - eps) / (4. * X);
151
152 phi33 = pol1.x() -
153 (X - 1.) * std::sqrt(1. - eps) * pol1.z() / std::sqrt(2. * X);
154 phi33add = -(X * X - 2. * X + 4. * pol0.x() + 5.) * (1. - eps) *
155 pol1.x() / (4. * X);
156 }
157 fPhi3 = G4ThreeVector(phi31 + phi31add, phi32, phi33 + phi33add);
158 }
159 fUnpXS *= diffXSFactor;
160 fPolXS *= diffXSFactor;
161 fPhi0 *= diffXSFactor;
162 fPhi2 *= diffXSFactor;
163 fPhi3 *= diffXSFactor;
164 }
165
166 G4double G4PolarizedComptonXS::XSection(const G4StokesVector& pol2,
167 const G4StokesVector& pol3)
168 {
169 G4bool gammaPol2 = !(pol2 == G4StokesVector::ZERO);
170 G4bool electronPol3 = !(pol3 == G4StokesVector::ZERO);
171
172 G4double phi = 0.;
173 // polarization independent part
174 phi += fPhi0;
175
176 if(gammaPol2)
177 {
178 // part depending on the polarization of the final photon
179 phi += fPhi2 * pol2;
180 }
181
182 if(electronPol3)
183 {
184 // part depending on the polarization of the final electron
185 phi += fPhi3 * pol3;
186 }
187
188 // return cross section.
189 return phi;
190 }
191
192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193 G4double G4PolarizedComptonXS::TotalXSection(G4double /*xmin*/,
194 G4double /*xmax*/, G4double k0,
195 const G4StokesVector& pol0,
196 const G4StokesVector& pol1)
197 {
198 G4double k1 = 1. + 2. * k0;
199
200 G4double unit = fZ * CLHEP::pi * CLHEP::classic_electr_radius *
201 CLHEP::classic_electr_radius;
202
203 G4double pre = unit / (sqr(k0) * sqr(1. + 2. * k0));
204
205 G4double xs_0 = ((k0 - 2.) * k0 - 2.) * sqr(k1) * std::log(k1) +
206 2. * k0 * (k0 * (k0 + 1.) * (k0 + 8.) + 2.);
207 G4double xs_pol = (k0 + 1.) * sqr(k1) * std::log(k1) -
208 2. * k0 * (5. * sqr(k0) + 4. * k0 + 1.);
209
210 return pre * (xs_0 / k0 + pol0.p3() * pol1.z() * xs_pol);
211 }
212
213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
214 G4StokesVector G4PolarizedComptonXS::GetPol2()
215 {
216 // Note, mean polarization can not contain correlation effects.
217 return G4StokesVector(1. / fPhi0 * fPhi2);
218 }
219
220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
221 G4StokesVector G4PolarizedComptonXS::GetPol3()
222 {
223 // Note, mean polarization can not contain correlation effects.
224 return G4StokesVector(1. / fPhi0 * fPhi3);
225 }
226
227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
228 void G4PolarizedComptonXS::DefineCoefficients(const G4StokesVector& pol0,
229 const G4StokesVector& pol1)
230 {
231 polxx = pol0.x() * pol1.x();
232 polyy = pol0.y() * pol1.y();
233 polzz = pol0.z() * pol1.z();
234
235 polxz = pol0.x() * pol1.z();
236 polzx = pol0.z() * pol1.x();
237
238 polyz = pol0.y() * pol1.z();
239 polzy = pol0.z() * pol1.y();
240
241 polxy = pol0.x() * pol1.y();
242 polyx = pol0.y() * pol1.x();
243 }
244