EvtBtoXsllUtil Class Reference

#include <EvtBtoXsllUtil.hh>

List of all members.

Public Member Functions
EvtComplex	GetC7Eff0 (double sh, bool nnlo=true)
EvtComplex	GetC7Eff1 (double sh, double mb, bool nnlo=true)
EvtComplex	GetC9Eff0 (double sh, double mb, bool nnlo=true, bool btod=false)
EvtComplex	GetC9Eff1 (double sh, double mb, bool nnlo=true, bool btod=false)
EvtComplex	GetC10Eff (double sh, bool nnlo=true)
double	dGdsProb (double mb, double ms, double ml, double s)
double	dGdsdupProb (double mb, double ms, double ml, double s, double u)
double	FermiMomentum (double pf)
double	FermiMomentumProb (double pb, double pf)

---

Detailed Description

Definition at line 30 of file EvtBtoXsllUtil.hh.

---

Member Function Documentation

double EvtBtoXsllUtil::dGdsdupProb	(	double	mb,
		double	ms,
		double	ml,
		double	s,
		double	u
	)

Definition at line 481 of file EvtBtoXsllUtil.cc.

References abs2(), conj(), ddilog_(), GetC10Eff(), GetC7Eff0(), GetC7Eff1(), GetC9Eff0(), GetC9Eff1(), and EvtConst::pi.

Referenced by EvtBtoXsll::decay(), and EvtBtoXsll::init().

{
  // Compute the decay probability density function given a value of s and u
  // according to Ali-Hiller-Handoko-Morozumi's 1997 paper
  // see Appendix E

  bool btod = false;
  bool nnlo = true;

  double prob;
  double f1sp, f2sp, f3sp;
  //double u_ext;
  double mbeff = 4.8;

  //  double sh = s / (mb*mb);
  double sh  = s  / (mbeff*mbeff);

  EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo);
  EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo);
  EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod);
  EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod);
  EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo);

  double alphas = 0.119/
     (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);

  double omega7 = -8.0/3.0*log(4.8/mb)
                  -4.0/3.0*ddilog_(sh) 
                  -2.0/9.0*EvtConst::pi*EvtConst::pi
                  -2.0/3.0*log(sh)*log(1.0-sh)
                  -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0 
    -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
    -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
  double eta7 = 1.0 + alphas*omega7/EvtConst::pi;

  double omega79 = -4.0/3.0*log(4.8/mb)
                   -4.0/3.0*ddilog_(sh) 
                   -2.0/9.0*EvtConst::pi*EvtConst::pi
                   -2.0/3.0*log(sh)*log(1.0-sh)
                   -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
                   -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2) 
                   +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
  double eta79 = 1.0 + alphas*omega79/EvtConst::pi;

  double omega9 = - 2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*ddilog_(sh)
                 - 2.0/3.0*log(sh)*log(1.0-sh)
                 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
                 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
                 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
                 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
  double eta9 = 1.0 + alphas*omega9/EvtConst::pi;

  EvtComplex c7eff = eta7*c7eff0 + c7eff1;
  EvtComplex c9eff = eta9*c9eff0 + c9eff1;
  c10eff *= eta9;

  double c7c7  = abs2(c7eff);
  double c7c9  = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1));
  double c7c10 = real((eta79*c7eff0 + c7eff1)*conj(eta9*c10eff));
  double c9c10 = real((eta9*c9eff0  + c9eff1)*conj(eta9*c10eff));
  double c9c9plusc10c10  = abs2(c9eff) + abs2(c10eff);
  //double c9c9minusc10c10 = abs2(c9eff) - abs2(c10eff);

  f1sp = ( pow(mb*mb-ms*ms,2) - s*s) * c9c9plusc10c10 
         + 4.0*( pow(mb,4) - ms*ms*mb*mb - pow(ms,4)*(1.0 - ms*ms/(mb*mb))
         - 8.0*s*ms*ms - s*s*(1.0 + ms*ms/(mb*mb) ))*mb*mb*c7c7/s
    // kludged mass term
         *(1.0 + 2.0*ml*ml/s)
         - 8.0*(s*(mb*mb + ms*ms) - pow(mb*mb-ms*ms,2)) * c7c9
    // kludged mass term
         *(1.0 + 2.0*ml*ml/s);

  f2sp = 4.0*s*c9c10 + 8.0*(mb*mb + ms*ms)*c7c10;
  f3sp = - (c9c9plusc10c10)
         + 4.0*(1.0 + pow(ms/mb,4)) * mb*mb*c7c7/s
    // kludged mass term
         *(1.0 + 2.0*ml*ml/s);

  prob = (f1sp + f2sp*u + f3sp*u*u)/ pow(mb,3);

  return prob;
}

