Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // Geant4 Class file 27 // 28 // File name: G4PolarizedAnnihilationXS 29 // 30 // Author: Andreas Schaelicke 31 // 32 // Class Description: 33 // * calculates the differential cross secti 34 35 #include "G4PolarizedAnnihilationXS.hh" 36 37 #include "G4PhysicalConstants.hh" 38 39 G4PolarizedAnnihilationXS::G4PolarizedAnnihila 40 : polxx(0.) 41 , polyy(0.) 42 , polzz(0.) 43 , polxz(0.) 44 , polzx(0.) 45 , polxy(0.) 46 , polyx(0.) 47 , polyz(0.) 48 , polzy(0.) 49 , fPhi0(0.) 50 { 51 fPhi2 = G4ThreeVector(0., 0., 0.); 52 fPhi3 = G4ThreeVector(0., 0., 0.); 53 fDice = 0.; 54 fPolXS = 0.; 55 fUnpXS = 0.; 56 ISPxx = ISPyy = ISPzz = ISPnd = 0.; 57 } 58 59 //....oooOO0OOooo........oooOO0OOooo........oo 60 G4PolarizedAnnihilationXS::~G4PolarizedAnnihil 61 62 //....oooOO0OOooo........oooOO0OOooo........oo 63 void G4PolarizedAnnihilationXS::TotalXS() 64 { 65 // total cross section is sum of 66 // + unpol xsec "sigma0" 67 // + longitudinal polarised cross section " 68 // pol_3(e-)*pol_3(e+) 69 // + transverse contribution "(sigma_xx+sig 70 // pol_T(e-)*pol_T(e+) 71 // (Note: if both beams are transversely 72 // pol_T(e+)!=0, and sigma_xx!=sigma_yy 73 // will exhibit a azimuthal asymmetry e 74 } 75 76 //....oooOO0OOooo........oooOO0OOooo........oo 77 void G4PolarizedAnnihilationXS::Initialize( 78 G4double eps, G4double gam, 79 G4double, // phi 80 const G4StokesVector& pol0, // positron pol 81 const G4StokesVector& pol1, // electron pol 82 G4int flag) 83 { 84 G4double diffXSFactor = re2 / (gam - 1.); 85 DefineCoefficients(pol0, pol1); 86 87 // prepare dicing 88 fDice = 0.; 89 G4double symmXS = 90 0.125 * ((-1. / sqr(gam + 1.)) / sqr(eps) 91 ((sqr(gam) + 4. * gam - 1.) / sqr 92 93 G4ThreeVector epsVector(1. / sqr(eps), 1. / 94 G4ThreeVector oneEpsVector(1. / sqr(1. - eps 95 G4ThreeVector sumEpsVector(epsVector + oneEp 96 G4ThreeVector difEpsVector(epsVector - oneEp 97 G4ThreeVector calcVector(0., 0., 0.); 98 99 // temporary variables 100 G4double helpVar2 = 0., helpVar1 = 0.; 101 102 // unpolarised contribution 103 helpVar1 = (gam * gam + 4. * gam + 1.) / s 104 helpVar2 = -1. / sqr(gam + 1.); 105 calcVector = G4ThreeVector(helpVar2, helpVar 106 fUnpXS = 0.125 * calcVector * sumEpsVect 107 108 // initial particles polarised contribution 109 helpVar2 = 1. / sqr(gam + 1.); 110 helpVar1 = -(gam * gam + 4. * gam + 1.) / 111 calcVector = G4ThreeVector(helpVar2, helpVar 112 ISPxx = 0.25 * (calcVector * sumEpsVect 113 114 helpVar1 = 1. / sqr(gam + 1.); 115 calcVector = G4ThreeVector(-helpVar1, 2. * g 116 ISPyy = 0.125 * calcVector * sumEpsVect 117 118 helpVar1 = 1. / (gam - 1.); 119 helpVar2 = 1. / sqr(gam + 1.); 120 calcVector = G4ThreeVector( 121 -(gam * gam + 1.) * helpVar2, 122 (gam * gam * (gam + 1.) + 7. * gam + 3.) * 123 ISPzz = 0.125 * helpVar1 * (calcVector * sum 124 125 helpVar1 = std::sqrt(std::fabs(eps * (1. - 126 calcVector = G4ThreeVector(-1. / (gam * gam 127 ISPnd = 0.125 * (calcVector * difEpsVec 128 129 fPolXS = ISPxx * polxx + ISPyy * polyy + ISP 130 fPolXS += ISPnd * (polzx + polxz); 131 fPhi0 = fUnpXS + fPolXS; 132 fDice = symmXS; 133 134 if(polzz != 0.) 135 { 136 fDice *= (1. + (polzz * ISPzz / fUnpXS)); 137 if(fDice < 0.) 138 fDice = 0.; 139 } 140 // prepare final state coefficients 141 if(flag == 2) 142 { 143 // circular polarisation 144 G4double circ1 = 0., circ2 = 0., circ3 = 0 145 helpVar1 = 8. * sqr(1. - eps) * sqr(eps) * 146 std::sqrt(gam * gam - 1.); 147 helpVar2 = sqr(gam + 1.) * sqr(eps) * (-2. 148 (gam * gam + 3. * gam + 2.) * e 149 circ1 = helpVar2 + gam; 150 circ1 /= helpVar1; 151 circ2 = helpVar2 + 1.; 152 circ2 /= helpVar1; 153 helpVar1 = std::sqrt(std::fabs(eps * (1. - 154 helpVar1 /= std::sqrt(gam * gam - 1.); 155 calcVector = G4ThreeVector(1., -2. * gam, 156 circ3 = 0.125 * (calcVector * sumEpsV 157 circ3 *= helpVar1; 158 159 fPhi2.setZ(circ2 * pol1.z() + circ1 * pol0 160 circ3 * (pol1.x() + pol0.x())); 161 fPhi3.setZ(-circ1 * pol1.z() - circ2 * pol 162 circ3 * (pol1.x() + pol0.x())); 163 164 // common to both linear polarisation 165 calcVector = G4ThreeVector(-1., 2 166 G4double linearZero = 0.125 * (calcVector 167 168 // Linear Polarisation #1 169 helpVar1 = std::sqrt(std::fabs(2. * (gam + 170 ((gam + 1.) * eps * (1. - eps)) 171 helpVar2 = helpVar1 * helpVar1; 172 173 // photon 1 174 G4double diagContrib = 0.125 * helpVar2 * 175 G4double nonDiagContrib = 176 0.125 * helpVar1 * (-polxz / (1. - eps) 177 178 fPhi2.setX(linearZero + diagContrib + nonD 179 180 // photon 2 181 nonDiagContrib = 0.125 * helpVar1 * (polxz 182 183 fPhi3.setX(linearZero + diagContrib + nonD 184 185 // Linear Polarisation #2 186 helpVar1 = 187 std::sqrt(gam * gam - 1.) * (2. * (gam + 188 helpVar1 /= 8. * sqr(1. - eps) * sqr(eps) 189 helpVar2 = std::sqrt((gam * gam - 1.) * 190 std::fabs(2. * (gam + 191 helpVar2 /= 8. * sqr(1. - eps) * sqr(eps) 192 193 G4double contrib21 = (-polxy + polyx) * he 194 G4double contrib32 = 195 -(eps * (gam + 1.) - 1.) * polyz + (eps 196 197 contrib32 *= helpVar2; 198 fPhi2.setY(contrib21 + contrib32); 199 200 contrib32 = 201 -(eps * (gam + 1.) - gam) * polyz + (eps 202 contrib32 *= helpVar2; 203 fPhi3.setY(contrib21 + contrib32); 204 } 205 fPhi0 *= diffXSFactor; 206 fPhi2 *= diffXSFactor; 207 fPhi3 *= diffXSFactor; 208 } 209 210 //....oooOO0OOooo........oooOO0OOooo........oo 211 G4double G4PolarizedAnnihilationXS::XSection(c 212 c 213 { 214 G4double xs = fPhi0 + pol2 * fPhi2 + pol3 * 215 return xs; 216 } 217 218 // calculate total cross section 219 //....oooOO0OOooo........oooOO0OOooo........oo 220 G4double G4PolarizedAnnihilationXS::TotalXSect 221 222 223 224 { 225 G4double totalXSFactor = pi * re2 / (gam + 1 226 DefineCoefficients(pol0, pol1); 227 228 G4double xs = 0.; 229 230 G4double gam2 = gam * gam; 231 G4double sqrtgam1 = std::sqrt(gam2 - 1.); 232 G4double logMEM = std::log(gam + sqrtgam1) 233 G4double unpME = (gam * (gam + 4.) + 1.) 234 unpME += -(gam + 3.) * sqrtgam1; 235 unpME /= 4. * (gam2 - 1.); 236 G4double longPart = (3 + gam * (gam * (gam + 237 longPart += -(5. + gam * (3 * gam + 4.)) * s 238 longPart /= 4. * sqr(gam - 1.) * (gam + 1.); 239 G4double tranPart = -(5 * gam + 1.) * logMEM 240 tranPart += (gam + 5.) * sqrtgam1; 241 tranPart /= 4. * sqr(gam - 1.) * (gam + 1.); 242 243 xs += unpME; 244 xs += polzz * longPart; 245 xs += (polxx + polyy) * tranPart; 246 247 return xs * totalXSFactor; 248 } 249 250 //....oooOO0OOooo........oooOO0OOooo........oo 251 G4StokesVector G4PolarizedAnnihilationXS::GetP 252 { 253 // Note, mean polarization can not contain c 254 return G4StokesVector(1. / fPhi0 * fPhi2); 255 } 256 257 //....oooOO0OOooo........oooOO0OOooo........oo 258 G4StokesVector G4PolarizedAnnihilationXS::GetP 259 { 260 // Note, mean polarization can not contain c 261 return G4StokesVector(1. / fPhi0 * fPhi3); 262 } 263 //....oooOO0OOooo........oooOO0OOooo........oo 264 void G4PolarizedAnnihilationXS::DefineCoeffici 265 266 { 267 polxx = pol0.x() * pol1.x(); 268 polyy = pol0.y() * pol1.y(); 269 polzz = pol0.z() * pol1.z(); 270 271 polxz = pol0.x() * pol1.z(); 272 polzx = pol0.z() * pol1.x(); 273 274 polyz = pol0.y() * pol1.z(); 275 polzy = pol0.z() * pol1.y(); 276 277 polxy = pol0.x() * pol1.y(); 278 polyx = pol0.y() * pol1.x(); 279 } 280 281 //....oooOO0OOooo........oooOO0OOooo........oo 282 G4double G4PolarizedAnnihilationXS::GetXmin(G4 283 { 284 return 0.5 * (1. - std::sqrt((y - 1.) / (y + 285 } 286 287 //....oooOO0OOooo........oooOO0OOooo........oo 288 G4double G4PolarizedAnnihilationXS::GetXmax(G4 289 { 290 return 0.5 * (1. + std::sqrt((y - 1.) / (y + 291 } 292 293 //....oooOO0OOooo........oooOO0OOooo........oo 294 G4double G4PolarizedAnnihilationXS::DiceEpsilo 295 296 //....oooOO0OOooo........oooOO0OOooo........oo 297 G4double G4PolarizedAnnihilationXS::getVar(G4i 298 { 299 if(choice == -1) 300 return fPolXS / fUnpXS; 301 if(choice == 0) 302 return fUnpXS; 303 if(choice == 1) 304 return ISPxx; 305 if(choice == 2) 306 return ISPyy; 307 if(choice == 3) 308 return ISPzz; 309 if(choice == 4) 310 return ISPnd; 311 return 0; 312 } 313