TEvtGen/EvtGenModels/EvtbTosllAmp.cpp

   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 C1;
 264   C1 = -0.697;
 265   double C2;
 266   C2 = 1.046;
 267
 268   double Lmu;
 269   Lmu = log(muscale/mbeff);
 270
 271   EvtComplex uniti(0.0,1.0);
 272
 273   EvtComplex c7eff;
 274   if (shat > 0.25)
 275   {
 276    c7eff = A7;
 277    return c7eff;
 278   }
 279
 280   // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
 281   muscale = 5.0;
 282   alphas = 0.215;
 283   A7 = -0.312 + 0.008;
 284   A8 = -0.148;
 285   C1 = -0.487;
 286   C2 = 1.024;
 287   Lmu = log(muscale/mbeff);
 288
 289   EvtComplex F71;
 290   EvtComplex f71;
 291   EvtComplex k7100(-0.68192,-0.074998);
 292   EvtComplex k7101(0.0,0.0);
 293   EvtComplex k7110(-0.23935,-0.12289);
 294   EvtComplex k7111(0.0027424,0.019676);
 295   EvtComplex k7120(-0.0018555,-0.175);
 296   EvtComplex k7121(0.022864,0.011456);
 297   EvtComplex k7130(0.28248,-0.12783);
 298   EvtComplex k7131(0.029027,-0.0082265);
 299   f71 = k7100 + k7101*logshat + shat*(k7110 + k7111*logshat) +
 300         shat*shat*(k7120 + k7121*logshat) +
 301         shat*shat*shat*(k7130 + k7131*logshat);
 302   F71 = (-208.0/243.0)*Lmu + f71;
 303
 304   EvtComplex F72;
 305   EvtComplex f72;
 306   EvtComplex k7200(4.0915,0.44999);
 307   EvtComplex k7201(0.0,0.0);
 308   EvtComplex k7210(1.4361,0.73732);
 309   EvtComplex k7211(-0.016454,-0.11806);
 310   EvtComplex k7220(0.011133,1.05);
 311   EvtComplex k7221(-0.13718,-0.068733);
 312   EvtComplex k7230(-1.6949,0.76698);
 313   EvtComplex k7231(-0.17416,0.049359);
 314   f72 = k7200 + k7201*logshat + shat*(k7210 + k7211*logshat) +
 315         shat*shat*(k7220 + k7221*logshat) +
 316         shat*shat*shat*(k7230 + k7231*logshat);
 317   F72 = (416.0/81.0)*Lmu + f72;
 318
 319   EvtComplex F78;
 320   F78 = (-32.0/9.0)*Lmu + 8.0*EvtConst::pi*EvtConst::pi/27.0 + (-44.0/9.0)
 321         + (-8.0*EvtConst::pi/9.0)*uniti +
 322         (4.0/3.0*EvtConst::pi*EvtConst::pi - 40.0/3.0)*shat +
 323         (32.0*EvtConst::pi*EvtConst::pi/9.0 - 316.0/9.0)*shat*shat +
 324         (200.0*EvtConst::pi*EvtConst::pi/27.0 - 658.0/9.0)*shat*shat*shat +
 325     (-8.0*logshat/9.0)*(shat + shat*shat + shat*shat*shat);
 326
 327   c7eff = A7 - alphas/(4.0*EvtConst::pi)*(C1*F71 + C2*F72 + A8*F78);
 328
 329   return c7eff;
 330 }
 331
 332
 333 EvtComplex EvtbTosllAmp::GetC9Eff(double q2, bool nnlo, bool btod)
 334 {
 335
 336   if (!nnlo) return 4.344;
 337   double mbeff = 4.8;
 338   double shat = q2/mbeff/mbeff;
 339   double logshat;
 340   logshat = log(shat);
 341   double mchat = 0.29;
 342
 343
 344   double muscale;
 345   muscale = 2.5;
 346   double alphas;
 347   alphas = 0.267;
 348   double A8;
 349   A8 = -0.164;
 350   double A9;
 351   A9 = 4.287 + (-0.218);
 352   double C1;
 353   C1 = -0.697;
 354   double C2;
 355   C2 = 1.046;
 356   double T9;
 357   T9 = 0.114 + 0.280;
 358   double U9;
 359   U9 = 0.045 + 0.023;
 360   double W9;
 361   W9 = 0.044 + 0.016;
 362
 363   double Lmu;
 364   Lmu = log(muscale/mbeff);
 365
 366
 367   EvtComplex uniti(0.0,1.0);
 368
 369   EvtComplex hc;
 370   double xarg;
 371   xarg = 4.0*mchat/shat;
 372   hc = -4.0/9.0*log(mchat*mchat) + 8.0/27.0 + 4.0*xarg/9.0;
 373
 374 if (xarg < 1.0)
 375   {
 376     hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 377       (log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
 378        uniti*EvtConst::pi);
 379   }
 380   else
 381   {
 382     hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 383       2.0*atan(1.0/sqrt(xarg - 1.0));
 384   }
 385
 386   EvtComplex h1;
 387   xarg = 4.0/shat;
 388   h1 = 8.0/27.0 + 4.0*xarg/9.0;
 389   if (xarg < 1.0)
 390   {
 391     h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 392       (log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
 393        uniti*EvtConst::pi);
 394   }
 395   else
 396   {
 397     h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
 398       2.0*atan(1.0/sqrt(xarg - 1.0));
 399   }
 400
 401
 402   EvtComplex h0;
 403   h0 = 8.0/27.0 - 4.0*log(2.0)/9.0 + 4.0*uniti*EvtConst::pi/9.0;
 404
 405
 406   // X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
 407   // h(\hat m_u^2, hat s))
 408   EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0);
 409   EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0);
 410   EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0);
 411   EvtComplex Vtb(1.0,0.0);
 412
 413   EvtComplex Xd;
 414   Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
 415
 416
 417   EvtComplex c9eff=4.344;
 418   if (shat > 0.25)
 419   {
 420    c9eff =  A9 + T9*hc + U9*h1 + W9*h0;
 421    if (btod)
 422    {
 423     c9eff += Xd;
 424    }
 425
 426    return c9eff;
 427   }
 428
 429   // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
 430   muscale = 5.0;
 431   alphas = 0.215;
 432   A9 = 4.174 + (-0.035);
 433   C1 = -0.487;
 434   C2 = 1.024;
 435   A8 = -0.148;
 436   T9 = 0.374 + 0.252;
 437   U9 = 0.033 + 0.015;
 438   W9 = 0.032 + 0.012;
 439   Lmu = log(muscale/mbeff);
 440
 441   EvtComplex F91;
 442   EvtComplex f91;
 443   EvtComplex k9100(-11.973,0.16371);
 444   EvtComplex k9101(-0.081271,-0.059691);
 445   EvtComplex k9110(-28.432,-0.25044);
 446   EvtComplex k9111(-0.040243,0.016442);
 447   EvtComplex k9120(-57.114,-0.86486);
 448   EvtComplex k9121(-0.035191,0.027909);
 449   EvtComplex k9130(-128.8,-2.5243);
 450   EvtComplex k9131(-0.017587,0.050639);
 451   f91 = k9100 + k9101*logshat + shat*(k9110 + k9111*logshat) +
 452         shat*shat*(k9120 + k9121*logshat) +
 453         shat*shat*shat*(k9130 + k9131*logshat);
 454   F91 = (-1424.0/729.0 + 16.0*uniti*EvtConst::pi/243.0
 455          + 64.0/27.0*log(mchat))*Lmu - 16.0*Lmu*logshat/243.0 +
 456         (16.0/1215.0 - 32.0/135.0/mchat/mchat)*Lmu*shat +
 457         (4.0/2835.0 - 8.0/315.0/mchat/mchat/mchat/mchat)*Lmu*shat*shat +
 458     (16.0/76545.0 - 32.0/8505.0/mchat/mchat/mchat/mchat/mchat/mchat)*
 459     Lmu*shat*shat*shat -256.0*Lmu*Lmu/243.0 + f91;
 460
 461   EvtComplex F92;
 462   EvtComplex f92;
 463   EvtComplex k9200(6.6338,-0.98225);
 464   EvtComplex k9201(0.48763,0.35815);
 465   EvtComplex k9210(3.3585,1.5026);
 466   EvtComplex k9211(0.24146,-0.098649);
 467   EvtComplex k9220(-1.1906,5.1892);
 468   EvtComplex k9221(0.21115,-0.16745);
 469   EvtComplex k9230(-17.12,15.146);
 470   EvtComplex k9231(0.10552,-0.30383);
 471   f92 = k9200 + k9201*logshat + shat*(k9210 + k9211*logshat) +
 472         shat*shat*(k9220 + k9221*logshat) +
 473         shat*shat*shat*(k9230 + k9231*logshat);
 474   F92 = (256.0/243.0 - 32.0*uniti*EvtConst::pi/81.0
 475          - 128.0/9.0*log(mchat))*Lmu + 32.0*Lmu*logshat/81.0 +
 476         (-32.0/405.0 + 64.0/45.0/mchat/mchat)*Lmu*shat +
 477         (-8.0/945.0 + 16.0/105.0/mchat/mchat/mchat/mchat)*Lmu*shat*shat +
 478     (-32.0/25515.0 + 64.0/2835.0/mchat/mchat/mchat/mchat/mchat/mchat)*
 479     Lmu*shat*shat*shat + 512.0*Lmu*Lmu/81.0 + f92;
 480
 481   EvtComplex F98;
 482   F98 = 104.0/9.0 - 32.0*EvtConst::pi*EvtConst::pi/27.0 +
 483         (1184.0/27.0 - 40.0*EvtConst::pi*EvtConst::pi/9.0)*shat +
 484         (14212.0/135.0 - 32.0*EvtConst::pi*EvtConst::pi/3.0)*shat*shat +
 485     (193444.0/945.0 - 560.0*EvtConst::pi*EvtConst::pi/27.0)*shat*shat*shat +
 486         16.0*logshat/9.0*(1.0 + shat + shat*shat + shat*shat*shat);
 487
 488   Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
 489
 490   c9eff = A9 + T9*hc + U9*h1 + W9*h0 -
 491     alphas/(4.0*EvtConst::pi)*(C1*F91 + C2*F92 + A8*F98);
 492   if (btod)
 493   {
 494    c9eff += Xd;
 495   }
 496
 497   return c9eff;
 498 }
 499
 500 EvtComplex EvtbTosllAmp::GetC10Eff(double /*q2*/, bool nnlo)
 501 {
 502
 503   if (!nnlo) return -4.669;
 504   double A10;
 505   A10 = -4.592 + 0.379;
 506
 507   EvtComplex c10eff;
 508   c10eff = A10;
 509
 510   return c10eff;
 511 }
 512
 513 double EvtbTosllAmp::dGdsProb(double mb, double ms, double ml,
 514                                 double s)
 515 {
 516   // Compute the decay probability density function given a value of s
 517   // according to Ali's paper
 518
 519
 520   double delta, lambda, prob;
 521   double f1, f2, f3, f4;
 522   double msh, mlh, sh;
 523
 524   mlh = ml / mb;
 525   msh = ms / mb;
 526   sh  = s  / (mb*mb);
 527
 528   EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
 529   EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
 530   EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);
 531
 532   double alphas = 0.119/
 533      (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
 534   double omega9 = -2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
 535                  - 2.0/3.0*log(sh)*log(1.0-sh)
 536                  - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
 537                  - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
 538                  /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
 539                  + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
 540   double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
 541   double omega7 = -8.0/3.0*log(4.8/mb)
 542                   -4.0/3.0*EvtDiLog::DiLog(sh)
 543                   -2.0/9.0*EvtConst::pi*EvtConst::pi
 544                   -2.0/3.0*log(sh)*log(1.0-sh)
 545                   -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
 546     -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
 547     -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
 548   double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
 549
 550   double omega79 = -4.0/3.0*log(4.8/mb)
 551                    -4.0/3.0*EvtDiLog::DiLog(sh)
 552                    -2.0/9.0*EvtConst::pi*EvtConst::pi
 553                    -2.0/3.0*log(sh)*log(1.0-sh)
 554                    -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
 555                    -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
 556                    +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
 557   double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
 558
 559   double c7c9 = abs(c7eff)*real(c9eff);
 560   c7c9 *= pow(eta79,2);
 561   double c7c7 = pow(abs(c7eff),2);
 562   c7c7 *= pow(eta7,2);
 563
 564   double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
 565   c9c9plusc10c10 *= pow(eta9,2);
 566   double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
 567   c9c9minusc10c10 *= pow(eta9,2);
 568
 569   lambda = 1.0 + sh*sh + pow(msh,4) - 2.0*(sh + sh*msh*msh + msh*msh);
 570
 571   f1 = pow(1.0-msh*msh,2) - sh*(1.0 + msh*msh);
 572   f2 = 2.0*(1.0 + msh*msh) * pow(1.0-msh*msh,2)
 573        - sh*(1.0 + 14.0*msh*msh + pow(msh,4)) - sh*sh*(1.0 + msh*msh);
 574   f3 = pow(1.0-msh*msh,2) + sh*(1.0 + msh*msh) - 2.0*sh*sh
 575        + lambda*2.0*mlh*mlh/sh;
 576   f4 = 1.0 - sh + msh*msh;
 577
 578   delta = (  12.0*c7c9*f1 + 4.0*c7c7*f2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
 579             + c9c9plusc10c10*f3
 580             + 6.0*mlh*mlh*c9c9minusc10c10*f4;
 581
 582   prob =  sqrt(lambda*(1.0 - 4.0*mlh*mlh/sh)) * delta;
 583
 584   return prob;
 585 }
 586
 587 double EvtbTosllAmp::dGdsdupProb(double mb, double ms, double ml,
 588                                    double s,  double u)
 589 {
 590   // Compute the decay probability density function given a value of s and u
 591   // according to Ali's paper
 592
 593   double prob;
 594   double f1sp, f2sp, f3sp;
 595
 596   double sh = s / (mb*mb);
 597
 598   EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
 599   EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
 600   EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);
 601
 602   double alphas = 0.119/
 603      (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
 604   double omega9 = - 2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
 605                  - 2.0/3.0*log(sh)*log(1.0-sh)
 606                  - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
 607                  - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
 608                  /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
 609                  + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
 610   double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
 611   double omega7 = -8.0/3.0*log(4.8/mb)
 612                   -4.0/3.0*EvtDiLog::DiLog(sh)
 613                   -2.0/9.0*EvtConst::pi*EvtConst::pi
 614                   -2.0/3.0*log(sh)*log(1.0-sh)
 615                   -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
 616     -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
 617     -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
 618   double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
 619
 620   double omega79 = -4.0/3.0*log(4.8/mb)
 621                    -4.0/3.0*EvtDiLog::DiLog(sh)
 622                    -2.0/9.0*EvtConst::pi*EvtConst::pi
 623                    -2.0/3.0*log(sh)*log(1.0-sh)
 624                    -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
 625                    -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
 626                    +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
 627   double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
 628
 629   double c7c9 = abs(c7eff)*real(c9eff);
 630   c7c9 *= pow(eta79,2);
 631   double c7c7 = pow(abs(c7eff),2);
 632   c7c7 *= pow(eta7,2);
 633
 634   double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
 635   c9c9plusc10c10 *= pow(eta9,2);
 636   double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
 637   c9c9minusc10c10 *= pow(eta9,2);
 638   double c7c10 = abs(c7eff)*real(c10eff);
 639   c7c10 *= eta7; c7c10 *= eta9;
 640   double c9c10 = real(c9eff)*real(c10eff);
 641   c9c10 *= pow(eta9,2);
 642
 643   f1sp = ( pow(mb*mb-ms*ms,2) - s*s) * c9c9plusc10c10
 644          + 4.0*( pow(mb,4) - ms*ms*mb*mb - pow(ms,4)*(1.0 - ms*ms/(mb*mb))
 645          - 8.0*s*ms*ms - s*s*(1.0 + ms*ms/(mb*mb) ))*mb*mb*c7c7/s
 646     // kludged mass term
 647          *(1.0 + 2.0*ml*ml/s)
 648          - 8.0*(s*(mb*mb + ms*ms) - pow(mb*mb-ms*ms,2)) * c7c9
 649     // kludged mass term
 650          *(1.0 + 2.0*ml*ml/s);
 651
 652   f2sp = 4.0*s*c9c10 + 8.0*(mb*mb + ms*ms)*c7c10;
 653   f3sp = - (c9c9plusc10c10)
 654          + 4.0*(1.0 + pow(ms/mb,4)) * mb*mb*c7c7/s
 655     // kludged mass term
 656          *(1.0 + 2.0*ml*ml/s);
 657
 658   prob = (f1sp + f2sp*u + f3sp*u*u)/ pow(mb,3);
 659
 660   return prob;
 661 }
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674