TEvtGen/EvtGenModels/EvtBtoXsllUtil.cxx

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Module: EvtBtoXsllUtil.cc
   9 //
  10 // Description: Routine to generate non-resonant B -> Xs l+ l- decays.
  11 // It generates a dilepton mass spectrum according to
  12 // F.Kruger and L.M.Sehgal, Phys. Lett. B380, 199 (1996)
  13 // and then generates the two lepton momenta according to
  14 // A.Ali, G.Hiller, L.T.Handoko and T.Morozumi, Phys. Rev. D55, 4105 (1997).
  15 // Expressions for Wilson coefficients and power corrections are taken
  16 // from A.Ali, E.Lunghi, C.Greub and G.Hiller, Phys. Rev. D66, 034002 (2002).
  17 // Detailed formulae for shat dependence of these coefficients are taken
  18 // from H.H.Asatryan, H.M.Asatrian, C.Greub and M.Walker, PRD65, 074004 (2002)
  19 // and C.Bobeth, M.Misiak and J.Urban, Nucl. Phys. B574, 291 (2000).
  20 // The resultant Xs particles may be decayed by JETSET.
  21 //
  22 // Modification history:
  23 //
  24 //    Stephane Willocq    Jan 19, 2001   Module created
  25 //    Stephane Willocq    Nov  6, 2003   Update Wilson Coeffs & dG's
  26 //    &Jeff Berryhill
  27 //
  28 //------------------------------------------------------------------------
  29 //
  30 #include "EvtGenBase/EvtPatches.hh"
  31 //
  32 #include <stdlib.h>
  33 #include "EvtGenBase/EvtRandom.hh"
  34 #include "EvtGenBase/EvtParticle.hh"
  35 #include "EvtGenBase/EvtGenKine.hh"
  36 #include "EvtGenBase/EvtPDL.hh"
  37 #include "EvtGenBase/EvtReport.hh"
  38 #include "EvtGenModels/EvtBtoXsllUtil.hh"
  39 #include "EvtGenBase/EvtComplex.hh"
  40 #include "EvtGenBase/EvtConst.hh"
  41 #include "EvtGenBase/EvtDiLog.hh"
  42
  43 EvtComplex EvtBtoXsllUtil::GetC7Eff0(double sh, bool nnlo)
  44 {
  45   // This function returns the zeroth-order alpha_s part of C7
  46
  47   if (!nnlo) return -0.313;
  48
  49   double A7;
  50
  51   // use energy scale of 2.5 GeV as a computational trick (G.Hiller)
  52   // at least for shat > 0.25
  53   A7 = -0.353 + 0.023;
  54
  55   EvtComplex c7eff;
  56   if (sh > 0.25)
  57   {
  58     c7eff = A7;
  59     return c7eff;
  60   }
  61
  62   // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
  63   A7 = -0.312 + 0.008;
  64   c7eff = A7;
  65
  66   return c7eff;
  67 }
  68
  69 EvtComplex EvtBtoXsllUtil::GetC7Eff1(double sh, double mbeff, bool nnlo)
  70 {
  71   // This function returns the first-order alpha_s part of C7
  72
  73   if (!nnlo) return 0.0;
  74   double logsh;
  75   logsh = log(sh);
  76
  77   EvtComplex uniti(0.0,1.0);
  78
  79   EvtComplex c7eff = 0.0;
  80   if (sh > 0.25)
  81   {
  82     return c7eff;
  83   }
  84
  85   // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
  86   double muscale = 5.0;
  87   double alphas = 0.215;
  88   //double A7 = -0.312 + 0.008;
  89   double A8 = -0.148;
  90   //double A9 = 4.174 + (-0.035);
  91   //double A10 = -4.592 + 0.379;
  92   double C1 = -0.487;
  93   double C2 = 1.024;
  94   //double T9 = 0.374 + 0.252;
  95   //double U9 = 0.033 + 0.015;
  96   //double W9 = 0.032 + 0.012;
  97   double Lmu = log(muscale/mbeff);
  98
  99   EvtComplex F71;
 100   EvtComplex f71;
 101   EvtComplex k7100(-0.68192,-0.074998);
 102   EvtComplex k7101(0.0,0.0);
 103   EvtComplex k7110(-0.23935,-0.12289);
 104   EvtComplex k7111(0.0027424,0.019676);
 105   EvtComplex k7120(-0.0018555,-0.175);
 106   EvtComplex k7121(0.022864,0.011456);
 107   EvtComplex k7130(0.28248,-0.12783);
 108   EvtComplex k7131(0.029027,-0.0082265);
 109   f71 = k7100 + k7101*logsh + sh*(k7110 + k7111*logsh) +
 110         sh*sh*(k7120 + k7121*logsh) +
 111         sh*sh*sh*(k7130 + k7131*logsh);
 112   F71 = (-208.0/243.0)*Lmu + f71;
 113
 114   EvtComplex F72;
 115   EvtComplex f72;
 116   EvtComplex k7200(4.0915,0.44999);
 117   EvtComplex k7201(0.0,0.0);
 118   EvtComplex k7210(1.4361,0.73732);
 119   EvtComplex k7211(-0.016454,-0.11806);
 120   EvtComplex k7220(0.011133,1.05);
 121   EvtComplex k7221(-0.13718,-0.068733);
 122   EvtComplex k7230(-1.6949,0.76698);
 123   EvtComplex k7231(-0.17416,0.049359);
 124   f72 = k7200 + k7201*logsh + sh*(k7210 + k7211*logsh) +
 125         sh*sh*(k7220 + k7221*logsh) +
 126         sh*sh*sh*(k7230 + k7231*logsh);
 127   F72 = (416.0/81.0)*Lmu + f72;
 128
 129   EvtComplex F78;
 130   F78 = (-32.0/9.0)*Lmu + 8.0*EvtConst::pi*EvtConst::pi/27.0 + (-44.0/9.0)
 131         + (-8.0*EvtConst::pi/9.0)*uniti +
 132         (4.0/3.0*EvtConst::pi*EvtConst::pi - 40.0/3.0)*sh +
 133         (32.0*EvtConst::pi*EvtConst::pi/9.0 - 316.0/9.0)*sh*sh +
 134         (200.0*EvtConst::pi*EvtConst::pi/27.0 - 658.0/9.0)*sh*sh*sh +
 135     (-8.0*logsh/9.0)*(sh + sh*sh + sh*sh*sh);
 136
 137   c7eff = - alphas/(4.0*EvtConst::pi)*(C1*F71 + C2*F72 + A8*F78);
 138
 139   return c7eff;
 140 }
 141
 142
 143 EvtComplex EvtBtoXsllUtil::GetC9Eff0(double sh, double mbeff,
 144                                      bool nnlo, bool btod)
 145 {
 146   // This function returns the zeroth-order alpha_s part of C9
 147
 148   if (!nnlo) return 4.344;
 149   double logsh;
 150   logsh = log(sh);
 151   double mch = 0.29;
 152
 153
 154   double muscale;
 155   muscale = 2.5;
 156   double alphas;
 157   alphas = 0.267;
 158   double A8;
 159   A8 = -0.164;
 160   double A9;
 161   A9 = 4.287 + (-0.218);
 162   double A10;
 163   A10 = -4.592 + 0.379;
 164   double C1;
 165   C1 = -0.697;
 166   double C2;
 167   C2 = 1.046;
 168   double T9;
 169   T9 = 0.114 + 0.280;
 170   double U9;
 171   U9 = 0.045 + 0.023;
 172   double W9;
 173   W9 = 0.044 + 0.016;
 174
 175   double Lmu;
 176   Lmu = log(muscale/mbeff);
 177
 178
 179   EvtComplex uniti(0.0,1.0);
 180
 181   EvtComplex hc;
 182   double xarg;
 183   xarg = 4.0*mch/sh;
 184
 185   hc = -4.0/9.0*log(mch*mch) + 8.0/27.0 + 4.0*xarg/9.0;
 186   if (xarg < 1.0)
 187   {
 188     hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 189       (log((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0)) -
 190        uniti*EvtConst::pi);
 191   }
 192   else
 193   {
 194     hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 195       2.0*atan(1.0/sqrt(xarg-1.0));
 196   }
 197
 198   EvtComplex h1;
 199   xarg = 4.0/sh;
 200   h1 = 8.0/27.0 + 4.0*xarg/9.0;
 201   if (xarg < 1.0)
 202   {
 203     h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 204       (log((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0)) -
 205        uniti*EvtConst::pi);
 206   }
 207   else
 208   {
 209     h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 210       2.0*atan(1.0/sqrt(xarg-1.0));
 211   }
 212
 213   EvtComplex h0;
 214   h0 = 8.0/27.0 - 4.0*log(2.0)/9.0 + 4.0*uniti*EvtConst::pi/9.0;
 215
 216
 217   // X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
 218   // h(\hat m_u^2, hat s))
 219   EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0);
 220   EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0);
 221   EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0);
 222   EvtComplex Vtb(1.0,0.0);
 223
 224   EvtComplex Xd;
 225   Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
 226
 227   EvtComplex c9eff = 4.344;
 228   if (sh > 0.25)
 229   {
 230     c9eff =  A9 + T9*hc + U9*h1 + W9*h0;
 231     if (btod)
 232     {
 233       c9eff += Xd;
 234     }
 235     return c9eff;
 236   }
 237
 238   // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
 239   muscale = 5.0;
 240   alphas = 0.215;
 241   A9 = 4.174 + (-0.035);
 242   C1 = -0.487;
 243   C2 = 1.024;
 244   A8 = -0.148;
 245   T9 = 0.374 + 0.252;
 246   U9 = 0.033 + 0.015;
 247   W9 = 0.032 + 0.012;
 248   Lmu = log(muscale/mbeff);
 249
 250   Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
 251
 252   c9eff = A9 + T9*hc + U9*h1 + W9*h0;
 253
 254   if (btod)
 255   {
 256     c9eff += Xd;
 257   }
 258
 259   return c9eff;
 260 }
 261
 262 EvtComplex EvtBtoXsllUtil::GetC9Eff1(double sh, double mbeff,
 263                                      bool nnlo, bool /*btod*/)
 264 {
 265   // This function returns the first-order alpha_s part of C9
 266
 267   if (!nnlo) return 0.0;
 268   double logsh;
 269   logsh = log(sh);
 270   double mch = 0.29;
 271
 272   EvtComplex uniti(0.0,1.0);
 273
 274   EvtComplex c9eff = 0.0;
 275   if (sh > 0.25)
 276   {
 277     return c9eff;
 278   }
 279
 280   // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
 281   double muscale = 5.0;
 282   double alphas = 0.215;
 283   double C1 = -0.487;
 284   double C2 = 1.024;
 285   double A8 = -0.148;
 286   double Lmu = log(muscale/mbeff);
 287
 288   EvtComplex F91;
 289   EvtComplex f91;
 290   EvtComplex k9100(-11.973,0.16371);
 291   EvtComplex k9101(-0.081271,-0.059691);
 292   EvtComplex k9110(-28.432,-0.25044);
 293   EvtComplex k9111(-0.040243,0.016442);
 294   EvtComplex k9120(-57.114,-0.86486);
 295   EvtComplex k9121(-0.035191,0.027909);
 296   EvtComplex k9130(-128.8,-2.5243);
 297   EvtComplex k9131(-0.017587,0.050639);
 298   f91 = k9100 + k9101*logsh + sh*(k9110 + k9111*logsh) +
 299         sh*sh*(k9120 + k9121*logsh) +
 300         sh*sh*sh*(k9130 + k9131*logsh);
 301   F91 = (-1424.0/729.0 + 16.0*uniti*EvtConst::pi/243.0
 302          + 64.0/27.0*log(mch))*Lmu - 16.0*Lmu*logsh/243.0 +
 303         (16.0/1215.0 - 32.0/135.0/mch/mch)*Lmu*sh +
 304         (4.0/2835.0 - 8.0/315.0/mch/mch/mch/mch)*Lmu*sh*sh +
 305     (16.0/76545.0 - 32.0/8505.0/mch/mch/mch/mch/mch/mch)*
 306     Lmu*sh*sh*sh -256.0*Lmu*Lmu/243.0 + f91;
 307
 308   EvtComplex F92;
 309   EvtComplex f92;
 310   EvtComplex k9200(6.6338,-0.98225);
 311   EvtComplex k9201(0.48763,0.35815);
 312   EvtComplex k9210(3.3585,1.5026);
 313   EvtComplex k9211(0.24146,-0.098649);
 314   EvtComplex k9220(-1.1906,5.1892);
 315   EvtComplex k9221(0.21115,-0.16745);
 316   EvtComplex k9230(-17.12,15.146);
 317   EvtComplex k9231(0.10552,-0.30383);
 318   f92 = k9200 + k9201*logsh + sh*(k9210 + k9211*logsh) +
 319         sh*sh*(k9220 + k9221*logsh) +
 320         sh*sh*sh*(k9230 + k9231*logsh);
 321   F92 = (256.0/243.0 - 32.0*uniti*EvtConst::pi/81.0
 322          - 128.0/9.0*log(mch))*Lmu + 32.0*Lmu*logsh/81.0 +
 323         (-32.0/405.0 + 64.0/45.0/mch/mch)*Lmu*sh +
 324         (-8.0/945.0 + 16.0/105.0/mch/mch/mch/mch)*Lmu*sh*sh +
 325     (-32.0/25515.0 + 64.0/2835.0/mch/mch/mch/mch/mch/mch)*
 326     Lmu*sh*sh*sh + 512.0*Lmu*Lmu/81.0 + f92;
 327
 328   EvtComplex F98;
 329   F98 = 104.0/9.0 - 32.0*EvtConst::pi*EvtConst::pi/27.0 +
 330         (1184.0/27.0 - 40.0*EvtConst::pi*EvtConst::pi/9.0)*sh +
 331         (14212.0/135.0 - 32.0*EvtConst::pi*EvtConst::pi/3.0)*sh*sh +
 332     (193444.0/945.0 - 560.0*EvtConst::pi*EvtConst::pi/27.0)*sh*sh*sh +
 333         16.0*logsh/9.0*(1.0 + sh + sh*sh + sh*sh*sh);
 334
 335   c9eff = - alphas/(4.0*EvtConst::pi)*(C1*F91 + C2*F92 + A8*F98);
 336
 337   return c9eff;
 338 }
 339
 340 EvtComplex EvtBtoXsllUtil::GetC10Eff(double /*sh*/, bool nnlo)
 341 {
 342
 343   if (!nnlo) return -4.669;
 344   double A10;
 345   A10 = -4.592 + 0.379;
 346
 347   EvtComplex c10eff;
 348   c10eff = A10;
 349
 350   return c10eff;
 351 }
 352
 353 double EvtBtoXsllUtil::dGdsProb(double mb, double ms, double ml,
 354                                 double s)
 355 {
 356   // Compute the decay probability density function given a value of s
 357   // according to Ali-Lunghi-Greub-Hiller's 2002 paper
 358   // Note that the form given below is taken from
 359   // F.Kruger and L.M.Sehgal, Phys. Lett. B380, 199 (1996)
 360   // but the differential rate as a function of dilepton mass
 361   // in this latter paper reduces to Eq.(12) in ALGH's 2002 paper
 362   // for ml = 0 and ms = 0.
 363
 364   bool btod = false;
 365   bool nnlo = true;
 366
 367   double delta, lambda, prob;
 368   double f1, f2, f3, f4;
 369   double msh, mlh, sh;
 370   double mbeff = 4.8;
 371
 372   mlh = ml / mb;
 373   msh = ms / mb;
 374   // set lepton and strange-quark masses to 0 if need to
 375   // be in strict agreement with ALGH 2002 paper
 376   //  mlh = 0.0; msh = 0.0;
 377   //  sh  = s  / (mb*mb);
 378   sh  = s  / (mbeff*mbeff);
 379
 380   // if sh >1.0 code will return a nan. so just skip it
 381   if ( sh > 1.0 ) return 0.0;
 382
 383
 384   EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo);
 385   EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo);
 386   EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod);
 387   EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod);
 388   EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo);
 389
 390   double alphas = 0.119/
 391      (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
 392
 393   double omega7 = -8.0/3.0*log(4.8/mb)
 394                   -4.0/3.0*EvtDiLog::DiLog(sh)
 395                   -2.0/9.0*EvtConst::pi*EvtConst::pi
 396                   -2.0/3.0*log(sh)*log(1.0-sh)
 397                   -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
 398     -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
 399     -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
 400   double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
 401
 402   double omega79 = -4.0/3.0*log(4.8/mb)
 403                    -4.0/3.0*EvtDiLog::DiLog(sh)
 404                    -2.0/9.0*EvtConst::pi*EvtConst::pi
 405                    -2.0/3.0*log(sh)*log(1.0-sh)
 406                    -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
 407                    -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
 408                    +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
 409   double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
 410
 411   double omega9 = -2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
 412                  - 2.0/3.0*log(sh)*log(1.0-sh)
 413                  - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
 414                  - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
 415                  /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
 416                  + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
 417   double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
 418
 419   EvtComplex c7eff = eta7*c7eff0 + c7eff1;
 420   EvtComplex c9eff = eta9*c9eff0 + c9eff1;
 421   c10eff *= eta9;
 422
 423   double c7c7 = abs2(c7eff);
 424   double c7c9 = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1));
 425   double c9c9plusc10c10  = abs2(c9eff) + abs2(c10eff);
 426   double c9c9minusc10c10 = abs2(c9eff) - abs2(c10eff);
 427
 428   // Power corrections according to ALGH 2002
 429   double lambda_1 = -0.2;
 430   double lambda_2 = 0.12;
 431   double C1 = -0.487;
 432   double C2 = 1.024;
 433   double mc = 0.29 * mb;
 434
 435   EvtComplex F;
 436   double r = s / (4.0 * mc * mc);
 437   EvtComplex uniti(0.0,1.0);
 438   F = 3.0 / (2.0 * r);
 439   if (r < 1)
 440   {
 441     F *= 1.0/sqrt(r*(1.0-r))*atan(sqrt(r/(1.0-r)))-1.0;
 442   }
 443   else
 444   {
 445     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)))
 446                               +uniti*EvtConst::pi)-1.0;
 447   }
 448
 449   double G1 = 1.0 + lambda_1 / (2.0 * mb * mb)
 450                   + 3.0 * (1.0 - 15.0*sh*sh + 10.0*sh*sh*sh)
 451                         / ((1.0 - sh)*(1.0 -sh)*(1.0 + 2.0*sh))
 452                         * lambda_2 / (2.0*mb*mb);
 453   double G2 = 1.0 + lambda_1 / (2.0 * mb * mb)
 454                   - 3.0 * (6.0 + 3.0*sh - 5.0*sh*sh*sh)
 455                         / ((1.0 - sh)*(1.0 -sh)*(2.0 + sh))
 456                         * lambda_2 / (2.0*mb*mb);
 457   double G3 = 1.0 + lambda_1 / (2.0 * mb * mb)
 458                   - (5.0 + 6.0*sh - 7.0*sh*sh)
 459                      / ((1.0 - sh)*(1.0 -sh))
 460                      * lambda_2 / (2.0*mb*mb);
 461   double Gc = -8.0/9.0 * (C2 - C1/6.0) * lambda_2/(mc*mc)
 462     * real(F*(conj(c9eff)*(2.0+sh)+conj(c7eff)*(1.0 + 6.0*sh - sh*sh)/sh));
 463
 464   // end of power corrections section
 465   // now back to Kruger & Sehgal expressions
 466
 467   double msh2=msh*msh;
 468   lambda = 1.0 + sh*sh + msh2*msh2 - 2.0*(sh + sh*msh2 + msh2);
 469   // negative lambda screw up sqrt below!
 470   if ( lambda < 0.0 ) return 0.0;
 471
 472   f1 = pow(1.0-msh2,2) - sh*(1.0 + msh2);
 473   f2 = 2.0*(1.0 + msh2) * pow(1.0-msh2,2)
 474        - sh*(1.0 + 14.0*msh2 + pow(msh,4)) - sh*sh*(1.0 + msh2);
 475   f3 = pow(1.0-msh2,2) + sh*(1.0 + msh2) - 2.0*sh*sh
 476        + lambda*2.0*mlh*mlh/sh;
 477   f4 = 1.0 - sh + msh2;
 478
 479   delta = (  12.0*c7c9*f1*G3 + 4.0*c7c7*f2*G2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
 480             + c9c9plusc10c10*f3*G1
 481             + 6.0*mlh*mlh*c9c9minusc10c10*f4
 482             + Gc;
 483
 484   // avoid negative probs
 485   if ( delta < 0.0 ) delta=0.;
 486   // negative when sh < 4*mlh*mlh
 487   //               s < 4*ml*ml
 488   ///  prob =  sqrt(lambda*(1.0 - 4.0*mlh*mlh/sh)) * delta;
 489   prob =  sqrt(lambda*(1.0 - 4.0*ml*ml/s)) * delta;
 490
 491   //   if ( !(prob>=0.0) && !(prob<=0.0) ) {
 492     //nan
 493      //     std::cout << lambda << " " << mlh << " " << sh << " " << delta << " " << mb << " " << mbeff << std::endl;
 494      // std::cout << 4.0*mlh*mlh/sh << " " << 4.0*ml*ml/s <<  " " << s-4.0*ml*ml << " " << ml << std::endl;
 495      // std::cout << sh << " " << sh*sh << " " << msh2*msh2 << " " << msh << std::endl;
 496      //std::cout <<  (  12.0*c7c9*f1*G3 + 4.0*c7c7*f2*G2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
 497      //       <<" " << c9c9plusc10c10*f3*G1
 498      //       << " "<< 6.0*mlh*mlh*c9c9minusc10c10*f4
 499      //       << " "<< Gc << std::endl;
 500      //std::cout << C2 << " " << C1 << " "<< lambda_2 << " " << mc <<  " " << real(F*(conj(c9eff)*(2.0+sh)+conj(c7eff)*(1.0 + 6.0*sh - sh*sh)/sh)) << " " << sh << " " << r << std::endl;
 501      //std::cout << c9eff << " " << eta9 << " " <<c9eff0 << " " <<  c9eff1 << " " << alphas << " " << omega9 << " " << sh << std::endl;
 502
 503      //}
 504 //  else{
 505 //    if ( sh > 1.0) std::cout << "not a nan \n";
 506 //  }
 507   return prob;
 508 }
 509
 510 double EvtBtoXsllUtil::dGdsdupProb(double mb, double ms, double ml,
 511                                    double s,  double u)
 512 {
 513   // Compute the decay probability density function given a value of s and u
 514   // according to Ali-Hiller-Handoko-Morozumi's 1997 paper
 515   // see Appendix E
 516
 517   bool btod = false;
 518   bool nnlo = true;
 519
 520   double prob;
 521   double f1sp, f2sp, f3sp;
 522   double mbeff = 4.8;
 523
 524   //  double sh = s / (mb*mb);
 525   double sh  = s  / (mbeff*mbeff);
 526
 527   // if sh >1.0 code will return a nan. so just skip it
 528   if ( sh > 1.0 ) return 0.0;
 529
 530   EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo);
 531   EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo);
 532   EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod);
 533   EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod);
 534   EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo);
 535
 536   double alphas = 0.119/
 537      (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
 538
 539   double omega7 = -8.0/3.0*log(4.8/mb)
 540                   -4.0/3.0*EvtDiLog::DiLog(sh)
 541                   -2.0/9.0*EvtConst::pi*EvtConst::pi
 542                   -2.0/3.0*log(sh)*log(1.0-sh)
 543                   -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
 544     -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
 545     -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
 546   double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
 547
 548   double omega79 = -4.0/3.0*log(4.8/mb)
 549                    -4.0/3.0*EvtDiLog::DiLog(sh)
 550                    -2.0/9.0*EvtConst::pi*EvtConst::pi
 551                    -2.0/3.0*log(sh)*log(1.0-sh)
 552                    -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
 553                    -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
 554                    +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
 555   double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
 556
 557   double omega9 = - 2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
 558                  - 2.0/3.0*log(sh)*log(1.0-sh)
 559                  - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
 560                  - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
 561                  /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
 562                  + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
 563   double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
 564
 565   EvtComplex c7eff = eta7*c7eff0 + c7eff1;
 566   EvtComplex c9eff = eta9*c9eff0 + c9eff1;
 567   c10eff *= eta9;
 568
 569   double c7c7  = abs2(c7eff);
 570   double c7c9  = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1));
 571   double c7c10 = real((eta79*c7eff0 + c7eff1)*conj(eta9*c10eff));
 572   double c9c10 = real((eta9*c9eff0  + c9eff1)*conj(eta9*c10eff));
 573   double c9c9plusc10c10  = abs2(c9eff) + abs2(c10eff);
 574
 575   f1sp = ( pow(mb*mb-ms*ms,2) - s*s) * c9c9plusc10c10
 576          + 4.0*( pow(mb,4) - ms*ms*mb*mb - pow(ms,4)*(1.0 - ms*ms/(mb*mb))
 577          - 8.0*s*ms*ms - s*s*(1.0 + ms*ms/(mb*mb) ))*mb*mb*c7c7/s
 578     // kludged mass term
 579          *(1.0 + 2.0*ml*ml/s)
 580          - 8.0*(s*(mb*mb + ms*ms) - pow(mb*mb-ms*ms,2)) * c7c9
 581     // kludged mass term
 582          *(1.0 + 2.0*ml*ml/s);
 583
 584   f2sp = 4.0*s*c9c10 + 8.0*(mb*mb + ms*ms)*c7c10;
 585   f3sp = - (c9c9plusc10c10)
 586          + 4.0*(1.0 + pow(ms/mb,4)) * mb*mb*c7c7/s
 587     // kludged mass term
 588          *(1.0 + 2.0*ml*ml/s);
 589
 590   prob = (f1sp + f2sp*u + f3sp*u*u)/ pow(mb,3);
 591   if ( prob < 0.0 ) prob=0.;
 592
 593   return prob;
 594 }
 595
 596 double EvtBtoXsllUtil::FermiMomentum(double pf)
 597 {
 598   // Pick a value for the b-quark Fermi motion momentum
 599   // according to Ali's Gaussian model
 600
 601   double pb, pbmax, xbox, ybox;
 602   pb    = 0.0;
 603   pbmax = 5.0 * pf;
 604
 605   while (pb == 0.0)
 606   {
 607     xbox = EvtRandom::Flat(pbmax);
 608     ybox = EvtRandom::Flat();
 609     if (ybox < FermiMomentumProb(xbox, pf)) { pb = xbox;}
 610   }
 611
 612   return pb;
 613 }
 614
 615 double EvtBtoXsllUtil::FermiMomentumProb(double pb, double pf)
 616 {
 617   // Compute probability according to Ali's Gaussian model
 618   // the function chosen has a convenient maximum value of 1 for pb = pf
 619
 620   double prsq = (pb*pb)/(pf*pf);
 621   double prob = prsq * exp(1.0 - prsq);
 622
 623   return prob;
 624 }
 625