double EvtBtoXsllUtil::dGdsProb	(	double	mb,
		double	ms,
		double	ml,
		double	s
	)

Definition at line 352 of file EvtBtoXsllUtil.cc.

References abs2(), conj(), ddilog_(), F, f1, f2, GetC10Eff(), GetC7Eff0(), GetC7Eff1(), GetC9Eff0(), GetC9Eff1(), lambda, and EvtConst::pi.

Referenced by EvtBtoXsll::init().

{
  // Compute the decay probability density function given a value of s
  // according to Ali-Lunghi-Greub-Hiller's 2002 paper
  // Note that the form given below is taken from
  // F.Kruger and L.M.Sehgal, Phys. Lett. B380, 199 (1996)
  // but the differential rate as a function of dilepton mass
  // in this latter paper reduces to Eq.(12) in ALGH's 2002 paper
  // for ml = 0 and ms = 0.

  bool btod = false;
  bool nnlo = true;

  double delta, lambda, prob;
  double f1, f2, f3, f4;
  double msh, mlh, sh;
  double mbeff = 4.8;

  mlh = ml / mb;
  msh = ms / mb;
  // set lepton and strange-quark masses to 0 if need to
  // be in strict agreement with ALGH 2002 paper
  //  mlh = 0.0; msh = 0.0;
  //  sh  = s  / (mb*mb);
  sh  = s  / (mbeff*mbeff);

  EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo);
  EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo);
  EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod);
  EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod);
  EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo);

  double alphas = 0.119/
     (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);

  double omega7 = -8.0/3.0*log(4.8/mb)
                  -4.0/3.0*ddilog_(sh) 
                  -2.0/9.0*EvtConst::pi*EvtConst::pi
                  -2.0/3.0*log(sh)*log(1.0-sh)
                  -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0 
    -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
    -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
  double eta7 = 1.0 + alphas*omega7/EvtConst::pi;

  double omega79 = -4.0/3.0*log(4.8/mb)
                   -4.0/3.0*ddilog_(sh) 
                   -2.0/9.0*EvtConst::pi*EvtConst::pi
                   -2.0/3.0*log(sh)*log(1.0-sh)
                   -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
                   -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2) 
                   +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
  double eta79 = 1.0 + alphas*omega79/EvtConst::pi;

  double omega9 = -2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*ddilog_(sh)
                 - 2.0/3.0*log(sh)*log(1.0-sh)
                 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
                 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
                 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
                 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
  double eta9 = 1.0 + alphas*omega9/EvtConst::pi;

  EvtComplex c7eff = eta7*c7eff0 + c7eff1;
  EvtComplex c9eff = eta9*c9eff0 + c9eff1;
  c10eff *= eta9;

  double c7c7 = abs2(c7eff);
  double c7c9 = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1));
  double c9c9plusc10c10  = abs2(c9eff) + abs2(c10eff);
  double c9c9minusc10c10 = abs2(c9eff) - abs2(c10eff);

  // Power corrections according to ALGH 2002
  double lambda_1 = -0.2;
  double lambda_2 = 0.12;
  double C1 = -0.487;
  double C2 = 1.024;
  double mc = 0.29 * mb;

  EvtComplex F;
  double r = s / (4.0 * mc * mc);
  EvtComplex uniti(0.0,1.0);
  F = 3.0 / (2.0 * r);
  if (r < 1)
  {
    F *= 1.0/sqrt(r*(1.0-r))*atan(sqrt(r/(1.0-r)))-1.0;
  }
  else
  {
    F *= 0.5/sqrt(r*(r-1.0))*(log((1.0-sqrt(1.0-1.0/r))/(1.0+sqrt(1.0-1.0/r)))
                              +uniti*EvtConst::pi)-1.0;
  }

  double G1 = 1.0 + lambda_1 / (2.0 * mb * mb)
                  + 3.0 * (1.0 - 15.0*sh*sh + 10.0*sh*sh*sh)
                        / ((1.0 - sh)*(1.0 -sh)*(1.0 + 2.0*sh))
                        * lambda_2 / (2.0*mb*mb);
  double G2 = 1.0 + lambda_1 / (2.0 * mb * mb)
                  - 3.0 * (6.0 + 3.0*sh - 5.0*sh*sh*sh)
                        / ((1.0 - sh)*(1.0 -sh)*(2.0 + sh))
                        * lambda_2 / (2.0*mb*mb);
  double G3 = 1.0 + lambda_1 / (2.0 * mb * mb)
                  - (5.0 + 6.0*sh - 7.0*sh*sh)
                     / ((1.0 - sh)*(1.0 -sh))
                     * lambda_2 / (2.0*mb*mb);
  double Gc = -8.0/9.0 * (C2 - C1/6.0) * lambda_2/(mc*mc) 
    * real(F*(conj(c9eff)*(2.0+sh)+conj(c7eff)*(1.0 + 6.0*sh - sh*sh)/sh));

  // end of power corrections section
  // now back to Kruger & Sehgal expressions

  lambda = 1.0 + sh*sh + pow(msh,4) - 2.0*(sh + sh*msh*msh + msh*msh);

  f1 = pow(1.0-msh*msh,2) - sh*(1.0 + msh*msh);
  f2 = 2.0*(1.0 + msh*msh) * pow(1.0-msh*msh,2)
       - sh*(1.0 + 14.0*msh*msh + pow(msh,4)) - sh*sh*(1.0 + msh*msh);
  f3 = pow(1.0-msh*msh,2) + sh*(1.0 + msh*msh) - 2.0*sh*sh
       + lambda*2.0*mlh*mlh/sh;
  f4 = 1.0 - sh + msh*msh;

  delta = (  12.0*c7c9*f1*G3 + 4.0*c7c7*f2*G2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
            + c9c9plusc10c10*f3*G1 
            + 6.0*mlh*mlh*c9c9minusc10c10*f4
            + Gc;

  prob =  sqrt(lambda*(1.0 - 4.0*mlh*mlh/sh)) * delta;

  return prob;
}

