00001 #include "TMuKFun.h"
00002 #include "TLi2.h"
00003 #include "TRadGlobal.h"
00004
00005 #include <fstream>
00006 #include <string>
00007 #include <sstream>
00008 #include <cstdlib>
00009
00010 double TMuKFun::Eval_even(const double &s, const double &c){
00011 double r = gGlobal->Get_MF2()*0.25*gGlobal->Get_s()/s;
00012 double rho =-log(r);
00013 double beta = sqrt(1-4*r);
00014 double beta1p = 1+beta;
00015 double beta1m = 1-beta;
00016 double rbeta = beta1m/beta1p;
00017 double ibeta = 1/beta;
00018 double beta2 = beta*beta;
00019 double lbeta = -log(rbeta);
00020 double k_even = -1 + rho*(1.5 + beta2 - beta)*0.5*ibeta +
00021 log(0.5*beta1p)*(1.5+beta2)*ibeta - 0.5*ibeta*(1-beta2)*lbeta/(2-beta2*(1-c*c))+
00022 0.5*ibeta*(1+beta2)*(gConst->Pi2()/6 + 2*TLi2::Eval(rbeta) + lbeta*log(0.5*ibeta*ibeta*beta1p));
00023
00024
00025 double iv = 0.5*ibeta*(1+beta2);
00026
00027 return k_even - 0.5*gConst->Pi2()*iv;
00028 }
00029
00030 double TMuKFun::Eval_Coloumb_Factor(const double &s){
00031
00032
00033
00034
00035 double beta = sqrt(1 - gGlobal->Get_MF2()*gGlobal->Get_s()/s);
00036
00037 double iv = (1+beta*beta)/(2*beta);
00038 double z = 2*gConst->Pi()*gConst->Alpha()*iv;
00039
00040 return z/(1. - exp(-z));
00041 }
00042
00043 double TMuKFun::Eval_even_ultra(const double &s){
00044 return 3./4.*log(s/(gGlobal->Get_MF2()*0.25*gGlobal->Get_s())) + gConst->Pi2()/6. - 1.;
00045 }
00046
00047 double TMuKFun::Eval_odd(const double &s, const double& c){
00048 double r = gGlobal->Get_MF2()*0.25*gGlobal->Get_s()/s;
00049 double rho =-log(r);
00050 double beta2 = 1 - 4*r;
00051 double beta = sqrt(beta2);
00052 double ibeta = 1/beta;
00053 double ibeta2 = 0.5*ibeta*ibeta;
00054
00055 double d = -rho*ibeta2 + (gConst->Pi2()/12 + 0.25*rho*rho)*
00056 (1 - ibeta - 0.5*beta + ibeta2*ibeta) + ibeta*(-1 - 0.5*beta2 + ibeta2)*
00057 (rho*log(0.5*(1 + beta)) - 2*TLi2::Eval(0.5*(1-beta)) - TLi2::Eval(-(1-beta)/(1+beta)));
00058
00059 double betac = beta*c;
00060 double betas2 = beta2*(1-c*c);
00061 double i2mbetas2 = 1/(2-betas2);
00062 double tY = 1 - beta2 - betas2;
00063
00064 double l_m = log(0.5*(1-betac));
00065 double b = rho + l_m;
00066 double t1 = (1-beta2)/(2*(1-betac));
00067 double a = 0.5*l_m*l_m - log(1 - t1)*b + TLi2::Eval(t1);
00068 double inv1 = 1/(1 + beta2 - 2*betac);
00069 double li2t2 = TLi2::Eval(betas2*inv1);
00070
00071 double res_plus = a + li2t2 + i2mbetas2*(-tY*b*inv1 - 2*r*a + betac*(d - a));
00072
00073
00074 betac = -betac;
00075 l_m = log(0.5*(1-betac));
00076 b = rho + l_m;
00077 t1 = (1-beta2)/(2*(1-betac));
00078 a = 0.5*l_m*l_m - log(1 - t1)*b + TLi2::Eval(t1);
00079 inv1 = 1/(1 + beta2 - 2*betac);
00080 li2t2 = TLi2::Eval(betas2*inv1);
00081
00082 double res_minus= a + li2t2 + i2mbetas2*(-tY*b*inv1 - 2*r*a + betac*(d - a));
00083
00084 return res_plus - res_minus + Eval_int_odd(beta,c);
00085 }
00086
00087 double TMuKFun::Eval_odd_ultra(const double& c){
00088
00089 double st_h2 = 0.5*(1-c);
00090 double st_h = sqrt(st_h2);
00091 double ct_h2 = 0.5*(1+c);
00092 double ct_h = sqrt(ct_h2);
00093 double lst_h = log(st_h);
00094 double lct_h = log(ct_h);
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 return 2*lst_h*lst_h - 2*lct_h*lct_h - TLi2::Eval(st_h2) + TLi2::Eval(ct_h2)
00106 + 2/(1+c*c)*(ct_h2*lst_h - st_h2*lct_h - c*(lst_h*lst_h +lct_h*lct_h ));
00107 }
00108
00109 double TMuKFun::Eval_even_odd(const double &s, const double &c){
00110
00111
00112 double r = gGlobal->Get_MF2()*0.25*gGlobal->Get_s()/s;
00113 double rho =-log(r);
00114 double beta2 = 1 - 4*r;
00115 double beta = sqrt(beta2);
00116 double beta1p = 1+beta;
00117 double beta1m = 1-beta;
00118 double rbeta = beta1m/beta1p;
00119 double ibeta = 1/beta;
00120 double lbeta = -log(rbeta);
00121 double betas2 = beta2*(1-c*c);
00122 double i2mbetas2 = 1/(2-betas2);
00123 double k_even = -1 + rho*(1.5 + beta2 - beta)*0.5*ibeta +
00124 log(0.5*beta1p)*(1.5+beta2)*ibeta - 0.5*ibeta*(1-beta2)*lbeta*i2mbetas2+
00125 0.5*ibeta*(1+beta2)*(gConst->Pi2()/6 + 2*TLi2::Eval(rbeta) + lbeta*log(0.5*ibeta*ibeta*beta1p));
00126
00127
00128 double iv = 0.5*ibeta*(1+beta2);
00129 k_even -= 0.5*gConst->Pi2()*iv;
00130
00131 double ibeta2 = 0.5*ibeta*ibeta;
00132 double d = -rho*ibeta2 + (gConst->Pi2()/12 + 0.25*rho*rho)*
00133 (1 - ibeta - 0.5*beta + ibeta2*ibeta) + ibeta*(-1 - 0.5*beta2 + ibeta2)*
00134 (rho*log(0.5*beta1p) - 2*TLi2::Eval(0.5*beta1m) - TLi2::Eval(-rbeta));
00135
00136 double betac = beta*c;
00137 double tY = 1 - beta2 - betas2;
00138
00139 double l_m = log(0.5*(1-betac));
00140 double b = rho + l_m;
00141 double t1 = (1-beta2)/(2*(1-betac));
00142 double a = 0.5*l_m*l_m - log(1 - t1)*b + TLi2::Eval(t1);
00143 double inv1 = 1/(1 + beta2 - 2*betac);
00144 double li2t2 = TLi2::Eval(betas2*inv1);
00145
00146 double res_plus = a + li2t2 + i2mbetas2*(-tY*b*inv1 - 2*r*a + betac*(d - a));
00147
00148
00149 betac = -betac;
00150 l_m = log(0.5*(1-betac));
00151 b = rho + l_m;
00152 t1 = (1-beta2)/(2*(1-betac));
00153 a = 0.5*l_m*l_m - log(1 - t1)*b + TLi2::Eval(t1);
00154 inv1 = 1/(1 + beta2 - 2*betac);
00155 li2t2 = TLi2::Eval(betas2*inv1);
00156
00157 double res_minus= a + li2t2 + i2mbetas2*(-tY*b*inv1 - 2*r*a + betac*(d - a));
00158
00159 return k_even + res_plus - res_minus + Eval_int_odd(beta,c);
00160 }
00161
00162 #define NX 50
00163 #define NY 101
00164 #define T(i,j) fm[NY*(i)+j]
00165
00166 void TMuKFun::Init(std::string intd, std::string intf){
00167 if(!fm){
00168 fm = new float[NX*NY];
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185 std::string fname;
00186 if(! intf.empty()) fname = intf;
00187 else fname = "integ.dat";
00188
00189 if(! intd.empty()) fname = intd + "/" + fname;
00190
00191 std::ifstream IN(fname.c_str());
00192 if(IN.is_open() == false ){
00193 std::cout<<"Can't read \""<<fname<<"\". "
00194 <<"Angular distribution will be partially incorrect!"
00195 <<std::endl;
00196 std::cout<<"Probably you should tune INTEG_DIR or INTEG_FNAME variables."<<std::endl;
00197 for(int i=0;i<NX;i++){
00198 for(int j=0;j<NY;j++)
00199 T(i,j) = 0;
00200 }
00201 return;
00202 }
00203 float t1,t2;
00204 for(int i=0;i<NX;i++){
00205 for(int j=0;j<NY-1;j++)
00206 IN>>t1>>t2>>T(i,j);
00207 T(i,0) = 0;
00208 T(i,NY-1) = 0;
00209 }
00210 IN.close();
00211 }
00212 }
00213
00214 double TMuKFun::Eval_int_odd(const double &beta, const double &c){
00215 int i = int(floor((beta - 0.01)*50));
00216 int j = int(floor((c + 1)*50));
00217 double wy = 50*c - (j - 50);
00218 double res = 0;
00219 if(i >= 0 && i < NX-1){
00220 double wx = 50*beta - i - 0.5;
00221 res =
00222 (T(i,j )*(1-wx) + T(i+1,j )*wx)*(1-wy) +
00223 (T(i,j+1)*(1-wx) + T(i+1,j+1)*wx)*wy;
00224 } else if (i == -1){
00225 double wx = 100*beta;
00226 res =
00227 (T(0,j )*(1-wy) + T(0, j+1)*wy)*wx;
00228 } else if (i == NX - 1){
00229 double wx = 1 - 100*(beta-0.99);
00230 res =
00231 (T(NX-1,j)*(1-wy)+ T(NX-1,j+1)*wy)*wx;
00232 }
00233 return -res;
00234 }
00235
00236 TMuKFun::TMuKFun(){
00237 fm = NULL;
00238
00239 }
00240
00241 TMuKFun::~TMuKFun(){
00242 delete fm;
00243 }