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 "G4PolarizedAnnihilationCrossSection.hh"
00050 #include "G4PhysicalConstants.hh"
00051
00052 G4PolarizedAnnihilationCrossSection::G4PolarizedAnnihilationCrossSection() :
00053 polxx(0.), polyy(0.), polzz(0.), polxz(0.), polzx(0.), polxy(0.),
00054 polyx(0.), polyz(0.), polzy(0.),
00055 re2(1.), diffXSFactor(1.), totalXSFactor(1.),
00056 phi0(0.)
00057 {
00058 re2 = classic_electr_radius * classic_electr_radius;
00059 phi2 = G4ThreeVector(0., 0., 0.);
00060 phi3 = G4ThreeVector(0., 0., 0.);
00061 dice = 0.;
00062 polXS= 0.;
00063 unpXS = 0.;
00064 ISPxx=ISPyy=ISPzz=ISPnd=0.;
00065 }
00066
00067
00068
00069 G4PolarizedAnnihilationCrossSection::~G4PolarizedAnnihilationCrossSection()
00070 {
00071 }
00072
00073
00074
00075 void G4PolarizedAnnihilationCrossSection::TotalXS()
00076 {
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 }
00087
00088
00089
00090 void G4PolarizedAnnihilationCrossSection::Initialize(
00091 G4double eps,
00092 G4double gam,
00093 G4double ,
00094 const G4StokesVector & pol0,
00095 const G4StokesVector & pol1,
00096 G4int flag)
00097 {
00098
00099 diffXSFactor=re2/(gam - 1.);
00100 DefineCoefficients(pol0,pol1);
00101
00102
00103
00104 dice = 0.;
00105 G4double symmXS = 0.125*((-1./sqr(gam + 1.))/sqr(eps) +
00106 ((sqr(gam) + 4.*gam - 1.)/sqr(gam + 1.))/eps - 1.);
00107
00108
00109
00110 G4ThreeVector epsVector(1./sqr(eps), 1./eps, 1.);
00111 G4ThreeVector oneEpsVector(1./sqr(1. - eps), 1./(1.-eps), 1.);
00112 G4ThreeVector sumEpsVector(epsVector + oneEpsVector);
00113 G4ThreeVector difEpsVector(epsVector - oneEpsVector);
00114 G4ThreeVector calcVector(0., 0., 0.);
00115
00116
00117
00118 G4double helpVar2 = 0., helpVar1 = 0.;
00119
00120
00121
00122 helpVar1 = (gam*gam + 4.*gam + 1.)/sqr(gam + 1.);
00123 helpVar2 = -1./sqr(gam + 1.);
00124 calcVector = G4ThreeVector(helpVar2, helpVar1, -1.);
00125 unpXS = 0.125 * calcVector * sumEpsVector;
00126
00127
00128 helpVar2 = 1./sqr(gam + 1.);
00129 helpVar1 = -(gam*gam + 4.*gam + 1.)/sqr(gam + 1.);
00130 calcVector = G4ThreeVector(helpVar2, helpVar1, 0.5*(gam + 3.));
00131 ISPxx = 0.25*(calcVector * sumEpsVector)/(gam - 1.);
00132
00133 helpVar1 = 1./sqr(gam + 1.);
00134 calcVector = G4ThreeVector(-helpVar1, 2.*gam*helpVar1, -1.);
00135 ISPyy = 0.125 * calcVector * sumEpsVector;
00136
00137 helpVar1 = 1./(gam - 1.);
00138 helpVar2 = 1./sqr(gam + 1.);
00139 calcVector = G4ThreeVector(-(gam*gam + 1.)*helpVar2,(gam*gam*(gam + 1.) + 7.*gam + 3.)*helpVar2, -(gam + 3.));
00140 ISPzz = 0.125*helpVar1*(calcVector * sumEpsVector);
00141
00142 helpVar1 = std::sqrt(std::fabs(eps*(1. - eps)*2.*(gam + 1.) - 1.));
00143 calcVector = G4ThreeVector(-1./(gam*gam - 1.), 2./(gam - 1.), 0.);
00144 ISPnd = 0.125*(calcVector * difEpsVector) * helpVar1;
00145
00146 polXS = 0.;
00147 polXS += ISPxx*polxx;
00148 polXS += ISPyy*polyy;
00149 polXS += ISPzz*polzz;
00150 polXS += ISPnd*(polzx + polxz);
00151 phi0 = unpXS + polXS;
00152 dice = symmXS;
00153
00154 if(polzz != 0.) {
00155 dice *= (1. + (polzz*ISPzz/unpXS));
00156 if (dice<0.) dice=0.;
00157 }
00158
00159 if (flag==2) {
00160
00161
00162
00163 G4double circ1 = 0., circ2 = 0., circ3 = 0.;
00164 helpVar1 = 8.*sqr(1. - eps)*sqr(eps)*(gam - 1.)*sqr(gam + 1.)/std::sqrt(gam*gam - 1.);
00165 helpVar2 = sqr(gam + 1.)*sqr(eps)*(-2.*eps + 3.) - (gam*gam + 3.*gam + 2.)*eps;
00166 circ1 = helpVar2 + gam;
00167 circ1 /= helpVar1;
00168 circ2 = helpVar2 + 1.;
00169 circ2 /= helpVar1;
00170 helpVar1 = std::sqrt(std::fabs(eps*(1. - eps)*2.*(gam + 1.) - 1.));
00171 helpVar1 /= std::sqrt(gam*gam - 1.);
00172 calcVector = G4ThreeVector(1., -2.*gam, 0.);
00173 circ3 = 0.125*(calcVector * sumEpsVector)/(gam + 1.);
00174 circ3 *= helpVar1;
00175
00176 phi2.setZ( circ2*pol1.z() + circ1*pol0.z() + circ3*(pol1.x() + pol0.x()));
00177 phi3.setZ(-circ1*pol1.z() - circ2*pol0.z() - circ3*(pol1.x() + pol0.x()));
00178
00179
00180
00181 calcVector = G4ThreeVector(-1., 2.*gam, 0.);
00182 G4double linearZero = 0.125*(calcVector * sumEpsVector)/sqr(gam + 1.);
00183
00184
00185
00186 helpVar1 = std::sqrt(std::fabs(2.*(gam + 1.)*(1. - eps)*eps - 1.))/((gam + 1.)*eps*(1. - eps));
00187 helpVar2 = helpVar1*helpVar1;
00188
00189
00190
00191 G4double diagContrib = 0.125*helpVar2*(polxx + polyy - polzz);
00192 G4double nonDiagContrib = 0.125*helpVar1*(-polxz/(1. - eps) + polzx/eps);
00193
00194 phi2.setX(linearZero + diagContrib + nonDiagContrib);
00195
00196
00197
00198 nonDiagContrib = 0.125*helpVar1*(polxz/eps - polzx/(1. - eps));
00199
00200
00201 phi3.setX(linearZero + diagContrib + nonDiagContrib);
00202
00203
00204
00205 helpVar1 = std::sqrt(gam*gam - 1.)*(2.*(gam + 1.)*eps*(1. - eps) - 1.);
00206 helpVar1 /= 8.*sqr(1. - eps)*sqr(eps)*sqr(gam + 1.)*(gam - 1.);
00207 helpVar2 = std::sqrt((gam*gam - 1.)*std::fabs(2.*(gam + 1.)*eps*(1. - eps) - 1.));
00208 helpVar2 /= 8.*sqr(1. - eps)*sqr(eps)*sqr(gam + 1.)*(gam - 1.);
00209
00210 G4double contrib21 = (-polxy + polyx)*helpVar1;
00211 G4double contrib32 = -(eps*(gam + 1.) - 1.)*polyz + (eps*(gam + 1.) - gam)*polzy;
00212
00213 contrib32 *=helpVar2;
00214 phi2.setY(contrib21 + contrib32);
00215
00216 contrib32 = -(eps*(gam + 1.) - gam)*polyz + (eps*(gam + 1.) - 1.)*polzy;
00217 contrib32 *=helpVar2;
00218 phi3.setY(contrib21 + contrib32);
00219
00220 }
00221 phi0 *= diffXSFactor;
00222 phi2 *= diffXSFactor;
00223 phi3 *= diffXSFactor;
00224 }
00225
00226
00227
00228
00229 G4double G4PolarizedAnnihilationCrossSection::XSection(const G4StokesVector & pol2,
00230 const G4StokesVector & pol3)
00231 {
00232 G4double xs=phi0+pol2*phi2+pol3*phi3;
00233 return xs;
00234 }
00235
00236
00237
00238 G4double G4PolarizedAnnihilationCrossSection::TotalXSection(
00239 G4double ,G4double ,G4double gam,
00240 const G4StokesVector & pol0,const G4StokesVector & pol1)
00241 {
00242 totalXSFactor =pi*re2/(gam + 1.);
00243 DefineCoefficients(pol0,pol1);
00244
00245 G4double xs = 0.;
00246
00247
00248 G4double gam2 = gam*gam;
00249 G4double sqrtgam1 = std::sqrt(gam2 - 1.);
00250 G4double logMEM = std::log(gam+sqrtgam1);
00251 G4double unpME = (gam*(gam + 4.) + 1.)*logMEM;
00252 unpME += -(gam + 3.)*sqrtgam1;
00253 unpME /= 4.*(gam2 - 1.);
00254
00255
00256
00257 G4double longPart = (3+gam*(gam*(gam + 1.) + 7.))*logMEM;
00258 longPart += - (5.+ gam*(3*gam + 4.))*sqrtgam1;
00259 longPart /= 4.*sqr(gam - 1.)*(gam + 1.);
00260 G4double tranPart = -(5*gam + 1.)*logMEM;
00261 tranPart += (gam + 5.)*sqrtgam1;
00262 tranPart /= 4.*sqr(gam - 1.)*(gam + 1.);
00263
00264 xs += unpME;
00265 xs += polzz*longPart;
00266 xs += (polxx + polyy)*tranPart;
00267
00268 return xs*totalXSFactor;
00269 }
00270
00271
00272 G4StokesVector G4PolarizedAnnihilationCrossSection::GetPol2()
00273 {
00274
00275
00276 return 1./phi0 * phi2;
00277 }
00278
00279
00280
00281 G4StokesVector G4PolarizedAnnihilationCrossSection::GetPol3()
00282 {
00283
00284
00285 return 1./phi0 * phi3;
00286 }
00287
00288 void G4PolarizedAnnihilationCrossSection::DefineCoefficients(const G4StokesVector & pol0,
00289 const G4StokesVector & pol1)
00290 {
00291 polxx=pol0.x()*pol1.x();
00292 polyy=pol0.y()*pol1.y();
00293 polzz=pol0.z()*pol1.z();
00294
00295 polxz=pol0.x()*pol1.z();
00296 polzx=pol0.z()*pol1.x();
00297
00298 polyz=pol0.y()*pol1.z();
00299 polzy=pol0.z()*pol1.y();
00300
00301 polxy=pol0.x()*pol1.y();
00302 polyx=pol0.y()*pol1.x();
00303 }
00304
00305
00306
00307 G4double G4PolarizedAnnihilationCrossSection::GetXmin(G4double y)
00308 {
00309 return 0.5*(1.-std::sqrt((y-1.)/(y+1.)));
00310 }
00311 G4double G4PolarizedAnnihilationCrossSection::GetXmax(G4double y)
00312 {
00313 return 0.5*(1.+std::sqrt((y-1.)/(y+1.)));
00314 }
00315
00316
00317 G4double G4PolarizedAnnihilationCrossSection::DiceEpsilon()
00318 {
00319 return dice;
00320 }
00321
00322 G4double G4PolarizedAnnihilationCrossSection::getVar(G4int choice)
00323 {
00324 if (choice == -1) return polXS/unpXS;
00325 if (choice == 0) return unpXS;
00326 if (choice == 1) return ISPxx;
00327 if (choice == 2) return ISPyy;
00328 if (choice == 3) return ISPzz;
00329 if (choice == 4) return ISPnd;
00330 return 0;
00331 }
00332