TEvtGen/EvtGenModels/EvtbTosllAmp.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 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtbTosllAmp.cc
  12 //
  13 // Description: Routine to implement semileptonic decays to pseudo-scalar
  14 //              mesons.
  15 //
  16 // Modification history:
  17 //
  18 //    DJL       April 17,1998       Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include "EvtGenBase/EvtPatches.hh"
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtGenKine.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenBase/EvtVector4C.hh"
  29 #include "EvtGenBase/EvtTensor4C.hh"
  30 #include "EvtGenBase/EvtDiracSpinor.hh"
  31 #include "EvtGenModels/EvtbTosllAmp.hh"
  32 #include "EvtGenBase/EvtId.hh"
  33 #include "EvtGenBase/EvtAmp.hh"
  34 #include "EvtGenBase/EvtScalarParticle.hh"
  35 #include "EvtGenBase/EvtVectorParticle.hh"
  36 #include "EvtGenBase/EvtDiLog.hh"
  37
  38 double EvtbTosllAmp::CalcMaxProb( EvtId parent, EvtId meson,
  39                                   EvtId lepton1, EvtId lepton2,
  40                                   EvtbTosllFF *FormFactors,
  41                                   double& poleSize) {
  42
  43   //This routine takes the arguements parent, meson, and lepton
  44   //number, and a form factor model, and returns a maximum
  45   //probability for this semileptonic form factor model.  A
  46   //brute force method is used.  The 2D cos theta lepton and
  47   //q2 phase space is probed.
  48
  49   //Start by declaring a particle at rest.
  50
  51   //It only makes sense to have a scalar parent.  For now.
  52   //This should be generalized later.
  53
  54   EvtScalarParticle *scalar_part;
  55   EvtParticle *root_part;
  56
  57   scalar_part=new EvtScalarParticle;
  58
  59   //cludge to avoid generating random numbers!
  60   scalar_part->noLifeTime();
  61
  62   EvtVector4R p_init;
  63
  64   p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
  65   scalar_part->init(parent,p_init);
  66   root_part=(EvtParticle *)scalar_part;
  67   root_part->setDiagonalSpinDensity();
  68
  69   EvtParticle *daughter, *lep1, *lep2;
  70
  71   EvtAmp amp;
  72
  73   EvtId listdaug[3];
  74   listdaug[0] = meson;
  75   listdaug[1] = lepton1;
  76   listdaug[2] = lepton2;
  77
  78   amp.init(parent,3,listdaug);
  79
  80   root_part->makeDaughters(3,listdaug);
  81   daughter=root_part->getDaug(0);
  82   lep1=root_part->getDaug(1);
  83   lep2=root_part->getDaug(2);
  84
  85   //cludge to avoid generating random numbers!
  86   daughter->noLifeTime();
  87   lep1->noLifeTime();
  88   lep2->noLifeTime();
  89
  90
  91   //Initial particle is unpolarized, well it is a scalar so it is
  92   //trivial
  93   EvtSpinDensity rho;
  94   rho.setDiag(root_part->getSpinStates());
  95
  96   double mass[3];
  97
  98   double m = root_part->mass();
  99
 100   EvtVector4R p4meson, p4lepton1, p4lepton2, p4w;
 101   double q2max;
 102
 103   double q2, elepton, plepton;
 104   int i,j;
 105   double erho,prho,costl;
 106
 107   double maxfoundprob = 0.0;
 108   double prob = -10.0;
 109   int massiter;
 110
 111   double maxpole=0;
 112
 113   for (massiter=0;massiter<3;massiter++){
 114
 115     mass[0] = EvtPDL::getMeanMass(meson);
 116     mass[1] = EvtPDL::getMeanMass(lepton1);
 117     mass[2] = EvtPDL::getMeanMass(lepton2);
 118     if ( massiter==1 ) {
 119       mass[0] = EvtPDL::getMinMass(meson);
 120     }
 121     if ( massiter==2 ) {
 122       mass[0] = EvtPDL::getMaxMass(meson);
 123       if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001;
 124     }
 125
 126     q2max = (m-mass[0])*(m-mass[0]);
 127
 128     //loop over q2
 129     //cout << "m " << m << "mass[0] " << mass[0] << " q2max "<< q2max << endl;
 130     for (i=0;i<25;i++) {
 131       //want to avoid picking up the tail of the photon propagator
 132       q2 = ((i+1.5)*q2max)/26.0;
 133
 134       if (i==0) q2=4*(mass[1]*mass[1]);
 135
 136       erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
 137
 138       prho = sqrt(erho*erho-mass[0]*mass[0]);
 139
 140       p4meson.set(erho,0.0,0.0,-1.0*prho);
 141       p4w.set(m-erho,0.0,0.0,prho);
 142
 143       //This is in the W rest frame
 144       elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
 145       plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
 146
 147       double probctl[3];
 148
 149       for (j=0;j<3;j++) {
 150
 151         costl = 0.99*(j - 1.0);
 152
 153         //These are in the W rest frame. Need to boost out into
 154         //the B frame.
 155         p4lepton1.set(elepton,0.0,
 156                   plepton*sqrt(1.0-costl*costl),plepton*costl);
 157         p4lepton2.set(elepton,0.0,
 158                  -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
 159
 160         EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
 161         p4lepton1=boostTo(p4lepton1,boost);
 162         p4lepton2=boostTo(p4lepton2,boost);
 163
 164         //Now initialize the daughters...
 165
 166         daughter->init(meson,p4meson);
 167         lep1->init(lepton1,p4lepton1);
 168         lep2->init(lepton2,p4lepton2);
 169
 170         CalcAmp(root_part,amp,FormFactors);
 171
 172         //Now find the probability at this q2 and cos theta lepton point
 173         //and compare to maxfoundprob.
 174
 175         //Do a little magic to get the probability!!
 176
 177         //cout <<"amp:"<<amp.getSpinDensity()<<endl;
 178
 179         prob = rho.normalizedProb(amp.getSpinDensity());
 180
 181         //cout << "prob:"<<q2<<" "<<costl<<" "<<prob<<endl;
 182
 183         probctl[j]=prob;
 184       }
 185
 186       //probclt contains prob at ctl=-1,0,1.
 187       //prob=a+b*ctl+c*ctl^2
 188
 189       double a=probctl[1];
 190       double b=0.5*(probctl[2]-probctl[0]);
 191       double c=0.5*(probctl[2]+probctl[0])-probctl[1];
 192
 193       prob=probctl[0];
 194       if (probctl[1]>prob) prob=probctl[1];
 195       if (probctl[2]>prob) prob=probctl[2];
 196
 197       if (fabs(c)>1e-20){
 198         double ctlx=-0.5*b/c;
 199         if (fabs(ctlx)<1.0){
 200           double probtmp=a+b*ctlx+c*ctlx*ctlx;
 201           if (probtmp>prob) prob=probtmp;
 202         }
 203
 204       }
 205
 206       //report(DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
 207       //                            << probctl[0]<<" "
 208       //                            << probctl[1]<<" "
 209       //                            << probctl[2]<<endl;
 210
 211       if (i==0) {
 212         maxpole=prob;
 213         continue;
 214       }
 215
 216       if ( prob > maxfoundprob ) {
 217         maxfoundprob = prob;
 218       }
 219
 220       //cout << "q2,maxfoundprob:"<<q2<<" "<<maxfoundprob<<endl;
 221
 222     }
 223     if ( EvtPDL::getWidth(meson) <= 0.0 ) {
 224       //if the particle is narrow dont bother with changing the mass.
 225       massiter = 4;
 226     }
 227
 228   }
 229
 230   root_part->deleteTree();
 231
 232   poleSize=0.04*(maxpole/maxfoundprob)*4*(mass[1]*mass[1]);
 233
 234   //poleSize=0.002;
 235
 236   //cout <<"maxfoundprob,maxpole,poleSize:"<<maxfoundprob<<" "
 237   //     <<maxpole<<" "<<poleSize<<endl;
 238
 239   maxfoundprob *=1.15;
 240
 241   return maxfoundprob;
 242
 243 }
 244
 245
 246 EvtComplex EvtbTosllAmp::GetC7Eff(double q2, bool nnlo)
 247 {
 248
 249   if (!nnlo) return -0.313;
 250   double mbeff = 4.8;
 251   double shat = q2/mbeff/mbeff;
 252   double logshat;
 253   logshat = log(shat);
 254
 255   double muscale;
 256   muscale = 2.5;
 257   double alphas;
 258   alphas = 0.267;
 259   double A7;
 260   A7 = -0.353 + 0.023;
 261   double A8;
 262   A8 = -0.164;
 263   double A9;
 264   A9 = 4.287 + (-0.218);
 265   double A10;
 266   A10 = -4.592 + 0.379;
 267   double C1;
 268   C1 = -0.697;
 269   double C2;
 270   C2 = 1.046;
 271   double T9;
 272   T9 = 0.114 + 0.280;
 273   double U9;
 274   U9 = 0.045 + 0.023;
 275   double W9;
 276   W9 = 0.044 + 0.016;
 277
 278   double Lmu;
 279   Lmu = log(muscale/mbeff);
 280
 281   EvtComplex uniti(0.0,1.0);
 282
 283   EvtComplex c7eff;
 284   if (shat > 0.25)
 285   {
 286    c7eff = A7;
 287    return c7eff;
 288   }
 289
 290
 291
 292
 293   // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
 294   muscale = 5.0;
 295   alphas = 0.215;
 296   A7 = -0.312 + 0.008;
 297   A8 = -0.148;
 298   A9 = 4.174 + (-0.035);
 299   A10 = -4.592 + 0.379;
 300   C1 = -0.487;
 301   C2 = 1.024;
 302   T9 = 0.374 + 0.252;
 303   U9 = 0.033 + 0.015;
 304   W9 = 0.032 + 0.012;
 305   Lmu = log(muscale/mbeff);
 306
 307   EvtComplex F71;
 308   EvtComplex f71;
 309   EvtComplex k7100(-0.68192,-0.074998);
 310   EvtComplex k7101(0.0,0.0);
 311   EvtComplex k7110(-0.23935,-0.12289);
 312   EvtComplex k7111(0.0027424,0.019676);
 313   EvtComplex k7120(-0.0018555,-0.175);
 314   EvtComplex k7121(0.022864,0.011456);
 315   EvtComplex k7130(0.28248,-0.12783);
 316   EvtComplex k7131(0.029027,-0.0082265);
 317   f71 = k7100 + k7101*logshat + shat*(k7110 + k7111*logshat) +
 318         shat*shat*(k7120 + k7121*logshat) +
 319         shat*shat*shat*(k7130 + k7131*logshat);
 320   F71 = (-208.0/243.0)*Lmu + f71;
 321
 322   EvtComplex F72;
 323   EvtComplex f72;
 324   EvtComplex k7200(4.0915,0.44999);
 325   EvtComplex k7201(0.0,0.0);
 326   EvtComplex k7210(1.4361,0.73732);
 327   EvtComplex k7211(-0.016454,-0.11806);
 328   EvtComplex k7220(0.011133,1.05);
 329   EvtComplex k7221(-0.13718,-0.068733);
 330   EvtComplex k7230(-1.6949,0.76698);
 331   EvtComplex k7231(-0.17416,0.049359);
 332   f72 = k7200 + k7201*logshat + shat*(k7210 + k7211*logshat) +
 333         shat*shat*(k7220 + k7221*logshat) +
 334         shat*shat*shat*(k7230 + k7231*logshat);
 335   F72 = (416.0/81.0)*Lmu + f72;
 336
 337   EvtComplex F78;
 338   F78 = (-32.0/9.0)*Lmu + 8.0*EvtConst::pi*EvtConst::pi/27.0 + (-44.0/9.0)
 339         + (-8.0*EvtConst::pi/9.0)*uniti +
 340         (4.0/3.0*EvtConst::pi*EvtConst::pi - 40.0/3.0)*shat +
 341         (32.0*EvtConst::pi*EvtConst::pi/9.0 - 316.0/9.0)*shat*shat +
 342         (200.0*EvtConst::pi*EvtConst::pi/27.0 - 658.0/9.0)*shat*shat*shat +
 343     (-8.0*logshat/9.0)*(shat + shat*shat + shat*shat*shat);
 344
 345   c7eff = A7 - alphas/(4.0*EvtConst::pi)*(C1*F71 + C2*F72 + A8*F78);
 346
 347   return c7eff;
 348 }
 349
 350
 351 EvtComplex EvtbTosllAmp::GetC9Eff(double q2, bool nnlo, bool btod)
 352 {
 353
 354   if (!nnlo) return 4.344;
 355   double mbeff = 4.8;
 356   double shat = q2/mbeff/mbeff;
 357   double logshat;
 358   logshat = log(shat);
 359   double mchat = 0.29;
 360
 361
 362   double muscale;
 363   muscale = 2.5;
 364   double alphas;
 365   alphas = 0.267;
 366   double A7;
 367   A7 = -0.353 + 0.023;
 368   double A8;
 369   A8 = -0.164;
 370   double A9;
 371   A9 = 4.287 + (-0.218);
 372   double A10;
 373   A10 = -4.592 + 0.379;
 374   double C1;
 375   C1 = -0.697;
 376   double C2;
 377   C2 = 1.046;
 378   double T9;
 379   T9 = 0.114 + 0.280;
 380   double U9;
 381   U9 = 0.045 + 0.023;
 382   double W9;
 383   W9 = 0.044 + 0.016;
 384
 385   double Lmu;
 386   Lmu = log(muscale/mbeff);
 387
 388
 389   EvtComplex uniti(0.0,1.0);
 390
 391   EvtComplex hc;
 392   double xarg;
 393   xarg = 4.0*mchat/shat;
 394   hc = -4.0/9.0*log(mchat*mchat) + 8.0/27.0 + 4.0*xarg/9.0;
 395
 396 if (xarg < 1.0)
 397   {
 398     hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 399       (log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
 400        uniti*EvtConst::pi);
 401   }
 402   else
 403   {
 404     hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 405       2.0*atan(1.0/sqrt(xarg - 1.0));
 406   }
 407
 408   EvtComplex h1;
 409   xarg = 4.0/shat;
 410   h1 = 8.0/27.0 + 4.0*xarg/9.0;
 411   if (xarg < 1.0)
 412   {
 413     h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 414       (log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
 415        uniti*EvtConst::pi);
 416   }
 417   else
 418   {
 419     h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 420       2.0*atan(1.0/sqrt(xarg - 1.0));
 421   }
 422
 423
 424   EvtComplex h0;
 425   h0 = 8.0/27.0 - 4.0*log(2.0)/9.0 + 4.0*uniti*EvtConst::pi/9.0;
 426
 427
 428   // X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
 429   // h(\hat m_u^2, hat s))
 430   EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0);
 431   EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0);
 432   EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0);
 433   EvtComplex Vtb(1.0,0.0);
 434
 435   EvtComplex Xd;
 436   Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
 437
 438
 439   EvtComplex c9eff=4.344;
 440   if (shat > 0.25)
 441   {
 442    c9eff =  A9 + T9*hc + U9*h1 + W9*h0;
 443    if (btod)
 444    {
 445     c9eff += Xd;
 446    }
 447
 448    return c9eff;
 449   }
 450
 451   // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
 452   muscale = 5.0;
 453   alphas = 0.215;
 454   A9 = 4.174 + (-0.035);
 455   C1 = -0.487;
 456   C2 = 1.024;
 457   A8 = -0.148;
 458   T9 = 0.374 + 0.252;
 459   U9 = 0.033 + 0.015;
 460   W9 = 0.032 + 0.012;
 461   Lmu = log(muscale/mbeff);
 462
 463   EvtComplex F91;
 464   EvtComplex f91;
 465   EvtComplex k9100(-11.973,0.16371);
 466   EvtComplex k9101(-0.081271,-0.059691);
 467   EvtComplex k9110(-28.432,-0.25044);
 468   EvtComplex k9111(-0.040243,0.016442);
 469   EvtComplex k9120(-57.114,-0.86486);
 470   EvtComplex k9121(-0.035191,0.027909);
 471   EvtComplex k9130(-128.8,-2.5243);
 472   EvtComplex k9131(-0.017587,0.050639);
 473   f91 = k9100 + k9101*logshat + shat*(k9110 + k9111*logshat) +
 474         shat*shat*(k9120 + k9121*logshat) +
 475         shat*shat*shat*(k9130 + k9131*logshat);
 476   F91 = (-1424.0/729.0 + 16.0*uniti*EvtConst::pi/243.0
 477          + 64.0/27.0*log(mchat))*Lmu - 16.0*Lmu*logshat/243.0 +
 478         (16.0/1215.0 - 32.0/135.0/mchat/mchat)*Lmu*shat +
 479         (4.0/2835.0 - 8.0/315.0/mchat/mchat/mchat/mchat)*Lmu*shat*shat +
 480     (16.0/76545.0 - 32.0/8505.0/mchat/mchat/mchat/mchat/mchat/mchat)*
 481     Lmu*shat*shat*shat -256.0*Lmu*Lmu/243.0 + f91;
 482
 483   EvtComplex F92;
 484   EvtComplex f92;
 485   EvtComplex k9200(6.6338,-0.98225);
 486   EvtComplex k9201(0.48763,0.35815);
 487   EvtComplex k9210(3.3585,1.5026);
 488   EvtComplex k9211(0.24146,-0.098649);
 489   EvtComplex k9220(-1.1906,5.1892);
 490   EvtComplex k9221(0.21115,-0.16745);
 491   EvtComplex k9230(-17.12,15.146);
 492   EvtComplex k9231(0.10552,-0.30383);
 493   f92 = k9200 + k9201*logshat + shat*(k9210 + k9211*logshat) +
 494         shat*shat*(k9220 + k9221*logshat) +
 495         shat*shat*shat*(k9230 + k9231*logshat);
 496   F92 = (256.0/243.0 - 32.0*uniti*EvtConst::pi/81.0
 497          - 128.0/9.0*log(mchat))*Lmu + 32.0*Lmu*logshat/81.0 +
 498         (-32.0/405.0 + 64.0/45.0/mchat/mchat)*Lmu*shat +
 499         (-8.0/945.0 + 16.0/105.0/mchat/mchat/mchat/mchat)*Lmu*shat*shat +
 500     (-32.0/25515.0 + 64.0/2835.0/mchat/mchat/mchat/mchat/mchat/mchat)*
 501     Lmu*shat*shat*shat + 512.0*Lmu*Lmu/81.0 + f92;
 502
 503   EvtComplex F98;
 504   F98 = 104.0/9.0 - 32.0*EvtConst::pi*EvtConst::pi/27.0 +
 505         (1184.0/27.0 - 40.0*EvtConst::pi*EvtConst::pi/9.0)*shat +
 506         (14212.0/135.0 - 32.0*EvtConst::pi*EvtConst::pi/3.0)*shat*shat +
 507     (193444.0/945.0 - 560.0*EvtConst::pi*EvtConst::pi/27.0)*shat*shat*shat +
 508         16.0*logshat/9.0*(1.0 + shat + shat*shat + shat*shat*shat);
 509
 510   Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
 511
 512   c9eff = A9 + T9*hc + U9*h1 + W9*h0 -
 513     alphas/(4.0*EvtConst::pi)*(C1*F91 + C2*F92 + A8*F98);
 514   if (btod)
 515   {
 516    c9eff += Xd;
 517   }
 518
 519   return c9eff;
 520 }
 521
 522 EvtComplex EvtbTosllAmp::GetC10Eff(double /*q2*/, bool nnlo)
 523 {
 524
 525   if (!nnlo) return -4.669;
 526   double A10;
 527   A10 = -4.592 + 0.379;
 528
 529   EvtComplex c10eff;
 530   c10eff = A10;
 531
 532   return c10eff;
 533 }
 534
 535 double EvtbTosllAmp::dGdsProb(double mb, double ms, double ml,
 536                                 double s)
 537 {
 538   // Compute the decay probability density function given a value of s
 539   // according to Ali's paper
 540
 541
 542   double delta, lambda, prob;
 543   double f1, f2, f3, f4;
 544   double msh, mlh, sh;
 545
 546   mlh = ml / mb;
 547   msh = ms / mb;
 548   sh  = s  / (mb*mb);
 549
 550   EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
 551   EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
 552   EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);
 553
 554   double alphas = 0.119/
 555      (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
 556   double omega9 = -2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
 557                  - 2.0/3.0*log(sh)*log(1.0-sh)
 558                  - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
 559                  - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
 560                  /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
 561                  + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
 562   double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
 563   double omega7 = -8.0/3.0*log(4.8/mb)
 564                   -4.0/3.0*EvtDiLog::DiLog(sh)
 565                   -2.0/9.0*EvtConst::pi*EvtConst::pi
 566                   -2.0/3.0*log(sh)*log(1.0-sh)
 567                   -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
 568     -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
 569     -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
 570   double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
 571
 572   double omega79 = -4.0/3.0*log(4.8/mb)
 573                    -4.0/3.0*EvtDiLog::DiLog(sh)
 574                    -2.0/9.0*EvtConst::pi*EvtConst::pi
 575                    -2.0/3.0*log(sh)*log(1.0-sh)
 576                    -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
 577                    -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
 578                    +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
 579   double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
 580
 581   double c7c9 = abs(c7eff)*real(c9eff);
 582   c7c9 *= pow(eta79,2);
 583   double c7c7 = pow(abs(c7eff),2);
 584   c7c7 *= pow(eta7,2);
 585
 586   double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
 587   c9c9plusc10c10 *= pow(eta9,2);
 588   double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
 589   c9c9minusc10c10 *= pow(eta9,2);
 590
 591   lambda = 1.0 + sh*sh + pow(msh,4) - 2.0*(sh + sh*msh*msh + msh*msh);
 592
 593   f1 = pow(1.0-msh*msh,2) - sh*(1.0 + msh*msh);
 594   f2 = 2.0*(1.0 + msh*msh) * pow(1.0-msh*msh,2)
 595        - sh*(1.0 + 14.0*msh*msh + pow(msh,4)) - sh*sh*(1.0 + msh*msh);
 596   f3 = pow(1.0-msh*msh,2) + sh*(1.0 + msh*msh) - 2.0*sh*sh
 597        + lambda*2.0*mlh*mlh/sh;
 598   f4 = 1.0 - sh + msh*msh;
 599
 600   delta = (  12.0*c7c9*f1 + 4.0*c7c7*f2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
 601             + c9c9plusc10c10*f3
 602             + 6.0*mlh*mlh*c9c9minusc10c10*f4;
 603
 604   prob =  sqrt(lambda*(1.0 - 4.0*mlh*mlh/sh)) * delta;
 605
 606   return prob;
 607 }
 608
 609 double EvtbTosllAmp::dGdsdupProb(double mb, double ms, double ml,
 610                                    double s,  double u)
 611 {
 612   // Compute the decay probability density function given a value of s and u
 613   // according to Ali's paper
 614
 615   double prob;
 616   double f1sp, f2sp, f3sp;
 617
 618   double sh = s / (mb*mb);
 619
 620   EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
 621   EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
 622   EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);
 623
 624   double alphas = 0.119/
 625      (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
 626   double omega9 = - 2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
 627                  - 2.0/3.0*log(sh)*log(1.0-sh)
 628                  - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
 629                  - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
 630                  /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
 631                  + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
 632   double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
 633   double omega7 = -8.0/3.0*log(4.8/mb)
 634                   -4.0/3.0*EvtDiLog::DiLog(sh)
 635                   -2.0/9.0*EvtConst::pi*EvtConst::pi
 636                   -2.0/3.0*log(sh)*log(1.0-sh)
 637                   -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
 638     -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
 639     -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
 640   double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
 641
 642   double omega79 = -4.0/3.0*log(4.8/mb)
 643                    -4.0/3.0*EvtDiLog::DiLog(sh)
 644                    -2.0/9.0*EvtConst::pi*EvtConst::pi
 645                    -2.0/3.0*log(sh)*log(1.0-sh)
 646                    -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
 647                    -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
 648                    +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
 649   double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
 650
 651   double c7c9 = abs(c7eff)*real(c9eff);
 652   c7c9 *= pow(eta79,2);
 653   double c7c7 = pow(abs(c7eff),2);
 654   c7c7 *= pow(eta7,2);
 655
 656   double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
 657   c9c9plusc10c10 *= pow(eta9,2);
 658   double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
 659   c9c9minusc10c10 *= pow(eta9,2);
 660   double c7c10 = abs(c7eff)*real(c10eff);
 661   c7c10 *= eta7; c7c10 *= eta9;
 662   double c9c10 = real(c9eff)*real(c10eff);
 663   c9c10 *= pow(eta9,2);
 664
 665   f1sp = ( pow(mb*mb-ms*ms,2) - s*s) * c9c9plusc10c10
 666          + 4.0*( pow(mb,4) - ms*ms*mb*mb - pow(ms,4)*(1.0 - ms*ms/(mb*mb))
 667          - 8.0*s*ms*ms - s*s*(1.0 + ms*ms/(mb*mb) ))*mb*mb*c7c7/s
 668     // kludged mass term
 669          *(1.0 + 2.0*ml*ml/s)
 670          - 8.0*(s*(mb*mb + ms*ms) - pow(mb*mb-ms*ms,2)) * c7c9
 671     // kludged mass term
 672          *(1.0 + 2.0*ml*ml/s);
 673
 674   f2sp = 4.0*s*c9c10 + 8.0*(mb*mb + ms*ms)*c7c10;
 675   f3sp = - (c9c9plusc10c10)
 676          + 4.0*(1.0 + pow(ms/mb,4)) * mb*mb*c7c7/s
 677     // kludged mass term
 678          *(1.0 + 2.0*ml*ml/s);
 679
 680   prob = (f1sp + f2sp*u + f3sp*u*u)/ pow(mb,3);
 681
 682   return prob;
 683 }
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696