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 // ------------------------------------------------------------------- 27 // 28 // Geant4 Class file 29 // 30 // File name: G4PolarizedIonisationMollerXS 31 // 32 // Author: Andreas Schaelicke 33 // 34 // Class Description: 35 // * calculates the differential cross section 36 // incoming electron K1(along positive z direction) scatters at an electron 37 // K2 at rest 38 // * phi denotes the angle between the scattering plane (defined by the 39 // outgoing electron) and X-axis 40 // * all stokes vectors refer to spins in the Global System (X,Y,Z) 41 // * cross section as calculated by P.Starovoitov 42 43 #include "G4PolarizedIonisationMollerXS.hh" 44 45 #include "G4PhysicalConstants.hh" 46 47 G4PolarizedIonisationMollerXS::G4PolarizedIonisationMollerXS() 48 : fPhi0(0.) 49 { 50 SetXmax(.5); 51 fPhi2 = G4ThreeVector(0., 0., 0.); 52 fPhi3 = G4ThreeVector(0., 0., 0.); 53 } 54 55 G4PolarizedIonisationMollerXS::~G4PolarizedIonisationMollerXS() {} 56 57 void G4PolarizedIonisationMollerXS::Initialize(G4double e, G4double gamma, 58 G4double /*phi*/, 59 const G4StokesVector& pol0, 60 const G4StokesVector& pol1, 61 G4int flag) 62 { 63 constexpr G4double re2 = classic_electr_radius * classic_electr_radius; 64 G4double gamma2 = gamma * gamma; 65 G4double gmo = (gamma - 1.); 66 G4double gmo2 = (gamma - 1.) * (gamma - 1.); 67 G4double gpo = (gamma + 1.); 68 G4double pref = gamma2 * re2 / (gmo2 * (gamma + 1.0)); 69 constexpr G4double sqrttwo = 1.41421356237309504880; // sqrt(2.) 70 G4double f = (-1. + e); 71 G4double e2 = e * e; 72 G4double f2 = f * f; 73 74 G4bool polarized = (!pol0.IsZero()) || (!pol1.IsZero()); 75 76 if(flag == 0) 77 polarized = false; 78 // Unpolarised part of XS 79 fPhi0 = gmo2 / gamma2; 80 fPhi0 += ((1. - 2. * gamma) / gamma2) * (1. / e + 1. / (1. - e)); 81 fPhi0 += 1. / (e * e) + 1. / ((1. - e) * (1. - e)); 82 fPhi0 *= 0.25; 83 // Initial state polarisarion dependence 84 if(polarized) 85 { 86 G4bool usephi = true; 87 if(flag <= 1) 88 usephi = false; 89 G4double xx = (gamma - f * e * gmo * (3. + gamma)) / (4. * f * e * gamma2); 90 G4double yy = (-1. + f * e * gmo2 + 2. * gamma) / (4. * f * e * gamma2); 91 G4double zz = (-(e * gmo * (3. + gamma)) + e2 * gmo * (3. + gamma) + 92 gamma * (-1. + 2. * gamma)) / 93 (4. * f * e * gamma2); 94 95 fPhi0 += xx * pol0.x() * pol1.x() + yy * pol0.y() * pol1.y() + 96 zz * pol0.z() * pol1.z(); 97 98 if(usephi) 99 { 100 G4double xy = 0.; 101 G4double xz = -((-1. + 2. * e) * gmo) / 102 (2. * sqrttwo * gamma2 * std::sqrt(-((f * e) / gpo))); 103 G4double yx = 0.; 104 G4double yz = 0.; 105 G4double zx = -((-1. + 2. * e) * gmo) / 106 (2. * sqrttwo * gamma2 * std::sqrt(-((f * e) / gpo))); 107 G4double zy = 0.; 108 fPhi0 += yx * pol0.y() * pol1.x() + xy * pol0.x() * pol1.y(); 109 fPhi0 += zx * pol0.z() * pol1.x() + xz * pol0.x() * pol1.z(); 110 fPhi0 += zy * pol0.z() * pol1.y() + yz * pol0.y() * pol1.z(); 111 } 112 } 113 // Final state polarisarion dependence 114 fPhi2 = G4ThreeVector(); 115 fPhi3 = G4ThreeVector(); 116 117 if(flag >= 1) 118 { 119 // Final Electron P1 120 // initial electron K1 121 if(!pol0.IsZero()) 122 { 123 G4double xxP1K1 = 124 (std::sqrt(gpo / (1. + e2 * gmo + gamma - 2. * e * gamma)) * 125 (gamma - e * gpo)) / 126 (4. * e2 * gamma); 127 G4double xyP1K1 = 0.; 128 G4double xzP1K1 = (-1. + 2. * e * gamma) / 129 (2. * sqrttwo * f * gamma * 130 std::sqrt(e * e2 * (1. + e + gamma - e * gamma))); 131 G4double yxP1K1 = 0.; 132 G4double yyP1K1 = 133 (-gamma2 + e * (-1. + gamma * (2. + gamma))) / (4. * f * e2 * gamma2); 134 G4double yzP1K1 = 0.; 135 G4double zxP1K1 = (1. + 2. * e2 * gmo - 2. * e * gamma) / 136 (2. * sqrttwo * f * e * gamma * 137 std::sqrt(e * (1. + e + gamma - e * gamma))); 138 G4double zyP1K1 = 0.; 139 G4double zzP1K1 = 140 (-gamma + e * (1. - 2. * e * gmo + gamma)) / 141 (4. * f * e2 * gamma * std::sqrt(1. - (2. * e) / (f * gpo))); 142 fPhi2[0] += xxP1K1 * pol0.x() + xyP1K1 * pol0.y() + xzP1K1 * pol0.z(); 143 fPhi2[1] += yxP1K1 * pol0.x() + yyP1K1 * pol0.y() + yzP1K1 * pol0.z(); 144 fPhi2[2] += zxP1K1 * pol0.x() + zyP1K1 * pol0.y() + zzP1K1 * pol0.z(); 145 } 146 // initial electron K2 147 if(!pol1.IsZero()) 148 { 149 G4double xxP1K2 = 150 ((1. + e * (-3. + gamma)) * 151 std::sqrt(gpo / (1. + e2 * gmo + gamma - 2. * e * gamma))) / 152 (4. * f * e * gamma); 153 G4double xyP1K2 = 0.; 154 G4double xzP1K2 = 155 (-2. + 2. * e + gamma) / (2. * sqrttwo * f2 * gamma * 156 std::sqrt(e * (1. + e + gamma - e * gamma))); 157 G4double yxP1K2 = 0.; 158 G4double yyP1K2 = (1. - 2. * gamma + e * (-1. + gamma * (2. + gamma))) / 159 (4. * f2 * e * gamma2); 160 G4double yzP1K2 = 0.; 161 G4double zxP1K2 = (2. * e * (1. + e * gmo - 2. * gamma) + gamma) / 162 (2. * sqrttwo * f2 * gamma * 163 std::sqrt(e * (1. + e + gamma - e * gamma))); 164 G4double zyP1K2 = 0.; 165 G4double zzP1K2 = 166 (1. - 2. * gamma + e * (-1. - 2. * e * gmo + 3. * gamma)) / 167 (4. * f2 * e * gamma * std::sqrt(1. - (2. * e) / (f * gpo))); 168 fPhi2[0] += xxP1K2 * pol1.x() + xyP1K2 * pol1.y() + xzP1K2 * pol1.z(); 169 fPhi2[1] += yxP1K2 * pol1.x() + yyP1K2 * pol1.y() + yzP1K2 * pol1.z(); 170 fPhi2[2] += zxP1K2 * pol1.x() + zyP1K2 * pol1.y() + zzP1K2 * pol1.z(); 171 } 172 173 // Final Electron P2 174 // initial electron K1 175 if(!pol0.IsZero()) 176 { 177 G4double xxP2K1 = 178 (-1. + e + e * gamma) / 179 (4. * f2 * gamma * std::sqrt((e * (2. + e * gmo)) / gpo)); 180 G4double xyP2K1 = 0.; 181 G4double xzP2K1 = 182 -((1. + 2. * f * gamma) * std::sqrt(f / (-2. + e - e * gamma))) / 183 (2. * sqrttwo * f2 * e * gamma); 184 G4double yxP2K1 = 0.; 185 G4double yyP2K1 = (1. - 2. * gamma + e * (-1. + gamma * (2. + gamma))) / 186 (4. * f2 * e * gamma2); 187 G4double yzP2K1 = 0.; 188 G4double zxP2K1 = 189 (1. + 2. * e * (-2. + e + gamma - e * gamma)) / 190 (2. * sqrttwo * f * e * std::sqrt(-(f * (2. + e * gmo))) * gamma); 191 G4double zyP2K1 = 0.; 192 G4double zzP2K1 = 193 (std::sqrt((e * gpo) / (2. + e * gmo)) * 194 (-3. + e * (5. + 2. * e * gmo - 3. * gamma) + 2. * gamma)) / 195 (4. * f2 * e * gamma); 196 197 fPhi3[0] += xxP2K1 * pol0.x() + xyP2K1 * pol0.y() + xzP2K1 * pol0.z(); 198 fPhi3[1] += yxP2K1 * pol0.x() + yyP2K1 * pol0.y() + yzP2K1 * pol0.z(); 199 fPhi3[2] += zxP2K1 * pol0.x() + zyP2K1 * pol0.y() + zzP2K1 * pol0.z(); 200 } 201 // initial electron K2 202 if(!pol1.IsZero()) 203 { 204 G4double xxP2K2 = 205 (-2. - e * (-3. + gamma) + gamma) / 206 (4. * f * e * gamma * std::sqrt((e * (2. + e * gmo)) / gpo)); 207 G4double xyP2K2 = 0.; 208 G4double xzP2K2 = 209 ((-2. * e + gamma) * std::sqrt(f / (-2. + e - e * gamma))) / 210 (2. * sqrttwo * f * e2 * gamma); 211 G4double yxP2K2 = 0.; 212 G4double yyP2K2 = 213 (-gamma2 + e * (-1. + gamma * (2. + gamma))) / (4. * f * e2 * gamma2); 214 G4double yzP2K2 = 0.; 215 G4double zxP2K2 = 216 (gamma + 2. * e * (-1. + e - e * gamma)) / 217 (2. * sqrttwo * e2 * std::sqrt(-(f * (2. + e * gmo))) * gamma); 218 G4double zyP2K2 = 0.; 219 G4double zzP2K2 = (std::sqrt((e * gpo) / (2. + e * gmo)) * 220 (-2. + e * (3. + 2. * e * gmo - gamma) + gamma)) / 221 (4. * f * e2 * gamma); 222 fPhi3[0] += xxP2K2 * pol1.x() + xyP2K2 * pol1.y() + xzP2K2 * pol1.z(); 223 fPhi3[1] += yxP2K2 * pol1.x() + yyP2K2 * pol1.y() + yzP2K2 * pol1.z(); 224 fPhi3[2] += zxP2K2 * pol1.x() + zyP2K2 * pol1.y() + zzP2K2 * pol1.z(); 225 } 226 } 227 fPhi0 *= pref; 228 fPhi2 *= pref; 229 fPhi3 *= pref; 230 } 231 232 G4double G4PolarizedIonisationMollerXS::XSection(const G4StokesVector& pol2, 233 const G4StokesVector& pol3) 234 { 235 G4double xs = fPhi0; 236 237 G4bool polarized = (!pol2.IsZero()) || (!pol3.IsZero()); 238 if(polarized) 239 { 240 xs += fPhi2 * pol2 + fPhi3 * pol3; 241 } 242 return xs; 243 } 244 245 G4double G4PolarizedIonisationMollerXS::TotalXSection( 246 G4double xmin, G4double xmax, G4double gamma, const G4StokesVector& pol0, 247 const G4StokesVector& pol1) 248 { 249 G4double xs = 0.; 250 G4double x = xmin; 251 252 if(xmax != 0.5) 253 { 254 G4ExceptionDescription ed; 255 ed << " warning xmax expected to be 1/2 but is " << xmax << "\n"; 256 G4Exception("G4PolarizedIonisationMollerXS::TotalXSection", "pol020", 257 JustWarning, ed); 258 } 259 260 constexpr G4double re2 = classic_electr_radius * classic_electr_radius; 261 G4double gamma2 = gamma * gamma; 262 G4double gmo2 = (gamma - 1.) * (gamma - 1.); 263 G4double logMEM = std::log(1. / x - 1.); 264 G4double pref = twopi * gamma2 * re2 / (gmo2 * (gamma + 1.0)); 265 // unpolarised XS 266 G4double sigma0 = (gmo2 / gamma2) * (0.5 - x); 267 sigma0 += ((1. - 2. * gamma) / gamma2) * logMEM; 268 sigma0 += 1. / x - 1. / (1. - x); 269 // longitudinal part 270 G4double sigma2 = ((gamma2 + 2. * gamma - 3.) / gamma2) * (0.5 - x); 271 sigma2 += (1. / gamma - 2.) * logMEM; 272 // transverse part 273 G4double sigma3 = (2. * (1. - gamma) / gamma2) * (0.5 - x); 274 sigma3 += (1. - 3. * gamma) / (2. * gamma2) * logMEM; 275 // total cross section 276 xs += pref * (sigma0 + sigma2 * pol0.z() * pol1.z() + 277 sigma3 * (pol0.x() * pol1.x() + pol0.y() * pol1.y())); 278 279 return xs; 280 } 281 282 G4StokesVector G4PolarizedIonisationMollerXS::GetPol2() 283 { 284 // Note, mean polarization can not contain correlation effects. 285 return G4StokesVector(1. / fPhi0 * fPhi2); 286 } 287 G4StokesVector G4PolarizedIonisationMollerXS::GetPol3() 288 { 289 // Note, mean polarization can not contain correlation effects. 290 return G4StokesVector(1. / fPhi0 * fPhi3); 291 } 292