TEvtGen/EvtGenModels/EvtWilsonCoeficients.cxx

   1 //---------------------------------------------------------------------------------
   2 //
   3 // Wilson coeficients according to A.J.Buras and M.Munz, Phys.Rev. D52, 186. (1995)
   4 // Thanks to N. Nikitine for example code for Pythia
   5 // Coeficient C8eff and C2 correction to C7eff taken from:
   6 //   A.J.Buras, M.Misiak, M.Munz, S.Pokorski, Nucl.Phys. B424, 374 (1994)
   7 //
   8 // Used constants come from PDG 2004
   9 //
  10 // P. Reznicek 18.02.2005
  11 //
  12 //  04/03/2005  PR   Added h-function
  13 //
  14 //---------------------------------------------------------------------------------
  15
  16 #include "EvtGenBase/EvtPatches.hh"
  17 #include "EvtGenBase/EvtConst.hh"
  18 #include "EvtGenBase/EvtReport.hh"
  19 #include "EvtGenModels/EvtWilsonCoeficients.hh"
  20 #include <stdlib.h>
  21
  22 // including EvtLi2Spence.F
  23 extern "C" {
  24   extern double li2spence_(double*);
  25 }
  26
  27 EvtWilsonCoeficients::EvtWilsonCoeficients(){
  28   int i,j;
  29   double tmpa[8]={14./23.,16./23.,6./23.,-12./23.,0.4086,-0.4230,-0.8994,0.1456};
  30   double tmph[8]={2.2996,-1.0880,-3./7.,-1./14.,-0.6494,-0.0380,-0.0186,-0.0057};
  31   double tmpp[8]={0,0,-80./203.,8./33.,0.0433,0.1384,0.1648,-0.0073};
  32   double tmps[8]={0,0,-0.2009,-0.3579,0.0490,-0.3616,-0.3554,0.0072};
  33   double tmpq[8]={0,0,0,0,0.0318,0.0918,-0.2700,0.0059};
  34   double tmpg[8]={313063./363036.,0,0,0,-0.9135,0.0873,-0.0571,-0.0209};
  35   double tmpk[6][8]={{0,0,1./2.,-1./2.,0,0,0,0},
  36                      {0,0,1./2.,+1./2.,0,0,0,0},
  37                      {0,0,-1./14.,+1./6.,0.0510,-0.1403,-0.0113,0.0054},
  38                      {0,0,-1./14.,-1./6.,0.0984,+0.1214,+0.0156,0.0026},
  39                      {0,0,0,0,-0.0397,0.0117,-0.0025,+0.0304},
  40                      {0,0,0,0,+0.0335,0.0239,-0.0462,-0.0112}};
  41   double tmpr[2][8]={{0,0,+0.8966,-0.1960,-0.2011,0.1328,-0.0292,-0.1858},
  42                      {0,0,-0.1193,+0.1003,-0.0473,0.2323,-0.0133,-0.1799}};
  43   for(i=0;i<8;i++){
  44     a[i]=tmpa[i];
  45     h[i]=tmph[i];
  46     p[i]=tmpp[i];
  47     s[i]=tmps[i];
  48     q[i]=tmpq[i];
  49     g[i]=tmpg[i];
  50     for(j=0;j<6;j++) k[j][i]=tmpk[j][i];
  51     for(j=0;j<2;j++) r[j][i]=tmpr[j][i];
  52   }
  53   m_n_f=5;
  54   m_Lambda=0.2167;
  55   m_alphaMZ=0.1187;
  56   m_mu=4.8;
  57   m_M_Z=91.1876;
  58   m_M_t=174.3;
  59   m_M_W=80.425;
  60   m_alphaS=0;
  61   m_eta=0;
  62   m_sin2W=0.23120;
  63   m_ialpha=137.036;
  64   m_C1=m_C2=m_C3=m_C4=m_C5=m_C6=m_C7=m_C7eff0=m_C8=m_C8eff0=m_C9=m_C9tilda=m_C10=m_C10tilda=m_P0=0;
  65   m_A=m_B=m_C=m_D=m_E=m_F=m_Y=m_Z=m_PE=0;
  66   m_ksi=0;
  67 }
  68
  69 void EvtWilsonCoeficients::SetRenormalizationScheme(std::string scheme){
  70   if(scheme=="NDR") m_ksi=0;
  71   else if(scheme=="HV") m_ksi=1;
  72   else{
  73     report(ERROR,"EvtGen") << "ERROR: EvtWilsonCoeficients knows only NDR and HV schemes !" << std::endl;
  74     ::abort();
  75   }
  76 }
  77
  78 double EvtWilsonCoeficients::alphaS(double mu=4.8,int n_f=5,double Lambda=0.2167){
  79 // calculate strong coupling constant for n_f flavours and scale mu
  80   double beta0=11.-2./3.*n_f;
  81   double beta1=51.-19./3.*n_f;
  82   double beta2=2857.-5033./9.*n_f+325./27.*n_f*n_f;
  83   double lnratio=log(mu*mu/Lambda/Lambda);
  84   double aS=4.*EvtConst::pi/beta0/lnratio*(1.-2*beta1/beta0/beta0*log(lnratio)/lnratio+
  85             4*beta1*beta1/beta0/beta0/beta0/beta0/lnratio/lnratio*((log(lnratio)-0.5)*(log(lnratio)-0.5)+beta2*beta0/8/beta1/beta1-5./4.));
  86   return aS;
  87 }
  88
  89 double EvtWilsonCoeficients::Lambda(double alpha=0.1187,int n_f=5,double mu=91.1876,double epsilon=0.00005,int maxstep=1000){
  90 // calculate Lambda matching alphaS using simple iterative method
  91   int i;
  92   double difference=0;
  93   double Lambda=mu*0.9999999999;
  94   double step=-mu/20;
  95   for(i=0;i<maxstep && (difference=fabs(alphaS(mu,n_f,Lambda)-alpha))>=epsilon;i++){
  96     report(DEBUG,"EvtGen") << " Difference of alpha_S from " << alpha << " is " << difference << " at Lambda = " << Lambda << std::endl;
  97     if(alphaS(mu,n_f,Lambda)>alpha){
  98       if(step>0) step*=-0.4;
  99       if(alphaS(mu,n_f,Lambda+step-epsilon)<alphaS(mu,n_f,Lambda+step)) Lambda+=step;
 100       else step*=0.4;
 101     }else{
 102       if(step<0) step*=-0.4;
 103       if(Lambda+step<mu) Lambda+=step;
 104       else step*=0.4;
 105     }
 106   }
 107   report(DEBUG,"EvtGen") << " Difference of alpha_S from " << alpha << " is " << difference << " at Lambda = " << Lambda << std::endl;
 108   if(difference>=epsilon){
 109     report(ERROR,"EvtGen") << " ERROR: Did not converge Lambda for alpha_s = " << alpha << " , difference " << difference << " >= " << epsilon << " after " << i << " steps !" << std::endl;
 110     ::abort();
 111     return -1;
 112   }else{
 113     report(INFO,"EvtGen") << " For alpha_s = " << alphaS(mu,n_f,Lambda) << " was found Lambda = " << Lambda << std::endl;
 114     return Lambda;
 115   }
 116 }
 117
 118 double EvtWilsonCoeficients::eta(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_W=80.425){
 119   return alphaS(M_W,n_f,Lambda)/alphaS(mu,n_f,Lambda);
 120 }
 121
 122
 123 EvtComplex EvtWilsonCoeficients::C1(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_W=80.425){
 124   int i;
 125   EvtComplex myC1(0,0);
 126   for(i=0;i<8;i++) myC1+=k[0][i]*pow(eta(mu,n_f,Lambda,M_W),a[i]);
 127   return myC1;
 128 }
 129
 130 EvtComplex EvtWilsonCoeficients::C2(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_W=80.425){
 131   int i;
 132   EvtComplex myC2(0,0);
 133   for(i=0;i<8;i++) myC2+=k[1][i]*pow(eta(mu,n_f,Lambda,M_W),a[i]);
 134   return myC2;
 135 }
 136
 137 EvtComplex EvtWilsonCoeficients::C3(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_W=80.425){
 138   int i;
 139   EvtComplex myC3(0,0);
 140   for(i=0;i<8;i++) myC3+=k[2][i]*pow(eta(mu,n_f,Lambda,M_W),a[i]);
 141   return myC3;
 142 }
 143
 144 EvtComplex EvtWilsonCoeficients::C4(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_W=80.425){
 145   int i;
 146   EvtComplex myC4(0,0);
 147   for(i=0;i<8;i++) myC4+=k[3][i]*pow(eta(mu,n_f,Lambda,M_W),a[i]);
 148   return myC4;
 149 }
 150
 151 EvtComplex EvtWilsonCoeficients::C5(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_W=80.425){
 152   int i;
 153   EvtComplex myC5(0,0);
 154   for(i=0;i<8;i++) myC5+=k[4][i]*pow(eta(mu,n_f,Lambda,M_W),a[i]);
 155   return myC5;
 156 }
 157
 158 EvtComplex EvtWilsonCoeficients::C6(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_W=80.425){
 159   int i;
 160   EvtComplex myC6(0,0);
 161   for(i=0;i<8;i++) myC6+=k[5][i]*pow(eta(mu,n_f,Lambda,M_W),a[i]);
 162   return myC6;
 163 }
 164
 165
 166 EvtComplex EvtWilsonCoeficients::C7(double M_t=174.3,double M_W=80.425){
 167   return EvtComplex(-0.5*A(M_t*M_t/M_W/M_W),0);
 168 }
 169
 170 EvtComplex EvtWilsonCoeficients::C8(double M_t=174.3,double M_W=80.425){
 171   return EvtComplex(-0.5*F(M_t*M_t/M_W/M_W),0);
 172 }
 173
 174 EvtComplex EvtWilsonCoeficients::C7eff0(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_t=174.3,double M_W=80.425){
 175   int i;
 176   EvtComplex myC7eff(0,0);
 177   for(i=0;i<8;i++) myC7eff+=h[i]*pow(eta(mu,n_f,Lambda,M_W),a[i]);
 178   myC7eff*=C2(mu,n_f,Lambda,M_W);
 179   myC7eff+=pow(eta(mu,n_f,Lambda,M_W),16./23.)*C7(M_t,M_W);
 180   myC7eff+=8./3.*(pow(eta(mu,n_f,Lambda,M_W),14./23.)-pow(eta(mu,n_f,Lambda,M_W),16./23.))*C8(M_t,M_W);
 181   return myC7eff;
 182 }
 183
 184 EvtComplex EvtWilsonCoeficients::C8eff0(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_t=174.3,double M_W=80.425){
 185   int i;
 186   EvtComplex myC8eff(0,0);
 187   for(i=0;i<8;i++) myC8eff+=g[i]*pow(eta(mu,n_f,Lambda,M_W),a[i]);
 188   myC8eff+=pow(eta(mu,n_f,Lambda,M_W),14./23.)*C8(M_t,M_W);
 189   return myC8eff;
 190 }
 191
 192
 193 EvtComplex EvtWilsonCoeficients::C10tilda(double sin2W=0.23120,double M_t=174.3,double M_W=80.425){
 194   return EvtComplex(-Y(M_t*M_t/M_W/M_W)/sin2W,0);
 195 }
 196
 197 EvtComplex EvtWilsonCoeficients::C10(double sin2W=0.23120,double M_t=174.3,double M_W=80.425,double ialpha=137.036){
 198   return ( 1./2/EvtConst::pi/ialpha*C10tilda(sin2W,M_t,M_W) );
 199 }
 200
 201
 202 double EvtWilsonCoeficients::A(double x){
 203   return ( x*(8*x*x+5*x-7)/12/pow(x-1,3) + x*x*(2-3*x)*log(x)/2/pow(x-1,4) );
 204 }
 205
 206 double EvtWilsonCoeficients::B(double x){
 207   return ( x/4/(1-x) + x/4/(x-1)/(x-1)*log(x) );
 208 }
 209
 210 double EvtWilsonCoeficients::C(double x){
 211   return ( x*(x-6)/8/(x-1) + x*(3*x+2)/8/(x-1)/(x-1)*log(x) );
 212 }
 213
 214 double EvtWilsonCoeficients::D(double x){
 215   return ( (-19*x*x*x+25*x*x)/36/pow(x-1,3) + x*x*(5*x*x-2*x-6)/18/pow(x-1,4)*log(x) - 4./9*log(x) );
 216 }
 217
 218 double EvtWilsonCoeficients::E(double x){
 219   return ( x*(18-11*x-x*x)/12/pow(1-x,3) + x*x*(15-16*x+4*x*x)/6/pow(1-x,4)*log(x) - 2./3*log(x) );
 220 }
 221
 222 double EvtWilsonCoeficients::F(double x){
 223   return ( x*(x*x-5*x-2)/4/pow(x-1,3) + 3*x*x/2/pow(x-1,4)*log(x) );
 224 }
 225
 226 double EvtWilsonCoeficients::Y(double x){
 227   return (C(x)-B(x));
 228 }
 229
 230 double EvtWilsonCoeficients::Z(double x){
 231   return (C(x)+1./4*D(x));
 232 }
 233
 234
 235 EvtComplex EvtWilsonCoeficients::C9(int ksi=0,double mu=4.8,int n_f=5,double Lambda=0.2167,double sin2W=0.23120,double M_t=174.3,double M_W=80.425,double ialpha=137.036){
 236   return ( 1./2/EvtConst::pi/ialpha*C9tilda(ksi,mu,n_f,Lambda,sin2W,M_t,M_W) );
 237 }
 238
 239 EvtComplex EvtWilsonCoeficients::C9tilda(int ksi=0,double mu=4.8,int n_f=5,double Lambda=0.2167,double sin2W=0.23120,double M_t=174.3,double M_W=80.425){
 240   return ( P0(ksi,mu,n_f,Lambda,M_W) + Y(M_t*M_t/M_W/M_W)/sin2W - 4*Z(M_t*M_t/M_W/M_W) + PE(mu,n_f,Lambda,M_W)*E(M_t*M_t/M_W/M_W) );
 241 }
 242
 243 EvtComplex EvtWilsonCoeficients::P0(int ksi=0,double mu=4.8,int n_f=5,double Lambda=0.2167,double M_W=80.425){
 244   int i;
 245   EvtComplex myP0(0,0);
 246   for(i=0;i<8;i++) myP0+=p[i]*pow(eta(mu,n_f,Lambda,M_W),a[i]+1);
 247   myP0=EvtConst::pi/alphaS(M_W,n_f,Lambda)*(-0.1875+myP0);
 248   myP0+=1.2468-ksi*4./9.*(3*C1(mu,n_f,Lambda,M_W)+C2(mu,n_f,Lambda,M_W)-C3(mu,n_f,Lambda,M_W)-3*C4(mu,n_f,Lambda,M_W));
 249   for(i=0;i<8;i++) myP0+=pow(eta(mu,n_f,Lambda,M_W),a[i])*(r[ksi][i]+s[i]*eta(mu,n_f,Lambda,M_W));
 250   return myP0;
 251 }
 252
 253 double EvtWilsonCoeficients::PE(double mu=4.8,int n_f=5,double Lambda=0.2167,double M_W=80.425){
 254   int i;
 255   double myPE=0.1405;
 256   for(i=0;i<8;i++) myPE+=q[i]*pow(eta(mu,n_f,Lambda,M_W),a[i]+1);
 257   return myPE;
 258 }
 259
 260
 261 void EvtWilsonCoeficients::CalculateAllCoeficients(){
 262   m_Lambda=Lambda(m_alphaMZ,m_n_f,m_M_Z);
 263   m_C1=C1(m_mu,m_n_f,m_Lambda,m_M_W);
 264   m_C2=C2(m_mu,m_n_f,m_Lambda,m_M_W);
 265   m_C3=C3(m_mu,m_n_f,m_Lambda,m_M_W);
 266   m_C4=C4(m_mu,m_n_f,m_Lambda,m_M_W);
 267   m_C5=C5(m_mu,m_n_f,m_Lambda,m_M_W);
 268   m_C6=C6(m_mu,m_n_f,m_Lambda,m_M_W);
 269   m_C7=C7(m_M_t,m_M_W);
 270   m_C8=C8(m_M_t,m_M_W);
 271   m_C7eff0=C7eff0(m_mu,m_n_f,m_Lambda,m_M_t,m_M_W);
 272   m_C8eff0=C8eff0(m_mu,m_n_f,m_Lambda,m_M_t,m_M_W);
 273   m_C10tilda=C10tilda(m_sin2W,m_M_t,m_M_W);
 274   m_C10=C10(m_sin2W,m_M_t,m_M_W,m_ialpha);
 275   m_A=A(m_M_t*m_M_t/m_M_W/m_M_W);
 276   m_B=B(m_M_t*m_M_t/m_M_W/m_M_W);
 277   m_C=C(m_M_t*m_M_t/m_M_W/m_M_W);
 278   m_D=D(m_M_t*m_M_t/m_M_W/m_M_W);
 279   m_E=E(m_M_t*m_M_t/m_M_W/m_M_W);
 280   m_F=F(m_M_t*m_M_t/m_M_W/m_M_W);
 281   m_Y=Y(m_M_t*m_M_t/m_M_W/m_M_W);
 282   m_Z=Z(m_M_t*m_M_t/m_M_W/m_M_W);
 283   m_C9=C9(m_ksi,m_mu,m_n_f,m_Lambda,m_sin2W,m_M_t,m_M_W,m_ialpha);
 284   m_C9tilda=C9tilda(m_ksi,m_mu,m_n_f,m_Lambda,m_sin2W,m_M_t,m_M_W);
 285   m_P0=P0(m_ksi,m_mu,m_n_f,m_Lambda,m_M_W);
 286   m_PE=PE(m_mu,m_n_f,m_Lambda,m_M_W);
 287   m_alphaS=alphaS(m_mu,m_n_f,m_Lambda);
 288   m_eta=eta(m_mu,m_n_f,m_Lambda,m_M_W);
 289   report(INFO,"EvtGen") << " +---------------------------------------" << std::endl;
 290   report(INFO,"EvtGen") << " | Table of Wilson coeficients:" << std::endl;
 291   report(INFO,"EvtGen") << " +---------------------------------------" << std::endl;
 292   report(INFO,"EvtGen") << " | C1     =  " << m_C1 << std::endl;
 293   report(INFO,"EvtGen") << " | C2     =  " << m_C2 << std::endl;
 294   report(INFO,"EvtGen") << " | C3     =  " << m_C3 << std::endl;
 295   report(INFO,"EvtGen") << " | C4     =  " << m_C4 << std::endl;
 296   report(INFO,"EvtGen") << " | C5     =  " << m_C5 << std::endl;
 297   report(INFO,"EvtGen") << " | C6     =  " << m_C6 << std::endl;
 298   report(INFO,"EvtGen") << " | C7     =  " << m_C7 << std::endl;
 299   report(INFO,"EvtGen") << " | C7eff0 =  " << m_C7eff0 << std::endl;
 300   report(INFO,"EvtGen") << " | C8     =  " << m_C8 << std::endl;
 301   report(INFO,"EvtGen") << " | C8eff0 =  " << m_C8eff0 << std::endl;
 302   report(INFO,"EvtGen") << " | C9     =  " << m_C9 << std::endl;
 303   report(INFO,"EvtGen") << " | C10    =  " << m_C10 << std::endl;
 304   report(INFO,"EvtGen") << " +---------------------------------------" << std::endl;
 305   report(INFO,"EvtGen") << " | Other constants:" << std::endl;
 306   report(INFO,"EvtGen") << " +---------------------------------------" << std::endl;
 307   report(INFO,"EvtGen") << " | Scale = " << m_mu << " GeV" << std::endl;
 308   report(INFO,"EvtGen") << " | Number of effective flavors = " << m_n_f << std::endl;
 309   report(INFO,"EvtGen") << " | Corresponding to aS(M_Z)" << "=" << m_alphaMZ << " Lambda = " << m_Lambda << " GeV" << std::endl;
 310   report(INFO,"EvtGen") << " | Strong coupling constant = " << m_alphaS << std::endl;
 311   report(INFO,"EvtGen") << " | Electromagnetic constant = 1/" << m_ialpha << std::endl;
 312   report(INFO,"EvtGen") << " | Top mass = " << m_M_t << " GeV" << std::endl;
 313   report(INFO,"EvtGen") << " | W-boson mass = " << m_M_W << " GeV" << std::endl;
 314   report(INFO,"EvtGen") << " | Z-boson mass = " << m_M_Z << " GeV" << std::endl;
 315   report(INFO,"EvtGen") << " | Sinus squared of Weinberg angle = " << m_sin2W << std::endl;
 316   report(INFO,"EvtGen") << " +---------------------------------------" << std::endl;
 317   report(DEBUG,"EvtGen") << " | Intermediate functions:" << std::endl;
 318   report(DEBUG,"EvtGen") << " +---------------------------------------" << std::endl;
 319   report(DEBUG,"EvtGen") << " | A    = " << m_A << std::endl;
 320   report(DEBUG,"EvtGen") << " | B    = " << m_B << std::endl;
 321   report(DEBUG,"EvtGen") << " | C    = " << m_C << std::endl;
 322   report(DEBUG,"EvtGen") << " | D    = " << m_D << std::endl;
 323   report(DEBUG,"EvtGen") << " | E    = " << m_E << std::endl;
 324   report(DEBUG,"EvtGen") << " | F    = " << m_F << std::endl;
 325   report(DEBUG,"EvtGen") << " | Y    = " << m_Y << std::endl;
 326   report(DEBUG,"EvtGen") << " | Z    = " << m_Z << std::endl;
 327   report(DEBUG,"EvtGen") << " | eta  = " << m_eta << std::endl;
 328   report(DEBUG,"EvtGen") << " | C9~  = " << m_C9tilda << std::endl;
 329   report(DEBUG,"EvtGen") << " | C10~ = " << m_C10tilda << std::endl;
 330   report(DEBUG,"EvtGen") << " | P0   = " << m_P0 << std::endl;
 331   report(DEBUG,"EvtGen") << " | PE   = " << m_PE << std::endl;
 332   report(DEBUG,"EvtGen") << " +--------------------------------------" << std::endl;
 333 }
 334
 335 EvtComplex EvtWilsonCoeficients::hzs(double z,double shat,double mu=4.8,double M_b=4.8){
 336   EvtComplex i1(0,1);
 337   double x=4.*z*z/shat;
 338   if(x==0)     return (8./27. - 8./9.*log(M_b/mu) - 4./9.*log(shat) + 4./9.*i1*EvtConst::pi);
 339   else if(x>1) return (8./27. - 8./9.*log(M_b/mu) - 8./9.*log(z) + 4./9.*x - 2./9.*(2.+x)*sqrt(x-1.) * 2*atan(1./sqrt(x-1.)));
 340   else         return (8./27. - 8./9.*log(M_b/mu) - 8./9.*log(z) + 4./9.*x - 2./9.*(2.+x)*sqrt(1.-x) * (log(fabs(sqrt(1.-x)+1)/fabs(sqrt(1.-x)-1))-i1*EvtConst::pi));
 341 }
 342
 343 double EvtWilsonCoeficients::fz(double z){
 344   return (1. - 8.*z*z + 8.*pow(z,6.) - pow(z,8.) - 24.*pow(z,4.)*log(z));
 345 }
 346
 347 double EvtWilsonCoeficients::kappa(double z,double alpha_S){
 348   return (1. - 2.*alpha_S/3./EvtConst::pi*((EvtConst::pi*EvtConst::pi-31./4.)*(1.-z)*(1.-z) + 1.5) );
 349 }
 350
 351 double EvtWilsonCoeficients::etatilda(double shat,double alpha_S){
 352   return (1. + alpha_S/EvtConst::pi*omega(shat));
 353 }
 354
 355 double EvtWilsonCoeficients::omega(double shat){
 356   double o=0;
 357   o -= (2./9.)*EvtConst::pi*EvtConst::pi;
 358   o -= (4./3.)*li2spence_(&shat);
 359   o -= (2./3.)*log(shat)*log(1.-shat);
 360   o -= log(1.-shat)*(5.+4.*shat)/(3.+6.*shat);
 361   o -= log(shat)*2.*shat*(1.+shat)*(1.-2.*shat)/3./(1.-shat)/(1.-shat)/(1.+2.*shat);
 362   o += (5.+9.*shat-6.*shat*shat)/6./(1.-shat)/(1.+2.*shat);
 363   return o;
 364 }
 365
 366 EvtComplex EvtWilsonCoeficients::C9efftilda(double z,double shat,double alpha_S,EvtComplex c1,EvtComplex c2,EvtComplex c3,EvtComplex c4,EvtComplex c5,EvtComplex c6,EvtComplex c9tilda,int ksi=0){
 367   EvtComplex c(0,0);
 368   c += (c9tilda+ksi*4./9.*(3.*c1+c2-c3-3.*c4))*etatilda(shat,alpha_S);
 369   c += hzs(z,shat)*(3.*c1+c2+3.*c3+c4+3.*c5+c6);
 370   c -= 0.5*hzs(1,shat)*(4.*c3+4.*c4+3.*c5+c6);
 371   c -= 0.5*hzs(0,shat)*(c3+3.*c4);
 372   c += 2./9.*(3.*c3+c4+3.*c5+c6);
 373   return c;
 374 }
 375
 376 EvtComplex EvtWilsonCoeficients::C7b2sg(double alpha_S,double et,EvtComplex c2,double M_t=174.3,double M_W=80.425){
 377   EvtComplex i1(0,1);
 378   return (i1*alpha_S*(2./9.*pow(et,14./23.)*(0.5*F(M_t*M_t/M_W/M_W)-0.1687)-0.03*c2));
 379 }
 380
 381 EvtComplex EvtWilsonCoeficients::Yld(double q2,double *ki,double *Gi,double *Mi,int ni,EvtComplex c1,EvtComplex c2,EvtComplex c3,EvtComplex c4,EvtComplex c5,EvtComplex c6,double ialpha=137.036){
 382   EvtComplex i1(0,1);
 383   EvtComplex y(0,0);
 384   int i;
 385   for(i=0;i<ni;i++) y+=ki[i]*Gi[i]*Mi[i]/(q2-Mi[i]*Mi[i]-i1*Mi[i]*Gi[i]);
 386   return (-3.*ialpha*ialpha*y*EvtConst::pi*(3.*c1+c2+3.*c3+c4+3.*c5+c6));
 387 }