00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 #include "G4PolarizedMollerCrossSection.hh"
00050 #include "G4PhysicalConstants.hh"
00051
00052 G4PolarizedMollerCrossSection::G4PolarizedMollerCrossSection() :
00053 phi0(0.)
00054 {
00055 SetXmax(.5);
00056 }
00057 G4PolarizedMollerCrossSection::~G4PolarizedMollerCrossSection() {}
00058 void G4PolarizedMollerCrossSection::Initialize(
00059 G4double e,
00060 G4double gamma,
00061 G4double ,
00062 const G4StokesVector & pol0,
00063 const G4StokesVector & pol1,
00064 G4int flag)
00065 {
00066 G4double re2 = classic_electr_radius * classic_electr_radius;
00067 G4double gamma2=gamma*gamma;
00068 G4double gmo = (gamma - 1.);
00069 G4double gmo2 = (gamma - 1.)*(gamma - 1.);
00070 G4double gpo = (gamma + 1.);
00071 G4double pref = gamma2*re2/(gmo2*(gamma + 1.0));
00072 G4double sqrttwo=std::sqrt(2.);
00073 G4double f = (-1. + e);
00074 G4double e2 = e*e;
00075 G4double f2 = f*f;
00076
00077
00078 G4bool polarized=(!pol0.IsZero())||(!pol1.IsZero());
00079
00080 if (flag==0) polarized=false;
00081
00082 phi0 = 0.;
00083 phi0+= gmo2/gamma2;
00084 phi0+= ((1. - 2.*gamma)/gamma2)*(1./e + 1./(1.-e));
00085 phi0+= 1./(e*e) + 1./((1. - e)*(1. - e));
00086 phi0*=0.25;
00087
00088 if (polarized) {
00089 G4double usephi=1.;
00090 if (flag<=1) usephi=0.;
00091
00092
00093
00094 G4double xx = (gamma - f*e*gmo*(3 + gamma))/(4*f*e*gamma2);
00095 G4double yy = (-1 + f*e*gmo2 + 2*gamma)/(4*f*e*gamma2);
00096 G4double zz = (-(e*gmo*(3 + gamma)) + e2*gmo*(3 + gamma) +
00097 gamma*(-1 + 2*gamma))/(4*f*e*gamma2);
00098
00099 phi0 += xx*pol0.x()*pol1.x() + yy*pol0.y()*pol1.y() + zz*pol0.z()*pol1.z();
00100
00101 if (usephi==1.) {
00102
00103 G4double xy = 0;
00104 G4double xz = -((-1 + 2*e)*gmo)/(2*sqrttwo*gamma2*
00105 std::sqrt(-((f*e)/gpo)));
00106 G4double yx = 0;
00107 G4double yz = 0;
00108 G4double zx = -((-1 + 2*e)*gmo)/(2*sqrttwo*gamma2*
00109 std::sqrt(-((f*e)/gpo)));
00110 G4double zy = 0;
00111 phi0+=yx*pol0.y()*pol1.x() + xy*pol0.x()*pol1.y();
00112 phi0+=zx*pol0.z()*pol1.x() + xz*pol0.x()*pol1.z();
00113 phi0+=zy*pol0.z()*pol1.y() + yz*pol0.y()*pol1.z();
00114 }
00115 }
00116
00117 phi2=G4ThreeVector();
00118 phi3=G4ThreeVector();
00119
00120 if (flag>=1) {
00121
00122
00123
00124
00125
00126 if (!pol0.IsZero()) {
00127 G4double xxP1K1 = (std::sqrt(gpo/(1 + e2*gmo + gamma - 2*e*gamma))*
00128 (gamma - e*gpo))/(4*e2*gamma);
00129 G4double xyP1K1 = 0;
00130 G4double xzP1K1 = (-1 + 2*e*gamma)/(2*sqrttwo*f*gamma*
00131 std::sqrt(e*e2*(1 + e + gamma - e*gamma)));
00132 G4double yxP1K1 = 0;
00133 G4double yyP1K1 = (-gamma2 + e*(-1 + gamma*(2 + gamma)))/(4*f*e2*gamma2);
00134 G4double yzP1K1 = 0;
00135 G4double zxP1K1 = (1 + 2*e2*gmo - 2*e*gamma)/(2*sqrttwo*f*e*gamma*
00136 std::sqrt(e*(1 + e + gamma - e*gamma)));
00137 G4double zyP1K1 = 0;
00138 G4double zzP1K1 = (-gamma + e*(1 - 2*e*gmo + gamma))/(4*f*e2*gamma*
00139 std::sqrt(1 - (2*e)/(f*gpo)));
00140 phi2[0] += xxP1K1*pol0.x() + xyP1K1*pol0.y() + xzP1K1*pol0.z();
00141 phi2[1] += yxP1K1*pol0.x() + yyP1K1*pol0.y() + yzP1K1*pol0.z();
00142 phi2[2] += zxP1K1*pol0.x() + zyP1K1*pol0.y() + zzP1K1*pol0.z();
00143 }
00144
00145 if (!pol1.IsZero()) {
00146 G4double xxP1K2 = ((1 + e*(-3 + gamma))*std::sqrt(gpo/(1 + e2*gmo + gamma -
00147 2*e*gamma)))/(4*f*e*gamma);
00148 G4double xyP1K2 = 0;
00149 G4double xzP1K2 = (-2 + 2*e + gamma)/(2*sqrttwo*f2*gamma*
00150 std::sqrt(e*(1 + e + gamma - e*gamma)));
00151 G4double yxP1K2 = 0;
00152 G4double yyP1K2 = (1 - 2*gamma + e*(-1 + gamma*(2 + gamma)))/(4*f2*e*gamma2);
00153 G4double yzP1K2 = 0;
00154 G4double zxP1K2 = (2*e*(1 + e*gmo - 2*gamma) + gamma)/(2*sqrttwo*f2*gamma*
00155 std::sqrt(e*(1 + e + gamma - e*gamma)));
00156 G4double zyP1K2 = 0;
00157 G4double zzP1K2 = (1 - 2*gamma + e*(-1 - 2*e*gmo + 3*gamma))/
00158 (4*f2*e*gamma*std::sqrt(1 - (2*e)/(f*gpo)));
00159 phi2[0] += xxP1K2*pol1.x() + xyP1K2*pol1.y() + xzP1K2*pol1.z();
00160 phi2[1] += yxP1K2*pol1.x() + yyP1K2*pol1.y() + yzP1K2*pol1.z();
00161 phi2[2] += zxP1K2*pol1.x() + zyP1K2*pol1.y() + zzP1K2*pol1.z();
00162 }
00163
00164
00165
00166
00167
00168 if (!pol0.IsZero()) {
00169
00170
00171 G4double xxP2K1 = (-1 + e + e*gamma)/(4*f2*gamma*
00172 std::sqrt((e*(2 + e*gmo))/gpo));
00173 G4double xyP2K1 = 0;
00174 G4double xzP2K1 = -((1 + 2*f*gamma)*std::sqrt(f/(-2 + e - e*gamma)))/
00175 (2*sqrttwo*f2*e*gamma);
00176 G4double yxP2K1 = 0;
00177 G4double yyP2K1 = (1 - 2*gamma + e*(-1 + gamma*(2 + gamma)))/(4*f2*e*gamma2);
00178 G4double yzP2K1 = 0;
00179 G4double zxP2K1 = (1 + 2*e*(-2 + e + gamma - e*gamma))/(2*sqrttwo*f*e*
00180 std::sqrt(-(f*(2 + e*gmo)))*gamma);
00181 G4double zyP2K1 = 0;
00182 G4double zzP2K1 = (std::sqrt((e*gpo)/(2 + e*gmo))*
00183 (-3 + e*(5 + 2*e*gmo - 3*gamma) + 2*gamma))/(4*f2*e*gamma);
00184
00185 phi3[0] += xxP2K1*pol0.x() + xyP2K1*pol0.y() + xzP2K1*pol0.z();
00186 phi3[1] += yxP2K1*pol0.x() + yyP2K1*pol0.y() + yzP2K1*pol0.z();
00187 phi3[2] += zxP2K1*pol0.x() + zyP2K1*pol0.y() + zzP2K1*pol0.z();
00188 }
00189
00190 if (!pol1.IsZero()) {
00191
00192 G4double xxP2K2 = (-2 - e*(-3 + gamma) + gamma)/
00193 (4*f*e*gamma* std::sqrt((e*(2 + e*gmo))/gpo));
00194 G4double xyP2K2 = 0;
00195 G4double xzP2K2 = ((-2*e + gamma)*std::sqrt(f/(-2 + e - e*gamma)))/
00196 (2*sqrttwo*f*e2*gamma);
00197 G4double yxP2K2 = 0;
00198 G4double yyP2K2 = (-gamma2 + e*(-1 + gamma*(2 + gamma)))/(4*f*e2*gamma2);
00199 G4double yzP2K2 = 0;
00200 G4double zxP2K2 = (gamma + 2*e*(-1 + e - e*gamma))/
00201 (2*sqrttwo*e2* std::sqrt(-(f*(2 + e*gmo)))*gamma);
00202 G4double zyP2K2 = 0;
00203 G4double zzP2K2 = (std::sqrt((e*gpo)/(2 + e*gmo))*
00204 (-2 + e*(3 + 2*e*gmo - gamma) + gamma))/(4*f*e2*gamma);
00205 phi3[0] += xxP2K2*pol1.x() + xyP2K2*pol1.y() + xzP2K2*pol1.z();
00206 phi3[1] += yxP2K2*pol1.x() + yyP2K2*pol1.y() + yzP2K2*pol1.z();
00207 phi3[2] += zxP2K2*pol1.x() + zyP2K2*pol1.y() + zzP2K2*pol1.z();
00208 }
00209 }
00210 phi0 *= pref;
00211 phi2 *= pref;
00212 phi3 *= pref;
00213 }
00214
00215 G4double G4PolarizedMollerCrossSection::XSection(const G4StokesVector & pol2,
00216 const G4StokesVector & pol3)
00217 {
00218 G4double xs=0.;
00219 xs+=phi0;
00220
00221 G4bool polarized=(!pol2.IsZero())||(!pol3.IsZero());
00222 if (polarized) {
00223 xs+=phi2*pol2 + phi3*pol3;
00224 }
00225 return xs;
00226 }
00227
00228 G4double G4PolarizedMollerCrossSection::TotalXSection(
00229 G4double xmin, G4double xmax, G4double gamma,
00230 const G4StokesVector & pol0,const G4StokesVector & pol1)
00231 {
00232 G4double xs=0.;
00233
00234 G4double x=xmin;
00235
00236 if (xmax != 1./2.) G4cout<<" warning xmax expected to be 1/2 but is "<<xmax<< G4endl;
00237
00238
00239 G4double re2 = classic_electr_radius * classic_electr_radius;
00240 G4double gamma2=gamma*gamma;
00241 G4double gmo2 = (gamma - 1.)*(gamma - 1.);
00242 G4double logMEM = std::log(1./x - 1.);
00243 G4double pref = twopi*gamma2*re2/(gmo2*(gamma + 1.0));
00244
00245 G4double sigma0 = 0.;
00246 sigma0 += (gmo2/gamma2)*(0.5 - x);
00247 sigma0 += ((1. - 2.*gamma)/gamma2)*logMEM;
00248 sigma0 += 1./x - 1./(1. - x);
00249
00250 G4double sigma2=0.;
00251 sigma2 += ((gamma2 + 2.*gamma - 3.)/gamma2)*(0.5 - x);
00252 sigma2 += (1./gamma - 2.)*logMEM;
00253
00254 G4double sigma3=0.;
00255 sigma3 += (2.*(1. - gamma)/gamma2)*(0.5 - x);
00256 sigma3 += (1. - 3.*gamma)/(2.*gamma2)*logMEM;
00257
00258 xs+=pref*(sigma0 + sigma2*pol0.z()*pol1.z() + sigma3*(pol0.x()*pol1.x()+pol0.y()*pol1.y()));
00259
00260 return xs;
00261 }
00262
00263
00264 G4StokesVector G4PolarizedMollerCrossSection::GetPol2()
00265 {
00266
00267
00268 return 1./phi0 * phi2;
00269 }
00270 G4StokesVector G4PolarizedMollerCrossSection::GetPol3()
00271 {
00272
00273
00274 return 1./phi0 * phi3;
00275 }