TEvtGen/EvtGen/EvtGenModels/EvtVubNLO.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:
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtVubNLO.cc
  12 //
  13 // Description: Routine to decay B->Xulnu according to Bosch, Lange, Neubert, and Paz hep-ph/0402094
  14 //              Equation numbers refer to this paper
  15 //
  16 // Modification history:
  17 //
  18 //    Riccardo Faccini       Feb. 11, 2004
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include <stdlib.h>
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtGenKine.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenModels/EvtVubNLO.hh"
  29 #include <string>
  30 #include "EvtGenBase/EvtVector4R.hh"
  31 #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
  32 #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
  33 #include "EvtGenModels/EvtItgPtrFunction.hh"
  34 #include "EvtGenModels/EvtPFermi.hh"
  35 #include "EvtGenBase/EvtRandom.hh"
  36 #include "EvtGenBase/EvtDiLog.hh"
  37
  38 using std::cout;
  39 using std::endl;
  40
  41 EvtVubNLO::~EvtVubNLO() {
  42   delete [] _masses;
  43   delete [] _weights;
  44   cout <<" max pdf : "<<_gmax<<endl;
  45   cout <<" efficiency : "<<(float)_ngood/(float)_ntot<<endl;
  46
  47 }
  48
  49
  50 std::string EvtVubNLO::getName(){
  51
  52   return "VUB_NLO";
  53
  54 }
  55
  56 EvtDecayBase* EvtVubNLO::clone(){
  57
  58   return new EvtVubNLO;
  59
  60 }
  61
  62
  63 void EvtVubNLO::init(){
  64
  65   // max pdf
  66   _gmax=0;
  67   _ntot=0;
  68   _ngood=0;
  69   _lbar=-1000;
  70   _mupi2=-1000;
  71
  72   // check that there are at least 6 arguments
  73   int npar = 8;
  74   if (getNArg()<npar) {
  75
  76     report(ERROR,"EvtGen") << "EvtVubNLO generator expected "
  77                            << " at least npar arguments  but found: "
  78                            <<getNArg()<<endl;
  79     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
  80     ::abort();
  81
  82   }
  83   // this is the shape function parameter
  84   _mb           = getArg(0);
  85   _b            = getArg(1);
  86   _lambdaSF     = getArg(2);// shape function lambda is different from lambda
  87   _mui          = 1.5;// GeV (scale)
  88   _kpar         = getArg(3);// 0
  89   _idSF         = abs((int)getArg(4));// type of shape function 1: exponential (from Neubert)
  90   _nbins        = abs((int)getArg(5));
  91   _masses       = new double[_nbins];
  92   _weights      = new double[_nbins];
  93
  94   // Shape function normalization
  95   _mB=5.28;// temporary B meson mass for normalization
  96
  97   std::vector<double> sCoeffs(11);
  98   sCoeffs[3] = _b;
  99   sCoeffs[4] = _mb;
 100   sCoeffs[5] = _mB;
 101   sCoeffs[6] = _idSF;
 102   sCoeffs[7] =  lambda_SF();
 103   sCoeffs[8] =  mu_h();
 104   sCoeffs[9] =  mu_i();
 105   sCoeffs[10] = 1.;
 106   _SFNorm = SFNorm(sCoeffs) ; // SF normalization;
 107
 108
 109   cout << " pdf 0.66, 1.32 , 4.32 "<<tripleDiff(0.66, 1.32 , 4.32)<<endl;
 110   cout << " pdf 0.23,0.37,3.76 "<<tripleDiff(0.23,0.37,3.76)<<endl;
 111   cout << " pdf 0.97,4.32,4.42 "<<tripleDiff(0.97,4.32,4.42)<<endl;
 112   cout << " pdf 0.52,1.02,2.01 "<<tripleDiff(0.52,1.02,2.01)<<endl;
 113   cout << " pdf 1.35,1.39,2.73 "<<tripleDiff(1.35,1.39,2.73)<<endl;
 114
 115
 116   if (getNArg()-npar+2 != 2*_nbins) {
 117     report(ERROR,"EvtGen") << "EvtVubNLO generator expected "
 118                            << _nbins << " masses and weights but found: "
 119                            <<(getNArg()-npar)/2 <<endl;
 120     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 121     ::abort();
 122   }
 123   int i,j = npar-2;
 124   double maxw = 0.;
 125   for (i=0;i<_nbins;i++) {
 126     _masses[i] = getArg(j++);
 127     if (i>0 && _masses[i] <= _masses[i-1]) {
 128       report(ERROR,"EvtGen") << "EvtVubNLO generator expected "
 129                              << " mass bins in ascending order!"
 130                              << "Will terminate execution!"<<endl;
 131       ::abort();
 132     }
 133     _weights[i] = getArg(j++);
 134     if (_weights[i] < 0) {
 135       report(ERROR,"EvtGen") << "EvtVubNLO generator expected "
 136                              << " weights >= 0, but found: "
 137                              <<_weights[i] <<endl;
 138       report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 139       ::abort();
 140     }
 141     if ( _weights[i] > maxw ) maxw = _weights[i];
 142   }
 143   if (maxw == 0) {
 144     report(ERROR,"EvtGen") << "EvtVubNLO generator expected at least one "
 145                            << " weight > 0, but found none! "
 146                            << "Will terminate execution!"<<endl;
 147     ::abort();
 148   }
 149   for (i=0;i<_nbins;i++) _weights[i]/=maxw;
 150
 151   // the maximum dGamma*p2 value depends on alpha_s only:
 152
 153
 154   //  _dGMax = 0.05;
 155    _dGMax = 150.;
 156
 157   // for the Fermi Motion we need a B-Meso\n mass - but it's not critical
 158   // to get an exact value; in order to stay in the phase space for
 159   // B+- and B0 use the smaller mass
 160
 161
 162   // check that there are 3 daughters
 163   checkNDaug(3);
 164 }
 165
 166 void EvtVubNLO::initProbMax(){
 167
 168   noProbMax();
 169
 170 }
 171
 172 void EvtVubNLO::decay( EvtParticle *p ){
 173
 174   int j;
 175   // B+ -> u-bar specflav l+ nu
 176
 177   EvtParticle *xuhad, *lepton, *neutrino;
 178   EvtVector4R p4;
 179
 180   double pp,pm,pl,ml,El(0.0),Eh(0.0),sh(0.0);
 181
 182
 183
 184   p->initializePhaseSpace(getNDaug(),getDaugs());
 185
 186   xuhad=p->getDaug(0);
 187   lepton=p->getDaug(1);
 188   neutrino=p->getDaug(2);
 189
 190   _mB = p->mass();
 191   ml = lepton->mass();
 192
 193   bool tryit = true;
 194
 195   while (tryit) {
 196     // pm=(E_H+P_H)
 197     pm= EvtRandom::Flat(0.,1);
 198     pm= pow(pm,1./3.)*_mB;
 199     // pl=mB-2*El
 200     pl = EvtRandom::Flat(0.,1);
 201     pl=sqrt(pl)*pm;
 202     // pp=(E_H-P_H)
 203     pp = EvtRandom::Flat(0.,pl);
 204
 205     _ntot++;
 206
 207     El = (_mB-pl)/2.;
 208     Eh = (pp+pm)/2;
 209     sh = pp*pm;
 210
 211     double pdf(0.);
 212     if (pp<pl && El>ml&& sh > _masses[0]*_masses[0]&& _mB*_mB + sh - 2*_mB*Eh > ml*ml) {
 213       double xran = EvtRandom::Flat(0,_dGMax);
 214       pdf = tripleDiff(pp,pl,pm); // triple differential distribution
 215       //      cout <<" P+,P-,Pl,Pdf= "<<pp <<" "<<pm<<" "<<pl<<" "<<pdf<<endl;
 216       if(pdf>_dGMax){
 217         report(ERROR,"EvtGen") << "EvtVubNLO pdf above maximum: " <<pdf
 218                                <<" P+,P-,Pl,Pdf= "<<pp <<" "<<pm<<" "<<pl<<" "<<pdf<<endl;
 219         //::abort();
 220
 221       }
 222       if ( pdf >= xran ) tryit = false;
 223
 224       if(pdf>_gmax)_gmax=pdf;
 225     } else {
 226       //      cout <<" EvtVubNLO incorrect kinematics  sh= "<<sh<<"EH "<<Eh<<endl;
 227     }
 228
 229
 230     // reweight the Mx distribution
 231     if(!tryit && _nbins>0){
 232       _ngood++;
 233       double xran1 = EvtRandom::Flat();
 234       double m = sqrt(sh);j=0;
 235       while ( j < _nbins && m > _masses[j] ) j++;
 236       double w = _weights[j-1];
 237       if ( w < xran1 ) tryit = true;// through away this candidate
 238     }
 239   }
 240
 241   //  cout <<" max prob "<<gmax<<" " << pp<<" "<<y<<" "<<x<<endl;
 242
 243   // o.k. we have the three kineamtic variables
 244   // now calculate a flat cos Theta_H [-1,1] distribution of the
 245   // hadron flight direction w.r.t the B flight direction
 246   // because the B is a scalar and should decay isotropic.
 247   // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
 248   // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
 249   // W flight direction.
 250
 251   double ctH = EvtRandom::Flat(-1,1);
 252   double phH = EvtRandom::Flat(0,2*M_PI);
 253   double phL = EvtRandom::Flat(0,2*M_PI);
 254
 255   // now compute the four vectors in the B Meson restframe
 256
 257   double ptmp,sttmp;
 258   // calculate the hadron 4 vector in the B Meson restframe
 259
 260   sttmp = sqrt(1-ctH*ctH);
 261   ptmp = sqrt(Eh*Eh-sh);
 262   double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
 263   p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
 264   xuhad->init( getDaug(0), p4);
 265
 266
 267   // calculate the W 4 vector in the B Meson restrframe
 268
 269   double apWB = ptmp;
 270   double pWB[4] = {_mB-Eh,-pHB[1],-pHB[2],-pHB[3]};
 271
 272   // first go in the W restframe and calculate the lepton and
 273   // the neutrino in the W frame
 274
 275   double mW2   = _mB*_mB + sh - 2*_mB*Eh;
 276   //  if(mW2<0.1){
 277   //  cout <<" low Q2! "<<pp<<" "<<epp<<" "<<x<<" "<<y<<endl;
 278   //}
 279   double beta  = ptmp/pWB[0];
 280   double gamma = pWB[0]/sqrt(mW2);
 281
 282   double pLW[4];
 283
 284   ptmp = (mW2-ml*ml)/2/sqrt(mW2);
 285   pLW[0] = sqrt(ml*ml + ptmp*ptmp);
 286
 287   double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
 288   if ( ctL < -1 ) ctL = -1;
 289   if ( ctL >  1 ) ctL =  1;
 290   sttmp = sqrt(1-ctL*ctL);
 291
 292   // eX' = eZ x eW
 293   double xW[3] = {-pWB[2],pWB[1],0};
 294   // eZ' = eW
 295   double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
 296
 297   double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
 298   for (j=0;j<2;j++)
 299     xW[j] /= lx;
 300
 301   // eY' = eZ' x eX'
 302   double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
 303   double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
 304   for (j=0;j<3;j++)
 305     yW[j] /= ly;
 306
 307   // p_lep = |p_lep| * (  sin(Theta) * cos(Phi) * eX'
 308   //                    + sin(Theta) * sin(Phi) * eY'
 309   //                    + cos(Theta) *            eZ')
 310   for (j=0;j<3;j++)
 311     pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
 312       +        sttmp*sin(phL)*ptmp*yW[j]
 313       +          ctL         *ptmp*zW[j];
 314
 315   double apLW = ptmp;
 316
 317   // boost them back in the B Meson restframe
 318
 319   double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
 320
 321   ptmp = sqrt(El*El-ml*ml);
 322   double ctLL = appLB/ptmp;
 323
 324   if ( ctLL >  1 ) ctLL =  1;
 325   if ( ctLL < -1 ) ctLL = -1;
 326
 327   double pLB[4] = {El,0,0,0};
 328   double pNB[8] = {pWB[0]-El,0,0,0};
 329
 330   for (j=1;j<4;j++) {
 331     pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
 332     pNB[j] = pWB[j] - pLB[j];
 333   }
 334
 335   p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
 336   lepton->init( getDaug(1), p4);
 337
 338   p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
 339   neutrino->init( getDaug(2), p4);
 340
 341   return ;
 342 }
 343
 344 double
 345 EvtVubNLO::tripleDiff (  double pp, double pl, double pm){
 346
 347   std::vector<double> sCoeffs(11);
 348   sCoeffs[0] = pp;
 349   sCoeffs[1] = pl;
 350   sCoeffs[2] = pm;
 351   sCoeffs[3] = _b;
 352   sCoeffs[4] =  _mb;
 353   sCoeffs[5] = _mB;
 354   sCoeffs[6] = _idSF;
 355   sCoeffs[7] =  lambda_SF();
 356   sCoeffs[8] =  mu_h();
 357   sCoeffs[9] =  mu_i();
 358   sCoeffs[10] =  _SFNorm; // SF normalization;
 359
 360
 361   double c1=(_mB+pl-pp-pm)*(pm-pl);
 362   double c2=2*(pl-pp)*(pm-pl);
 363   double c3=(_mB-pm)*(pm-pp);
 364   double aF1=F10(sCoeffs);
 365   double aF2=F20(sCoeffs);
 366   double aF3=F30(sCoeffs);
 367   double td0=c1*aF1+c2*aF2+c3*aF3;
 368
 369
 370   EvtItgPtrFunction *func = new EvtItgPtrFunction(&integrand, 0., _mB, sCoeffs);
 371   EvtItgAbsIntegrator *jetSF = new EvtItgSimpsonIntegrator(*func,0.01,25);
 372   double smallfrac=0.000001;// stop a bit before the end to avoid problems with numerical integration
 373   double tdInt = jetSF->evaluate(0,pp*(1-smallfrac));
 374   delete jetSF;
 375
 376   double SU=U1lo(mu_h(),mu_i())*pow((pm-pp)/(_mB-pp),alo(mu_h(),mu_i()));
 377   double TD=(_mB-pp)*SU*(td0+tdInt);
 378
 379   return TD;
 380
 381 }
 382
 383 double
 384 EvtVubNLO::integrand(double omega, const std::vector<double> &coeffs){
 385   //double pp=coeffs[0];
 386   double c1=(coeffs[5]+coeffs[1]-coeffs[0]-coeffs[2])*(coeffs[2]-coeffs[1]);
 387   double c2=2*(coeffs[1]-coeffs[0])*(coeffs[2]-coeffs[1]);
 388   double c3=(coeffs[5]-coeffs[2])*(coeffs[2]-coeffs[0]);
 389
 390   return  c1*F1Int(omega,coeffs)+c2*F2Int(omega,coeffs)+c3*F3Int(omega,coeffs);
 391 }
 392
 393 double
 394 EvtVubNLO::F10(const std::vector<double> &coeffs){
 395   double pp=coeffs[0];
 396   double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
 397   double mui=coeffs[9];
 398   double muh=coeffs[8];
 399   double z=1-y;
 400   double result= U1nlo(muh,mui)/ U1lo(muh,mui);
 401
 402   result += anlo(muh,mui)*log(y);
 403
 404   result += C_F(muh)*(-4*pow(log(y*coeffs[4]/muh),2)+10*log(y*coeffs[4]/muh)-4*log(y)-2*log(y)/(1-y)-4.0*EvtDiLog::DiLog(z)-pow(EvtConst::pi,2)/6.-12  );
 405
 406   result += C_F(mui)*(2*pow(log(y*coeffs[4]*pp/pow(mui,2)),2)-3*log(y*coeffs[4]*pp/pow(mui,2))+7-pow(EvtConst::pi,2) );
 407   result *=shapeFunction(pp,coeffs);
 408   // changes due to SSF
 409   result +=  (-subS(coeffs)+2*subT(coeffs)+(subU(coeffs)-subV(coeffs))*(1/y-1.))/(coeffs[5]-pp);
 410   result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0]),2)*(-5*(lambda1()+3*lambda2())/6+2*(2*lambda1()/3-lambda2())/pow(y,2));
 411   //  result +=  (subS(coeffs)+subT(coeffs)+(subU(coeffs)-subV(coeffs))/y)/(coeffs[5]-pp);
 412   // this part has been added after Feb '05
 413
 414   //result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0]),2)*((lambda1()+3*lambda2())/6+2*(2*lambda1()/3-lambda2())/pow(y,2));
 415   return result;
 416 }
 417
 418 double
 419 EvtVubNLO::F1Int(double omega,const std::vector<double> &coeffs){
 420   double pp=coeffs[0];
 421   double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
 422   // mubar == mui
 423   return C_F(coeffs[9])*(
 424                          (shapeFunction(omega,coeffs)-shapeFunction(pp,coeffs))*(4*log(y*coeffs[4]*(pp-omega)/pow(coeffs[9],2))-3)/(pp-omega)+
 425                          (g1(y,(pp-omega)/(coeffs[5]-coeffs[0]))/(coeffs[5]-pp)*shapeFunction(omega,coeffs))
 426                          );
 427 }
 428
 429 double
 430 EvtVubNLO::F20(const std::vector<double> &coeffs){
 431   double pp=coeffs[0];
 432   double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
 433   double result= C_F(coeffs[8])*log(y)/(1-y)*shapeFunction(pp,coeffs)-
 434     1/y*(subS(coeffs)+2*subT(coeffs)-(subT(coeffs)+subV(coeffs))/y)/(coeffs[5]-pp);
 435   // added after Feb '05
 436   result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0])*y,2)*(2*lambda1()/3+4*lambda2()-y*(7/6*lambda1()+3*lambda2()));
 437   return result;
 438 }
 439
 440 double
 441 EvtVubNLO::F2Int(double omega,const std::vector<double> &coeffs){
 442   double pp=coeffs[0];
 443   double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
 444   return C_F(coeffs[9])*g3(y,(pp-omega)/(coeffs[5]-coeffs[0]))*shapeFunction(omega,coeffs)/(coeffs[5]-pp);
 445 }
 446
 447 double
 448 EvtVubNLO::F30(const std::vector<double> &coeffs){
 449   double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
 450   return shapeFunction(coeffs[0],coeffs)/pow((coeffs[5]-coeffs[0])*y,2)*(-2*lambda1()/3+lambda2());
 451 }
 452
 453 double
 454 EvtVubNLO::F3Int(double omega,const std::vector<double> &coeffs){
 455   double pp=coeffs[0];
 456   double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
 457   return C_F(coeffs[9])*g3(y,(pp-omega)/(coeffs[5]-coeffs[0]))/2*shapeFunction(omega,coeffs)/(coeffs[2]-coeffs[0]);
 458 }
 459
 460 double
 461 EvtVubNLO::g1(double y, double x){
 462   double result=(y*(-9+10*y)+x*x*(-12+13*y)+2*x*(-8+6*y+3*y*y))/y/pow(1+x,2)/(x+y);
 463   result -= 4*log((1+1/x)*y)/x;
 464   result -=2*log(1+y/x)*(3*pow(x,4)*(-2+y)-2*pow(y,3)-4*pow(x,3)*(2+y)-2*x*y*y*(4+y)-x*x*y*(12+4*y+y*y))/x/pow((1+x)*y,2)/(x+y);
 465   return result;
 466 }
 467
 468 double
 469 EvtVubNLO::g2(double y, double x){
 470   double result=y*(10*pow(x,4)+y*y+3*x*x*y*(10+y)+pow(x,3)*(12+19*y)+x*y*(8+4*y+y*y));
 471   result -= 2*x*log(1+y/x)*(5*pow(x,4)+2*y*y+6*pow(x,3)*(1+2*y)+4*y*x*(1+2*y)+x*x*y*(18+5*y));
 472   result *= 2/(pow(y*(1+x),2)*y*(x+y));
 473   return result;
 474 }
 475
 476 double
 477 EvtVubNLO::g3(double y, double x){
 478   double result=(2*pow(y,3)*(-11+2*y)-10*pow(x,4)*(6-6*y+y*y)+x*y*y*(-94+29*y+2*y*y)+2*x*x*y*(-72+18*y+13*y*y)-pow(x,3)*(72+42*y-70*y*y+3*pow(y,3)))/(pow(y*(1+x),2)*y*(x+y));
 479   result += 2*log(1+y/x)*(-6*x*pow(y,3)*(-5+y)+4*pow(y,4)+5*pow(x,5)*(6-6*y+y*y)-4*pow(x*y,2)*(-20+6*y+y*y)+pow(x,3)*y*(90-10*y-28*y*y+pow(y,3))+pow(x,4)*(36+36*y-50*y*y+4*pow(y,3)))/(pow((1+x)*y*y,2)*(x+y));
 480   return result;
 481 }
 482
 483 /* old version (before Feb 05 notebook from NNeubert
 484
 485 double
 486 EvtVubNLO::F1Int(double omega,const std::vector<double> &coeffs){
 487   double pp=coeffs[0];
 488   double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
 489   // mubar == mui
 490   return C_F(coeffs[9])*(
 491                          (shapeFunction(omega,coeffs)-shapeFunction(pp,coeffs))*(4*log(y*coeffs[4]*(pp-omega)/pow(coeffs[9],2))-3)/(pp-omega)-
 492                          (1./y/(coeffs[5]-pp)*shapeFunction(omega,coeffs)*(5-6*y+4*(3-y)*log((pp-omega)/y/coeffs[4])))
 493                          );
 494 }
 495
 496
 497 double
 498 EvtVubNLO::F2Int(double omega,const std::vector<double> &coeffs){
 499   double pp=coeffs[0];
 500   double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
 501   return C_F(coeffs[9])*shapeFunction(omega,coeffs)*(2-11/y-4/y*log((pp-omega)/y/coeffs[4]))/(coeffs[5]-pp);
 502 }
 503
 504 double
 505 EvtVubNLO::F3(const std::vector<double> &coeffs){
 506   return C_F(coeffs[9])*shapeFunction(omega,coeffs)/(coeffs[2]-coeffs[0]);
 507 }
 508 */
 509
 510 double EvtVubNLO::SFNorm(  const std::vector<double> &/*coeffs*/){
 511
 512   double omega0=1.68;//normalization scale (mB-2*1.8)
 513   if(_idSF==1){ // exponential SF
 514     double omega0=1.68;//normalization scale (mB-2*1.8)
 515     return M0(mu_i(),omega0)*pow(_b,_b)/lambda_SF()/ (Gamma(_b)-Gamma(_b,_b*omega0/lambda_SF()));
 516   } else if(_idSF==2){ // Gaussian SF
 517     double c=cGaus(_b);
 518     return M0(mu_i(),omega0)*2/lambda_SF()/pow(c,-(1+_b)/2.)/
 519       (Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c));
 520   } else {
 521     report(ERROR,"EvtGen")  << "unknown SF "<<_idSF<<endl;
 522     return -1;
 523   }
 524 }
 525
 526 double
 527 EvtVubNLO::shapeFunction ( double omega, const std::vector<double> &sCoeffs){
 528   if( sCoeffs[6]==1){
 529     return sCoeffs[10]*expShapeFunction(omega, sCoeffs);
 530   } else if( sCoeffs[6]==2) {
 531     return sCoeffs[10]*gausShapeFunction(omega, sCoeffs);
 532   } else {
 533  report(ERROR,"EvtGen") << "EvtVubNLO : unknown shape function # "
 534                            <<sCoeffs[6]<<endl;
 535   }
 536   return -1.;
 537 }
 538
 539
 540 // SSF
 541 double
 542 EvtVubNLO::subS(const std::vector<double> &c){ return (lambda_bar(1.68)-c[0])*shapeFunction(c[0],c);}
 543 double
 544 EvtVubNLO::subT(const std::vector<double> &c){  return -3*lambda2()*subS(c)/mu_pi2(1.68);}
 545 double
 546 EvtVubNLO::subU(const std::vector<double> &c){ return -2*subS(c);}
 547 double
 548 EvtVubNLO::subV(const std::vector<double> &c){ return -subT(c);}
 549
 550
 551 double
 552 EvtVubNLO::lambda_bar(double omega0){
 553   if(_lbar<0){
 554     if(_idSF==1){ // exponential SF
 555       double rat=omega0*_b/lambda_SF();
 556       _lbar=lambda_SF()/_b*(Gamma(1+_b)-Gamma(1+_b,rat))/(Gamma(_b)-Gamma(_b,rat));
 557     } else if(_idSF==2){ // Gaussian SF
 558       double c=cGaus(_b);
 559       _lbar=lambda_SF()*(Gamma(1+_b/2)-Gamma(1+_b/2,pow(omega0/lambda_SF(),2)*c))/(Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c))/sqrt(c);
 560     }
 561   }
 562   return _lbar;
 563 }
 564
 565
 566 double
 567 EvtVubNLO::mu_pi2(double omega0){
 568   if(_mupi2<0){
 569     if(_idSF==1){ // exponential SF
 570       double rat=omega0*_b/lambda_SF();
 571       _mupi2= 3*(pow(lambda_SF()/_b,2)*(Gamma(2+_b)-Gamma(2+_b,rat))/(Gamma(_b)-Gamma(_b,rat))-pow(lambda_bar(omega0),2));
 572     } else if(_idSF==2){ // Gaussian SF
 573       double c=cGaus(_b);
 574       double m1=Gamma((3+_b)/2)-Gamma((3+_b)/2,pow(omega0/lambda_SF(),2)*c);
 575       double m2=Gamma(1+_b/2)-Gamma(1+_b/2,pow(omega0/lambda_SF(),2)*c);
 576       double m3=Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c);
 577       _mupi2= 3*pow(lambda_SF(),2)*(m1/m3-pow(m2/m3,2))/c;
 578     }
 579   }
 580   return _mupi2;
 581 }
 582
 583 double
 584 EvtVubNLO::M0(double mui,double omega0){
 585   double mf=omega0-lambda_bar(omega0);
 586   return 1+4*C_F(mui)*(-pow(log(mf/mui),2)-log(mf/mui)-pow(EvtConst::pi/2,2)/6.+mu_pi2(omega0)/3/pow(mf,2)*(log(mf/mui)-0.5));
 587 }
 588
 589 double
 590 EvtVubNLO::alphas(double mu){
 591   double Lambda4=0.302932;
 592   double lg=2*log(mu/Lambda4);
 593   return 4*EvtConst::pi/lg/beta0()*(1-beta1()*log(lg)/pow(beta0(),2)/lg+pow(beta1()/lg,2)/pow(beta0(),4)*(pow(log(lg)-0.5,2)-1.25+beta2()*beta0()/pow(beta1(),2)));
 594 }
 595
 596 double
 597 EvtVubNLO::gausShapeFunction ( double omega, const std::vector<double> &sCoeffs){
 598   double b=sCoeffs[3];
 599   double l=sCoeffs[7];
 600   double wL=omega/l;
 601
 602   return pow(wL,b)*exp(-cGaus(b)*wL*wL);
 603 }
 604
 605 double
 606 EvtVubNLO::expShapeFunction ( double omega, const std::vector<double> &sCoeffs){
 607   double b=sCoeffs[3];
 608   double l=sCoeffs[7];
 609   double wL=omega/l;
 610
 611   return pow(wL,b-1)*exp(-b*wL);
 612 }
 613
 614 double
 615 EvtVubNLO::Gamma(double z) {
 616
 617   std::vector<double> gammaCoeffs(6);
 618   gammaCoeffs[0]=76.18009172947146;
 619   gammaCoeffs[1]=-86.50532032941677;
 620   gammaCoeffs[2]=24.01409824083091;
 621   gammaCoeffs[3]=-1.231739572450155;
 622   gammaCoeffs[4]=0.1208650973866179e-2;
 623   gammaCoeffs[5]=-0.5395239384953e-5;
 624
 625   //Lifted from Numerical Recipies in C
 626   double x, y, tmp, ser;
 627
 628   int j;
 629   y = z;
 630   x = z;
 631
 632   tmp = x + 5.5;
 633   tmp = tmp - (x+0.5)*log(tmp);
 634   ser=1.000000000190015;
 635
 636   for (j=0;j<6;j++) {
 637     y = y +1.0;
 638     ser = ser + gammaCoeffs[j]/y;
 639   }
 640
 641   return exp(-tmp+log(2.5066282746310005*ser/x));
 642
 643 }
 644
 645
 646
 647 double
 648 EvtVubNLO::Gamma(double z, double tmin) {
 649   std::vector<double> c(1);
 650   c[0]=z;
 651   EvtItgPtrFunction *func = new EvtItgPtrFunction(&dgamma, tmin, 100., c);
 652   EvtItgAbsIntegrator *jetSF = new EvtItgSimpsonIntegrator(*func,0.001);
 653   return jetSF->evaluate(tmin,100.);
 654 }