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