double EvtBtoXsllUtil::FermiMomentum

(

double

pf

)

Definition at line 565 of file EvtBtoXsllUtil.cc.

References FermiMomentumProb(), and EvtRandom::Flat().

Referenced by EvtBtoXsll::decay().

{
  // Pick a value for the b-quark Fermi motion momentum
  // according to Ali's Gaussian model

  double pb, pbmax, xbox, ybox;
  pb    = 0.0;
  pbmax = 5.0 * pf;

  while (pb == 0.0)
  {
    xbox = EvtRandom::Flat(pbmax);
    ybox = EvtRandom::Flat();
    if (ybox < FermiMomentumProb(xbox, pf)) { pb = xbox;}
  }

  return pb;
}

double EvtBtoXsllUtil::FermiMomentumProb	(	double	pb,
		double	pf
	)

Definition at line 584 of file EvtBtoXsllUtil.cc.

References exp().

Referenced by FermiMomentum().

{
  // Compute probability according to Ali's Gaussian model
  // the function chosen has a convenient maximum value of 1 for pb = pf

  double prsq = (pb*pb)/(pf*pf);
  double prob = prsq * exp(1.0 - prsq);

  return prob;
}

EvtComplex EvtBtoXsllUtil::GetC10Eff	(	double	sh,
		bool	nnlo = true
	)

