| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Geant4 Class file 26 // Geant4 Class file
27 // 27 //
28 // File name: G4PolarizedComptonXS 28 // File name: G4PolarizedComptonXS
29 // 29 //
30 // Author: Andreas Schaelicke 30 // Author: Andreas Schaelicke
31 // 31 //
32 // Class Description: 32 // Class Description:
33 // determine the polarization of the final 33 // determine the polarization of the final state in a Compton scattering
34 // process employing the differential cross 34 // process employing the differential cross section by F.W.Lipps & H.A.Tolhoek
35 // ( Physica 20 (1954) 395 ) 35 // ( Physica 20 (1954) 395 )
36 // recalculated by P.Starovoitov 36 // recalculated by P.Starovoitov
37 37
38 #include "G4PolarizedComptonXS.hh" 38 #include "G4PolarizedComptonXS.hh"
39 39
40 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41 G4PolarizedComptonXS::G4PolarizedComptonXS() 41 G4PolarizedComptonXS::G4PolarizedComptonXS()
42 { 42 {
43 SetYmin(0.); 43 SetYmin(0.);
44 44
45 fPhi0 = 0.; 45 fPhi0 = 0.;
46 fPolXS = 0.; 46 fPolXS = 0.;
47 fUnpXS = 0.; 47 fUnpXS = 0.;
48 fPhi2 = G4ThreeVector(0., 0., 0.); 48 fPhi2 = G4ThreeVector(0., 0., 0.);
49 fPhi3 = G4ThreeVector(0., 0., 0.); 49 fPhi3 = G4ThreeVector(0., 0., 0.);
50 polxx = polyy = polzz = polxz = polzx = poly 50 polxx = polyy = polzz = polxz = polzx = polyz = polzy = polxy = polyx = 0.;
51 } 51 }
52 52
53 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 G4PolarizedComptonXS::~G4PolarizedComptonXS() 54 G4PolarizedComptonXS::~G4PolarizedComptonXS() {}
55 55
56 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 void G4PolarizedComptonXS::Initialize(G4double 57 void G4PolarizedComptonXS::Initialize(G4double eps, G4double X,
58 G4double 58 G4double, // phi
59 const G4 59 const G4StokesVector& pol0,
60 const G4 60 const G4StokesVector& pol1, G4int flag)
61 { 61 {
62 G4double cosT = 1. - (1. / eps - 1.) / X; 62 G4double cosT = 1. - (1. / eps - 1.) / X;
63 if(cosT > 1. + 1.e-8) 63 if(cosT > 1. + 1.e-8)
64 cosT = 1.; 64 cosT = 1.;
65 else if(cosT < -1. - 1.e-8) 65 else if(cosT < -1. - 1.e-8)
66 cosT = -1.; 66 cosT = -1.;
67 G4double cosT2 = cosT * cosT; 67 G4double cosT2 = cosT * cosT;
68 G4double cosT3 = cosT2 * cosT; 68 G4double cosT3 = cosT2 * cosT;
69 G4double sinT2 = 1. - cosT2; 69 G4double sinT2 = 1. - cosT2;
70 if(sinT2 > 1. + 1.e-8) 70 if(sinT2 > 1. + 1.e-8)
71 sinT2 = 1.; 71 sinT2 = 1.;
72 else if(sinT2 < 0.) 72 else if(sinT2 < 0.)
73 sinT2 = 0.; 73 sinT2 = 0.;
74 G4double sinT = std::sqrt(sinT2); 74 G4double sinT = std::sqrt(sinT2);
75 G4double cos2T = 2. * cosT2 - 1.; 75 G4double cos2T = 2. * cosT2 - 1.;
76 G4double sin2T = 2. * sinT * cosT; 76 G4double sin2T = 2. * sinT * cosT;
77 G4double eps2 = sqr(eps); 77 G4double eps2 = sqr(eps);
78 DefineCoefficients(pol0, pol1); 78 DefineCoefficients(pol0, pol1);
79 G4double diffXSFactor = re2 / (4. * X); 79 G4double diffXSFactor = re2 / (4. * X);
80 80
81 // unpolarized Cross Section 81 // unpolarized Cross Section
82 fUnpXS = (eps2 + 1. - eps * sinT2) / (2. * e 82 fUnpXS = (eps2 + 1. - eps * sinT2) / (2. * eps);
83 // initial polarization dependence 83 // initial polarization dependence
84 fPolXS = -sinT2 * pol0.x() + (1. - eps) * si 84 fPolXS = -sinT2 * pol0.x() + (1. - eps) * sinT * polzx +
85 ((eps2 - 1.) / eps) * cosT * polzz; 85 ((eps2 - 1.) / eps) * cosT * polzz;
86 fPolXS *= 0.5; 86 fPolXS *= 0.5;
87 87
88 fPhi0 = fUnpXS + fPolXS; 88 fPhi0 = fUnpXS + fPolXS;
89 89
90 if(flag == 2) 90 if(flag == 2)
91 { 91 {
92 // polarization of outgoing photon 92 // polarization of outgoing photon
93 G4double phi21 = -sinT2 + 0.5 * (cos2T + 3 93 G4double phi21 = -sinT2 + 0.5 * (cos2T + 3.) * pol0.x() -
94 ((1. - eps) / eps) * sinT 94 ((1. - eps) / eps) * sinT * polzx;
95 phi21 *= 0.5; 95 phi21 *= 0.5;
96 G4double phi22 = cosT * pol0.y() + ((1. - 96 G4double phi22 = cosT * pol0.y() + ((1. - eps) / (2. * eps)) * sinT * polzy;
97 G4double phi23 = ((eps2 + 1.) / eps) * cos 97 G4double phi23 = ((eps2 + 1.) / eps) * cosT * pol0.z() -
98 ((1. - eps) / eps) * (eps 98 ((1. - eps) / eps) * (eps * cosT2 + 1.) * pol1.z();
99 phi23 += 0.5 * (1. - eps) * sin2T * pol1.x 99 phi23 += 0.5 * (1. - eps) * sin2T * pol1.x();
100 phi23 += (eps - 1.) * (-sinT2 * polxz + si 100 phi23 += (eps - 1.) * (-sinT2 * polxz + sinT * polyy - 0.5 * sin2T * polxx);
101 phi23 *= 0.5; 101 phi23 *= 0.5;
102 102
103 fPhi2 = G4ThreeVector(phi21, phi22, phi23) 103 fPhi2 = G4ThreeVector(phi21, phi22, phi23);
104 104
105 // polarization of outgoing electron 105 // polarization of outgoing electron
106 G4double phi32 = -sinT2 * polxy + ((1. - e 106 G4double phi32 = -sinT2 * polxy + ((1. - eps) / eps) * sinT * polyz +
107 0.5 * (cos2T + 3.) * pol1 107 0.5 * (cos2T + 3.) * pol1.y();
108 phi32 *= 0.5; 108 phi32 *= 0.5;
109 109
110 G4double phi31 = 0.; 110 G4double phi31 = 0.;
111 G4double phi31add = 0.; 111 G4double phi31add = 0.;
112 G4double phi33 = 0.; 112 G4double phi33 = 0.;
113 G4double phi33add = 0.; 113 G4double phi33add = 0.;
114 114
115 if((1. - eps) > 1.e-12) 115 if((1. - eps) > 1.e-12)
116 { 116 {
117 G4double helpVar = std::sqrt(eps2 - 2. * 117 G4double helpVar = std::sqrt(eps2 - 2. * cosT * eps + 1.);
118 118
119 phi31 = (1. - eps) * (1. + cosT) * sinT 119 phi31 = (1. - eps) * (1. + cosT) * sinT * pol0.z();
120 phi31 += 120 phi31 +=
121 (-eps * cosT3 + eps * cosT2 + (eps - 2 121 (-eps * cosT3 + eps * cosT2 + (eps - 2.) * cosT + eps) * pol1.x();
122 phi31 += -(eps * cosT2 - eps * cosT + co 122 phi31 += -(eps * cosT2 - eps * cosT + cosT + 1.) * sinT * pol1.z();
123 phi31 /= 2. * helpVar; 123 phi31 /= 2. * helpVar;
124 124
125 phi31add = -eps * sqr(1. - cosT) * (1. + 125 phi31add = -eps * sqr(1. - cosT) * (1. + cosT) * polxx;
126 phi31add += (1. - eps) * sinT2 * polyy; 126 phi31add += (1. - eps) * sinT2 * polyy;
127 phi31add += -(-eps2 + cosT * (cosT * eps 127 phi31add += -(-eps2 + cosT * (cosT * eps - eps + 1.) * eps + eps - 1.) *
128 sinT * polxz / eps; 128 sinT * polxz / eps;
129 phi31add /= 2. * helpVar; 129 phi31add /= 2. * helpVar;
130 130
131 phi33 = ((1. - eps) / eps) * 131 phi33 = ((1. - eps) / eps) *
132 (-eps * cosT2 + eps * (eps + 1.) 132 (-eps * cosT2 + eps * (eps + 1.) * cosT - 1.) * pol0.z();
133 phi33 += -(eps * cosT2 + (1. - eps) * ep 133 phi33 += -(eps * cosT2 + (1. - eps) * eps * cosT + 1.) * sinT * pol1.x();
134 phi33 += 134 phi33 +=
135 -(-eps2 * cosT3 + eps * (eps2 - eps + 135 -(-eps2 * cosT3 + eps * (eps2 - eps + 1.) * cosT2 - cosT + eps2) *
136 pol1.z() / eps; 136 pol1.z() / eps;
137 phi33 /= -2. * helpVar; 137 phi33 /= -2. * helpVar;
138 138
139 phi33add = (eps * (eps - cosT - 1.) * co 139 phi33add = (eps * (eps - cosT - 1.) * cosT + 1.) * sinT * polxx;
140 phi33add += -(-eps2 + cosT * eps + eps - 140 phi33add += -(-eps2 + cosT * eps + eps - 1.) * sinT2 * polxz;
141 phi33add += (eps - 1.) * (cosT - eps) * 141 phi33add += (eps - 1.) * (cosT - eps) * sinT * polyy;
142 phi33add /= -2. * helpVar; 142 phi33add /= -2. * helpVar;
143 } 143 }
144 else 144 else
145 { 145 {
146 phi31 = -pol1.z() - 146 phi31 = -pol1.z() -
147 (X - 1.) * std::sqrt(1. - eps) * 147 (X - 1.) * std::sqrt(1. - eps) * pol1.x() / std::sqrt(2. * X);
148 phi31add = -(-X * X * pol1.z() - 2. * X 148 phi31add = -(-X * X * pol1.z() - 2. * X * (2. * pol0.z() - pol1.z()) -
149 (4. * pol0.x() + 5.) * pol1 149 (4. * pol0.x() + 5.) * pol1.z()) *
150 (1. - eps) / (4. * X); 150 (1. - eps) / (4. * X);
151 151
152 phi33 = pol1.x() - 152 phi33 = pol1.x() -
153 (X - 1.) * std::sqrt(1. - eps) * 153 (X - 1.) * std::sqrt(1. - eps) * pol1.z() / std::sqrt(2. * X);
154 phi33add = -(X * X - 2. * X + 4. * pol0. 154 phi33add = -(X * X - 2. * X + 4. * pol0.x() + 5.) * (1. - eps) *
155 pol1.x() / (4. * X); 155 pol1.x() / (4. * X);
156 } 156 }
157 fPhi3 = G4ThreeVector(phi31 + phi31add, ph 157 fPhi3 = G4ThreeVector(phi31 + phi31add, phi32, phi33 + phi33add);
158 } 158 }
159 fUnpXS *= diffXSFactor; 159 fUnpXS *= diffXSFactor;
160 fPolXS *= diffXSFactor; 160 fPolXS *= diffXSFactor;
161 fPhi0 *= diffXSFactor; 161 fPhi0 *= diffXSFactor;
162 fPhi2 *= diffXSFactor; 162 fPhi2 *= diffXSFactor;
163 fPhi3 *= diffXSFactor; 163 fPhi3 *= diffXSFactor;
164 } 164 }
165 165
166 G4double G4PolarizedComptonXS::XSection(const 166 G4double G4PolarizedComptonXS::XSection(const G4StokesVector& pol2,
167 const 167 const G4StokesVector& pol3)
168 { 168 {
169 G4bool gammaPol2 = !(pol2 == G4StokesVect 169 G4bool gammaPol2 = !(pol2 == G4StokesVector::ZERO);
170 G4bool electronPol3 = !(pol3 == G4StokesVect 170 G4bool electronPol3 = !(pol3 == G4StokesVector::ZERO);
171 171
172 G4double phi = 0.; 172 G4double phi = 0.;
173 // polarization independent part 173 // polarization independent part
174 phi += fPhi0; 174 phi += fPhi0;
175 175
176 if(gammaPol2) 176 if(gammaPol2)
177 { 177 {
178 // part depending on the polarization of t 178 // part depending on the polarization of the final photon
179 phi += fPhi2 * pol2; 179 phi += fPhi2 * pol2;
180 } 180 }
181 181
182 if(electronPol3) 182 if(electronPol3)
183 { 183 {
184 // part depending on the polarization of t 184 // part depending on the polarization of the final electron
185 phi += fPhi3 * pol3; 185 phi += fPhi3 * pol3;
186 } 186 }
187 187
188 // return cross section. 188 // return cross section.
189 return phi; 189 return phi;
190 } 190 }
191 191
192 //....oooOO0OOooo........oooOO0OOooo........oo 192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193 G4double G4PolarizedComptonXS::TotalXSection(G 193 G4double G4PolarizedComptonXS::TotalXSection(G4double /*xmin*/,
194 G 194 G4double /*xmax*/, G4double k0,
195 c 195 const G4StokesVector& pol0,
196 c 196 const G4StokesVector& pol1)
197 { 197 {
198 G4double k1 = 1. + 2. * k0; 198 G4double k1 = 1. + 2. * k0;
199 199
200 G4double unit = fZ * CLHEP::pi * CLHEP::clas 200 G4double unit = fZ * CLHEP::pi * CLHEP::classic_electr_radius *
201 CLHEP::classic_electr_radius 201 CLHEP::classic_electr_radius;
202 202
203 G4double pre = unit / (sqr(k0) * sqr(1. + 2. 203 G4double pre = unit / (sqr(k0) * sqr(1. + 2. * k0));
204 204
205 G4double xs_0 = ((k0 - 2.) * k0 - 2.) * sqr( 205 G4double xs_0 = ((k0 - 2.) * k0 - 2.) * sqr(k1) * std::log(k1) +
206 2. * k0 * (k0 * (k0 + 1.) * 206 2. * k0 * (k0 * (k0 + 1.) * (k0 + 8.) + 2.);
207 G4double xs_pol = (k0 + 1.) * sqr(k1) * std: 207 G4double xs_pol = (k0 + 1.) * sqr(k1) * std::log(k1) -
208 2. * k0 * (5. * sqr(k0) + 208 2. * k0 * (5. * sqr(k0) + 4. * k0 + 1.);
209 209
210 return pre * (xs_0 / k0 + pol0.p3() * pol1.z 210 return pre * (xs_0 / k0 + pol0.p3() * pol1.z() * xs_pol);
211 } 211 }
212 212
213 //....oooOO0OOooo........oooOO0OOooo........oo 213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
214 G4StokesVector G4PolarizedComptonXS::GetPol2() 214 G4StokesVector G4PolarizedComptonXS::GetPol2()
215 { 215 {
216 // Note, mean polarization can not contain c 216 // Note, mean polarization can not contain correlation effects.
217 return G4StokesVector(1. / fPhi0 * fPhi2); 217 return G4StokesVector(1. / fPhi0 * fPhi2);
218 } 218 }
219 219
220 //....oooOO0OOooo........oooOO0OOooo........oo 220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
221 G4StokesVector G4PolarizedComptonXS::GetPol3() 221 G4StokesVector G4PolarizedComptonXS::GetPol3()
222 { 222 {
223 // Note, mean polarization can not contain c 223 // Note, mean polarization can not contain correlation effects.
224 return G4StokesVector(1. / fPhi0 * fPhi3); 224 return G4StokesVector(1. / fPhi0 * fPhi3);
225 } 225 }
226 226
227 //....oooOO0OOooo........oooOO0OOooo........oo 227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
228 void G4PolarizedComptonXS::DefineCoefficients( 228 void G4PolarizedComptonXS::DefineCoefficients(const G4StokesVector& pol0,
229 229 const G4StokesVector& pol1)
230 { 230 {
231 polxx = pol0.x() * pol1.x(); 231 polxx = pol0.x() * pol1.x();
232 polyy = pol0.y() * pol1.y(); 232 polyy = pol0.y() * pol1.y();
233 polzz = pol0.z() * pol1.z(); 233 polzz = pol0.z() * pol1.z();
234 234
235 polxz = pol0.x() * pol1.z(); 235 polxz = pol0.x() * pol1.z();
236 polzx = pol0.z() * pol1.x(); 236 polzx = pol0.z() * pol1.x();
237 237
238 polyz = pol0.y() * pol1.z(); 238 polyz = pol0.y() * pol1.z();
239 polzy = pol0.z() * pol1.y(); 239 polzy = pol0.z() * pol1.y();
240 240
241 polxy = pol0.x() * pol1.y(); 241 polxy = pol0.x() * pol1.y();
242 polyx = pol0.y() * pol1.x(); 242 polyx = pol0.y() * pol1.x();
243 } 243 }
244 244