Definition at line 339 of file EvtBtoXsllUtil.cc.

Referenced by dGdsdupProb(), and dGdsProb().

{

  if (!nnlo) return -4.669;
  double A10;
  A10 = -4.592 + 0.379;

  EvtComplex c10eff;
  c10eff = A10;

  return c10eff;
}

EvtComplex EvtBtoXsllUtil::GetC7Eff0	(	double	sh,
		bool	nnlo = true
	)

Definition at line 43 of file EvtBtoXsllUtil.cc.

Referenced by dGdsdupProb(), and dGdsProb().

{
  // This function returns the zeroth-order alpha_s part of C7

  if (!nnlo) return -0.313;

  double A7;

  // use energy scale of 2.5 GeV as a computational trick (G.Hiller)
  // at least for shat > 0.25
  A7 = -0.353 + 0.023;

  EvtComplex c7eff;
  if (sh > 0.25)
  { 
    c7eff = A7;
    return c7eff;
  }

  // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
  A7 = -0.312 + 0.008;
  c7eff = A7;

  return c7eff;
}

EvtComplex EvtBtoXsllUtil::GetC7Eff1	(	double	sh,
		double	mb,
		bool	nnlo = true
	)

Definition at line 69 of file EvtBtoXsllUtil.cc.

References EvtConst::pi.

Referenced by dGdsdupProb(), and dGdsProb().

{
  // This function returns the first-order alpha_s part of C7

  if (!nnlo) return 0.0;
  double logsh;
  logsh = log(sh);

  EvtComplex uniti(0.0,1.0);

  EvtComplex c7eff = 0.0;
  if (sh > 0.25)
  { 
    return c7eff;
  }

  // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
  double muscale = 5.0;
  double alphas = 0.215;
  //double A7 = -0.312 + 0.008;
  double A8 = -0.148;
  //double A9 = 4.174 + (-0.035);
  //double A10 = -4.592 + 0.379;
  double C1 = -0.487;
  double C2 = 1.024;
  //double T9 = 0.374 + 0.252;
  //double U9 = 0.033 + 0.015;
  //double W9 = 0.032 + 0.012;
  double Lmu = log(muscale/mbeff);

  EvtComplex F71;
  EvtComplex f71;
  EvtComplex k7100(-0.68192,-0.074998);
  EvtComplex k7101(0.0,0.0);
  EvtComplex k7110(-0.23935,-0.12289);
  EvtComplex k7111(0.0027424,0.019676);
  EvtComplex k7120(-0.0018555,-0.175);
  EvtComplex k7121(0.022864,0.011456);
  EvtComplex k7130(0.28248,-0.12783);
  EvtComplex k7131(0.029027,-0.0082265);
  f71 = k7100 + k7101*logsh + sh*(k7110 + k7111*logsh) +
        sh*sh*(k7120 + k7121*logsh) + 
        sh*sh*sh*(k7130 + k7131*logsh); 
  F71 = (-208.0/243.0)*Lmu + f71;

  EvtComplex F72;
  EvtComplex f72;
  EvtComplex k7200(4.0915,0.44999);
  EvtComplex k7201(0.0,0.0);
  EvtComplex k7210(1.4361,0.73732);
  EvtComplex k7211(-0.016454,-0.11806);
  EvtComplex k7220(0.011133,1.05);
  EvtComplex k7221(-0.13718,-0.068733);
  EvtComplex k7230(-1.6949,0.76698);
  EvtComplex k7231(-0.17416,0.049359);
  f72 = k7200 + k7201*logsh + sh*(k7210 + k7211*logsh) +
        sh*sh*(k7220 + k7221*logsh) + 
        sh*sh*sh*(k7230 + k7231*logsh); 
  F72 = (416.0/81.0)*Lmu + f72;
  
  EvtComplex F78;
  F78 = (-32.0/9.0)*Lmu + 8.0*EvtConst::pi*EvtConst::pi/27.0 + (-44.0/9.0) 
        + (-8.0*EvtConst::pi/9.0)*uniti +
        (4.0/3.0*EvtConst::pi*EvtConst::pi - 40.0/3.0)*sh +
        (32.0*EvtConst::pi*EvtConst::pi/9.0 - 316.0/9.0)*sh*sh +
        (200.0*EvtConst::pi*EvtConst::pi/27.0 - 658.0/9.0)*sh*sh*sh +
    (-8.0*logsh/9.0)*(sh + sh*sh + sh*sh*sh);
        
  c7eff = - alphas/(4.0*EvtConst::pi)*(C1*F71 + C2*F72 + A8*F78);

  return c7eff;
}

EvtComplex EvtBtoXsllUtil::GetC9Eff0	(	double	sh,
		double	mb,
		bool	nnlo = true,
		bool	btod = false
	)

Definition at line 143 of file EvtBtoXsllUtil.cc.

References rb::hc, and EvtConst::pi.

Referenced by dGdsdupProb(), and dGdsProb().

{
  // This function returns the zeroth-order alpha_s part of C9

  if (!nnlo) return 4.344;
  double logsh;
  logsh = log(sh);
  double mch = 0.29;

  
  double muscale;
  muscale = 2.5;
  double alphas;
  alphas = 0.267;
  double A8;
  A8 = -0.164;
  double A9;
  A9 = 4.287 + (-0.218);
  double A10;
  A10 = -4.592 + 0.379;
  double C1;
  C1 = -0.697;
  double C2;
  C2 = 1.046;
  double T9;
  T9 = 0.114 + 0.280;
  double U9;
  U9 = 0.045 + 0.023;
  double W9;
  W9 = 0.044 + 0.016;
  
  double Lmu;
  Lmu = log(muscale/mbeff);


  EvtComplex uniti(0.0,1.0);

  EvtComplex hc;
  double xarg;
  xarg = 4.0*mch/sh;
  hc = -4.0/9.0*log(mch*mch) + 8.0/27.0 + 4.0*xarg/9.0;
  if (xarg < 1.0)
  { 
    hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
      (log((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0)) - 
       uniti*EvtConst::pi);
  } 
  else
  {
    hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
      2.0*atan(1.0/sqrt(xarg - 1.0));
  }

  EvtComplex h1;
  xarg = 4.0/sh;
  h1 = 8.0/27.0 + 4.0*xarg/9.0;
  if (xarg < 1.0)
  { 
    h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
      (log((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0)) - 
       uniti*EvtConst::pi);
  } 
  else
  {
    h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
      2.0*atan(1.0/sqrt(xarg - 1.0));
  }

  EvtComplex h0;
  h0 = 8.0/27.0 - 4.0*log(2.0)/9.0 + 4.0*uniti*EvtConst::pi/9.0;


  // X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
  // h(\hat m_u^2, hat s))
  EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0);
  EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0);
  EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0);
  EvtComplex Vtb(1.0,0.0);

  EvtComplex Xd;
  Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);

  EvtComplex c9eff = 4.344;
  if (sh > 0.25)
  { 
    c9eff =  A9 + T9*hc + U9*h1 + W9*h0;
    if (btod)
    {
      c9eff += Xd; 
    }
    return c9eff;
  }

  // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
  muscale = 5.0;
  alphas = 0.215;
  A9 = 4.174 + (-0.035);
  C1 = -0.487;
  C2 = 1.024;
  A8 = -0.148;
  T9 = 0.374 + 0.252;
  U9 = 0.033 + 0.015;
  W9 = 0.032 + 0.012;
  Lmu = log(muscale/mbeff);

  Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);

  c9eff = A9 + T9*hc + U9*h1 + W9*h0;

  if (btod)
  {
    c9eff += Xd; 
  }

  return c9eff;
}

EvtComplex EvtBtoXsllUtil::GetC9Eff1	(	double	sh,
		double	mb,
		bool	nnlo = true,
		bool	btod = false
	)

Definition at line 261 of file EvtBtoXsllUtil.cc.

References EvtConst::pi.

Referenced by dGdsdupProb(), and dGdsProb().

{
  // This function returns the first-order alpha_s part of C9

  if (!nnlo) return 0.0;
  double logsh;
  logsh = log(sh);
  double mch = 0.29;

  EvtComplex uniti(0.0,1.0);

  EvtComplex c9eff = 0.0;
  if (sh > 0.25)
  { 
    return c9eff;
  }

  // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
  double muscale = 5.0;
  double alphas = 0.215;
  double C1 = -0.487;
  double C2 = 1.024;
  double A8 = -0.148;
  double Lmu = log(muscale/mbeff);

  EvtComplex F91;
  EvtComplex f91;
  EvtComplex k9100(-11.973,0.16371);
  EvtComplex k9101(-0.081271,-0.059691);
  EvtComplex k9110(-28.432,-0.25044);
  EvtComplex k9111(-0.040243,0.016442);
  EvtComplex k9120(-57.114,-0.86486);
  EvtComplex k9121(-0.035191,0.027909);
  EvtComplex k9130(-128.8,-2.5243);
  EvtComplex k9131(-0.017587,0.050639);
  f91 = k9100 + k9101*logsh + sh*(k9110 + k9111*logsh) +
        sh*sh*(k9120 + k9121*logsh) + 
        sh*sh*sh*(k9130 + k9131*logsh); 
  F91 = (-1424.0/729.0 + 16.0*uniti*EvtConst::pi/243.0 
         + 64.0/27.0*log(mch))*Lmu - 16.0*Lmu*logsh/243.0 +
        (16.0/1215.0 - 32.0/135.0/mch/mch)*Lmu*sh +
        (4.0/2835.0 - 8.0/315.0/mch/mch/mch/mch)*Lmu*sh*sh +
    (16.0/76545.0 - 32.0/8505.0/mch/mch/mch/mch/mch/mch)*
    Lmu*sh*sh*sh -256.0*Lmu*Lmu/243.0 + f91;

  EvtComplex F92;
  EvtComplex f92;
  EvtComplex k9200(6.6338,-0.98225);
  EvtComplex k9201(0.48763,0.35815);
  EvtComplex k9210(3.3585,1.5026);
  EvtComplex k9211(0.24146,-0.098649);
  EvtComplex k9220(-1.1906,5.1892);
  EvtComplex k9221(0.21115,-0.16745);
  EvtComplex k9230(-17.12,15.146);
  EvtComplex k9231(0.10552,-0.30383);
  f92 = k9200 + k9201*logsh + sh*(k9210 + k9211*logsh) +
        sh*sh*(k9220 + k9221*logsh) + 
        sh*sh*sh*(k9230 + k9231*logsh); 
  F92 = (256.0/243.0 - 32.0*uniti*EvtConst::pi/81.0 
         - 128.0/9.0*log(mch))*Lmu + 32.0*Lmu*logsh/81.0 +
        (-32.0/405.0 + 64.0/45.0/mch/mch)*Lmu*sh +
        (-8.0/945.0 + 16.0/105.0/mch/mch/mch/mch)*Lmu*sh*sh +
    (-32.0/25515.0 + 64.0/2835.0/mch/mch/mch/mch/mch/mch)*
    Lmu*sh*sh*sh + 512.0*Lmu*Lmu/81.0 + f92;
  
  EvtComplex F98;
  F98 = 104.0/9.0 - 32.0*EvtConst::pi*EvtConst::pi/27.0 + 
        (1184.0/27.0 - 40.0*EvtConst::pi*EvtConst::pi/9.0)*sh +
        (14212.0/135.0 - 32.0*EvtConst::pi*EvtConst::pi/3.0)*sh*sh +
    (193444.0/945.0 - 560.0*EvtConst::pi*EvtConst::pi/27.0)*sh*sh*sh +
        16.0*logsh/9.0*(1.0 + sh + sh*sh + sh*sh*sh);

  c9eff = - alphas/(4.0*EvtConst::pi)*(C1*F91 + C2*F92 + A8*F98);

  return c9eff;
}

---