EvtVubNLO.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information:
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtVubNLO.cc
//
// Description: Routine to decay B->Xulnu according to Bosch, Lange, Neubert, and Paz
// hep-ph/0402094
//              Equation numbers refer to this paper
//
// Modification history:
//
//    Riccardo Faccini       Feb. 11, 2004
//
//------------------------------------------------------------------------
//
#include "EvtVubNLO.hh"
#include "../EvtGenBase/EvtGenKine.hh"
#include "../EvtGenBase/EvtPDL.hh"
#include "../EvtGenBase/EvtParticle.hh"
#include "../EvtGenBase/EvtPatches.hh"
#include "../EvtGenBase/EvtRandom.hh"
#include "../EvtGenBase/EvtReport.hh"
#include "../EvtGenBase/EvtVector4R.hh"
#include "EvtBtoXsgammaFermiUtil.hh"
#include "EvtItgPtrFunction.hh"
#include "EvtItgSimpsonIntegrator.hh"
#include "EvtPFermi.hh"
#include <stdlib.h>
#include <string>
using std::cout;
using std::endl;
 
EvtVubNLO::~EvtVubNLO() {
  delete[] _masses;
  delete[] _weights;
  cout << " max pdf : " << _gmax << endl;
  cout << " efficiency : " << (float)_ngood / (float)_ntot << endl;
}
 
extern "C" {
double ddilog_( const double* );
}
 
void EvtVubNLO::getName( std::string& model_name ) { model_name = "VUB_NLO"; }
 
EvtDecayBase* EvtVubNLO::clone() { return new EvtVubNLO; }
 
void EvtVubNLO::init() {
 
  // max pdf
  _gmax  = 0;
  _ntot  = 0;
  _ngood = 0;
  _lbar  = -1000;
  _mupi2 = -1000;
 
  // check that there are at least 6 arguments
  int npar = 8;
  if ( getNArg() < npar )
  {
 
    report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected "
                              << " at least npar arguments  but found: " << getNArg() << endl;
    report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
    ::abort();
  }
  // this is the shape function parameter
  _mb       = getArg( 0 );
  _b        = getArg( 1 );
  _lambdaSF = getArg( 2 );             // shape function lambda is different from lambda
  _mui      = 1.5;                     // GeV (scale)
  _kpar     = getArg( 3 );             // 0
  _idSF     = abs( (int)getArg( 4 ) ); // type of shape function 1: exponential (from Neubert)
  _nbins    = abs( (int)getArg( 5 ) );
  _masses   = new double[_nbins];
  _weights  = new double[_nbins];
 
  // Shape function normalization
  _mB = 5.28; // temporary B meson mass for normalization
 
  std::vector<double> sCoeffs( 11 );
  sCoeffs[3]  = _b;
  sCoeffs[4]  = _mb;
  sCoeffs[5]  = _mB;
  sCoeffs[6]  = _idSF;
  sCoeffs[7]  = lambda_SF();
  sCoeffs[8]  = mu_h();
  sCoeffs[9]  = mu_i();
  sCoeffs[10] = 1.;
  _SFNorm     = SFNorm( sCoeffs ); // SF normalization;
 
  cout << " pdf 0.66, 1.32 , 4.32 " << tripleDiff( 0.66, 1.32, 4.32 ) << endl;
  cout << " pdf 0.23,0.37,3.76 " << tripleDiff( 0.23, 0.37, 3.76 ) << endl;
  cout << " pdf 0.97,4.32,4.42 " << tripleDiff( 0.97, 4.32, 4.42 ) << endl;
  cout << " pdf 0.52,1.02,2.01 " << tripleDiff( 0.52, 1.02, 2.01 ) << endl;
  cout << " pdf 1.35,1.39,2.73 " << tripleDiff( 1.35, 1.39, 2.73 ) << endl;
 
  if ( getNArg() - npar + 2 != 2 * _nbins )
  {
    report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected " << _nbins
                              << " masses and weights but found: " << ( getNArg() - npar ) / 2
                              << endl;
    report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
    ::abort();
  }
  int    i, j = npar - 2;
  double maxw = 0.;
  for ( i = 0; i < _nbins; i++ )
  {
    _masses[i] = getArg( j++ );
    if ( i > 0 && _masses[i] <= _masses[i - 1] )
    {
      report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected "
                                << " mass bins in ascending order!"
                                << "Will terminate execution!" << endl;
      ::abort();
    }
    _weights[i] = getArg( j++ );
    if ( _weights[i] < 0 )
    {
      report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected "
                                << " weights >= 0, but found: " << _weights[i] << endl;
      report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
      ::abort();
    }
    if ( _weights[i] > maxw ) maxw = _weights[i];
  }
  if ( maxw == 0 )
  {
    report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected at least one "
                              << " weight > 0, but found none! "
                              << "Will terminate execution!" << endl;
    ::abort();
  }
  for ( i = 0; i < _nbins; i++ ) _weights[i] /= maxw;
 
  // the maximum dGamma*p2 value depends on alpha_s only:
 
  //  _dGMax = 0.05;
  _dGMax = 150.;
 
  // for the Fermi Motion we need a B-Meso\n mass - but it's not critical
  // to get an exact value; in order to stay in the phase space for
  // B+- and B0 use the smaller mass
 
  // check that there are 3 daughters
  checkNDaug( 3 );
}
 
void EvtVubNLO::initProbMax() { noProbMax(); }
 
void EvtVubNLO::decay( EvtParticle* p ) {
 
  int j;
  // B+ -> u-bar specflav l+ nu
 
  EvtParticle *xuhad, *lepton, *neutrino;
  EvtVector4R  p4;
 
  double pp, pm, pl, ml, El( 0.0 ), Eh( 0.0 ), sh( 0.0 );
 
  p->initializePhaseSpace( getNDaug(), getDaugs() );
 
  xuhad    = p->getDaug( 0 );
  lepton   = p->getDaug( 1 );
  neutrino = p->getDaug( 2 );
 
  _mB = p->mass();
  ml  = lepton->mass();
 
  bool tryit = true;
 
  while ( tryit )
  {
    // pm=(E_H+P_H)
    pm = EvtRandom::Flat( 0., 1 );
    pm = pow( pm, 1. / 3. ) * _mB;
    // pl=mB-2*El
    pl = EvtRandom::Flat( 0., 1 );
    pl = sqrt( pl ) * pm;
    // pp=(E_H-P_H)
    pp = EvtRandom::Flat( 0., pl );
 
    _ntot++;
 
    El = ( _mB - pl ) / 2.;
    Eh = ( pp + pm ) / 2;
    sh = pp * pm;
 
    double pdf( 0. );
    if ( pp < pl && El > ml && sh > _masses[0] * _masses[0] &&
         _mB * _mB + sh - 2 * _mB * Eh > ml * ml )
    {
      double xran = EvtRandom::Flat( 0, _dGMax );
      pdf         = tripleDiff( pp, pl, pm ); // triple differential distribution
      //      cout <<" P+,P-,Pl,Pdf= "<<pp <<" "<<pm<<" "<<pl<<" "<<pdf<<endl;
      if ( pdf > _dGMax )
      {
        report( ERROR, "EvtGen" )
            << "EvtVubNLO pdf above maximum: " << pdf << " P+,P-,Pl,Pdf= " << pp << " " << pm
            << " " << pl << " " << pdf << endl;
        //::abort();
      }
      if ( pdf >= xran ) tryit = false;
 
      if ( pdf > _gmax ) _gmax = pdf;
    }
    else
    {
      //      cout <<" EvtVubNLO incorrect kinematics  sh= "<<sh<<"EH "<<Eh<<endl;
    }
 
    // reweight the Mx distribution
    if ( !tryit && _nbins > 0 )
    {
      _ngood++;
      double xran1 = EvtRandom::Flat();
      double m     = sqrt( sh );
      j            = 0;
      while ( j < _nbins && m > _masses[j] ) j++;
      double w = _weights[j - 1];
      if ( w < xran1 ) tryit = true; // through away this candidate
    }
  }
 
  //  cout <<" max prob "<<gmax<<" " << pp<<" "<<y<<" "<<x<<endl;
 
  // o.k. we have the three kineamtic variables
  // now calculate a flat cos Theta_H [-1,1] distribution of the
  // hadron flight direction w.r.t the B flight direction
  // because the B is a scalar and should decay isotropic.
  // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
  // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
  // W flight direction.
 
  double ctH = EvtRandom::Flat( -1, 1 );
  double phH = EvtRandom::Flat( 0, 2 * M_PI );
  double phL = EvtRandom::Flat( 0, 2 * M_PI );
 
  // now compute the four vectors in the B Meson restframe
 
  double ptmp, sttmp;
  // calculate the hadron 4 vector in the B Meson restframe
 
  sttmp         = sqrt( 1 - ctH * ctH );
  ptmp          = sqrt( Eh * Eh - sh );
  double pHB[4] = { Eh, ptmp * sttmp * cos( phH ), ptmp * sttmp * sin( phH ), ptmp * ctH };
  p4.set( pHB[0], pHB[1], pHB[2], pHB[3] );
  xuhad->init( getDaug( 0 ), p4 );
 
  // calculate the W 4 vector in the B Meson restrframe
 
  double apWB   = ptmp;
  double pWB[4] = { _mB - Eh, -pHB[1], -pHB[2], -pHB[3] };
 
  // first go in the W restframe and calculate the lepton and
  // the neutrino in the W frame
 
  double mW2 = _mB * _mB + sh - 2 * _mB * Eh;
  //  if(mW2<0.1){
  //  cout <<" low Q2! "<<pp<<" "<<epp<<" "<<x<<" "<<y<<endl;
  //}
  double beta  = ptmp / pWB[0];
  double gamma = pWB[0] / sqrt( mW2 );
 
  double pLW[4];
 
  ptmp   = ( mW2 - ml * ml ) / 2 / sqrt( mW2 );
  pLW[0] = sqrt( ml * ml + ptmp * ptmp );
 
  double ctL = ( El - gamma * pLW[0] ) / beta / gamma / ptmp;
  if ( ctL < -1 ) ctL = -1;
  if ( ctL > 1 ) ctL = 1;
  sttmp = sqrt( 1 - ctL * ctL );
 
  // eX' = eZ x eW
  double xW[3] = { -pWB[2], pWB[1], 0 };
  // eZ' = eW
  double zW[3] = { pWB[1] / apWB, pWB[2] / apWB, pWB[3] / apWB };
 
  double lx = sqrt( xW[0] * xW[0] + xW[1] * xW[1] );
  for ( j = 0; j < 2; j++ ) xW[j] /= lx;
 
  // eY' = eZ' x eX'
  double yW[3] = { -pWB[1] * pWB[3], -pWB[2] * pWB[3], pWB[1] * pWB[1] + pWB[2] * pWB[2] };
  double ly    = sqrt( yW[0] * yW[0] + yW[1] * yW[1] + yW[2] * yW[2] );
  for ( j = 0; j < 3; j++ ) yW[j] /= ly;
 
  // p_lep = |p_lep| * (  sin(Theta) * cos(Phi) * eX'
  //                    + sin(Theta) * sin(Phi) * eY'
  //                    + cos(Theta) *            eZ')
  for ( j = 0; j < 3; j++ )
    pLW[j + 1] = sttmp * cos( phL ) * ptmp * xW[j] + sttmp * sin( phL ) * ptmp * yW[j] +
                 ctL * ptmp * zW[j];
 
  double apLW = ptmp;
 
  // boost them back in the B Meson restframe
 
  double appLB = beta * gamma * pLW[0] + gamma * ctL * apLW;
 
  ptmp        = sqrt( El * El - ml * ml );
  double ctLL = appLB / ptmp;
 
  if ( ctLL > 1 ) ctLL = 1;
  if ( ctLL < -1 ) ctLL = -1;
 
  double pLB[4] = { El, 0, 0, 0 };
  double pNB[8] = { pWB[0] - El, 0, 0, 0 };
 
  for ( j = 1; j < 4; j++ )
  {
    pLB[j] = pLW[j] + ( ctLL * ptmp - ctL * apLW ) / apWB * pWB[j];
    pNB[j] = pWB[j] - pLB[j];
  }
 
  p4.set( pLB[0], pLB[1], pLB[2], pLB[3] );
  lepton->init( getDaug( 1 ), p4 );
 
  p4.set( pNB[0], pNB[1], pNB[2], pNB[3] );
  neutrino->init( getDaug( 2 ), p4 );
 
  return;
}
 
double EvtVubNLO::tripleDiff( double pp, double pl, double pm ) {
 
  std::vector<double> sCoeffs( 11 );
  sCoeffs[0]  = pp;
  sCoeffs[1]  = pl;
  sCoeffs[2]  = pm;
  sCoeffs[3]  = _b;
  sCoeffs[4]  = _mb;
  sCoeffs[5]  = _mB;
  sCoeffs[6]  = _idSF;
  sCoeffs[7]  = lambda_SF();
  sCoeffs[8]  = mu_h();
  sCoeffs[9]  = mu_i();
  sCoeffs[10] = _SFNorm; // SF normalization;
 
  double c1  = ( _mB + pl - pp - pm ) * ( pm - pl );
  double c2  = 2 * ( pl - pp ) * ( pm - pl );
  double c3  = ( _mB - pm ) * ( pm - pp );
  double aF1 = F10( sCoeffs );
  double aF2 = F20( sCoeffs );
  double aF3 = F30( sCoeffs );
  double td0 = c1 * aF1 + c2 * aF2 + c3 * aF3;
 
  EvtItgPtrFunction*   func  = new EvtItgPtrFunction( &integrand, 0., _mB, sCoeffs );
  EvtItgAbsIntegrator* jetSF = new EvtItgSimpsonIntegrator( *func, 0.01, 25 );
  double               smallfrac =
      0.000001; // stop a bit before the end to avoid problems with numerical integration
  double tdInt = jetSF->evaluate( 0, pp * ( 1 - smallfrac ) );
  delete jetSF;
 
  double SU =
      U1lo( mu_h(), mu_i() ) * pow( ( pm - pp ) / ( _mB - pp ), alo( mu_h(), mu_i() ) );
  double TD = ( _mB - pp ) * SU * ( td0 + tdInt );
 
  return TD;
}
 
double EvtVubNLO::integrand( double omega, const std::vector<double>& coeffs ) {
  // double pp=coeffs[0];
  double c1 = ( coeffs[5] + coeffs[1] - coeffs[0] - coeffs[2] ) * ( coeffs[2] - coeffs[1] );
  double c2 = 2 * ( coeffs[1] - coeffs[0] ) * ( coeffs[2] - coeffs[1] );
  double c3 = ( coeffs[5] - coeffs[2] ) * ( coeffs[2] - coeffs[0] );
 
  return c1 * F1Int( omega, coeffs ) + c2 * F2Int( omega, coeffs ) +
         c3 * F3Int( omega, coeffs );
}
 
double EvtVubNLO::F10( const std::vector<double>& coeffs ) {
  double pp     = coeffs[0];
  double y      = ( coeffs[2] - coeffs[0] ) / ( coeffs[5] - coeffs[0] );
  double mui    = coeffs[9];
  double muh    = coeffs[8];
  double z      = 1 - y;
  double result = U1nlo( muh, mui ) / U1lo( muh, mui );
 
  result += anlo( muh, mui ) * log( y );
 
  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 * ddilog_( &z ) -
              pow( EvtConst::pi, 2 ) / 6. - 12 );
 
  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 ) );
  result *= shapeFunction( pp, coeffs );
  // changes due to SSF
  result += ( -subS( coeffs ) + 2 * subT( coeffs ) +
              ( subU( coeffs ) - subV( coeffs ) ) * ( 1 / y - 1. ) ) /
            ( coeffs[5] - pp );
  result += shapeFunction( pp, coeffs ) / pow( ( coeffs[5] - coeffs[0] ), 2 ) *
            ( -5 * ( lambda1() + 3 * lambda2() ) / 6 +
              2 * ( 2 * lambda1() / 3 - lambda2() ) / pow( y, 2 ) );
  //  result +=  (subS(coeffs)+subT(coeffs)+(subU(coeffs)-subV(coeffs))/y)/(coeffs[5]-pp);
  // this part has been added after Feb '05
 
  // result +=
  // shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0]),2)*((lambda1()+3*lambda2())/6+2*(2*lambda1()/3-lambda2())/pow(y,2));
  return result;
}
 
double EvtVubNLO::F1Int( double omega, const std::vector<double>& coeffs ) {
  double pp = coeffs[0];
  double y  = ( coeffs[2] - coeffs[0] ) / ( coeffs[5] - coeffs[0] );
  // mubar == mui
  return C_F( coeffs[9] ) *
         ( ( shapeFunction( omega, coeffs ) - shapeFunction( pp, coeffs ) ) *
               ( 4 * log( y * coeffs[4] * ( pp - omega ) / pow( coeffs[9], 2 ) ) - 3 ) /
               ( pp - omega ) +
           ( g1( y, ( pp - omega ) / ( coeffs[5] - coeffs[0] ) ) / ( coeffs[5] - pp ) *
             shapeFunction( omega, coeffs ) ) );
}
 
double EvtVubNLO::F20( const std::vector<double>& coeffs ) {
  double pp = coeffs[0];
  double y  = ( coeffs[2] - coeffs[0] ) / ( coeffs[5] - coeffs[0] );
  double result =
      C_F( coeffs[8] ) * log( y ) / ( 1 - y ) * shapeFunction( pp, coeffs ) -
      1 / y *
          ( subS( coeffs ) + 2 * subT( coeffs ) - ( subT( coeffs ) + subV( coeffs ) ) / y ) /
          ( coeffs[5] - pp );
  // added after Feb '05
  result += shapeFunction( pp, coeffs ) / pow( ( coeffs[5] - coeffs[0] ) * y, 2 ) *
            ( 2 * lambda1() / 3 + 4 * lambda2() - y * ( 7 / 6 * lambda1() + 3 * lambda2() ) );
  return result;
}
 
double EvtVubNLO::F2Int( double omega, const std::vector<double>& coeffs ) {
  double pp = coeffs[0];
  double y  = ( coeffs[2] - coeffs[0] ) / ( coeffs[5] - coeffs[0] );
  return C_F( coeffs[9] ) * g3( y, ( pp - omega ) / ( coeffs[5] - coeffs[0] ) ) *
         shapeFunction( omega, coeffs ) / ( coeffs[5] - pp );
}
 
double EvtVubNLO::F30( const std::vector<double>& coeffs ) {
  double y = ( coeffs[2] - coeffs[0] ) / ( coeffs[5] - coeffs[0] );
  return shapeFunction( coeffs[0], coeffs ) / pow( ( coeffs[5] - coeffs[0] ) * y, 2 ) *
         ( -2 * lambda1() / 3 + lambda2() );
}
 
double EvtVubNLO::F3Int( double omega, const std::vector<double>& coeffs ) {
  double pp = coeffs[0];
  double y  = ( coeffs[2] - coeffs[0] ) / ( coeffs[5] - coeffs[0] );
  return C_F( coeffs[9] ) * g3( y, ( pp - omega ) / ( coeffs[5] - coeffs[0] ) ) / 2 *
         shapeFunction( omega, coeffs ) / ( coeffs[2] - coeffs[0] );
}
 
double EvtVubNLO::g1( double y, double x ) {
  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 );
  result -= 4 * log( ( 1 + 1 / x ) * y ) / x;
  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 );
  return result;
}
 
double EvtVubNLO::g2( double y, double x ) {
  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 ) );
  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 ) );
  result *= 2 / ( pow( y * ( 1 + x ), 2 ) * y * ( x + y ) );
  return result;
}
 
double EvtVubNLO::g3( double y, double x ) {
  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 ) );
  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 ) );
  return result;
}
 
/* old version (before Feb 05 notebook from NNeubert
 
double
EvtVubNLO::F1Int(double omega,const std::vector<double> &coeffs){
  double pp=coeffs[0];
  double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
  // mubar == mui
  return C_F(coeffs[9])*(
                         (shapeFunction(omega,coeffs)-shapeFunction(pp,coeffs))*(4*log(y*coeffs[4]*(pp-omega)/pow(coeffs[9],2))-3)/(pp-omega)-
                         (1./y/(coeffs[5]-pp)*shapeFunction(omega,coeffs)*(5-6*y+4*(3-y)*log((pp-omega)/y/coeffs[4])))
                         );
}
 
 
double
EvtVubNLO::F2Int(double omega,const std::vector<double> &coeffs){
  double pp=coeffs[0];
  double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
  return
C_F(coeffs[9])*shapeFunction(omega,coeffs)*(2-11/y-4/y*log((pp-omega)/y/coeffs[4]))/(coeffs[5]-pp);
}
 
double
EvtVubNLO::F3(const std::vector<double> &coeffs){
  return C_F(coeffs[9])*shapeFunction(omega,coeffs)/(coeffs[2]-coeffs[0]);
}
*/
 
double EvtVubNLO::SFNorm( const std::vector<double>& coeffs ) {
 
  double omega0 = 1.68; // normalization scale (mB-2*1.8)
  if ( _idSF == 1 )
  {                       // exponential SF
    double omega0 = 1.68; // normalization scale (mB-2*1.8)
    return M0( mu_i(), omega0 ) * pow( _b, _b ) / lambda_SF() /
           ( Gamma( _b ) - Gamma( _b, _b * omega0 / lambda_SF() ) );
  }
  else if ( _idSF == 2 )
  { // Gaussian SF
    double c = cGaus( _b );
    return M0( mu_i(), omega0 ) * 2 / lambda_SF() / pow( c, -( 1 + _b ) / 2. ) /
           ( Gamma( ( 1 + _b ) / 2 ) -
             Gamma( ( 1 + _b ) / 2, pow( omega0 / lambda_SF(), 2 ) * c ) );
  }
  else
  {
    report( ERROR, "EvtGen" ) << "unknown SF " << _idSF << endl;
    return -1;
  }
}
 
double EvtVubNLO::shapeFunction( double omega, const std::vector<double>& sCoeffs ) {
  if ( sCoeffs[6] == 1 ) { return sCoeffs[10] * expShapeFunction( omega, sCoeffs ); }
  else if ( sCoeffs[6] == 2 ) { return sCoeffs[10] * gausShapeFunction( omega, sCoeffs ); }
  else
  {
    report( ERROR, "EvtGen" ) << "EvtVubNLO : unknown shape function # " << sCoeffs[6] << endl;
  }
  return -1.;
}
 
// SSF
double EvtVubNLO::subS( const std::vector<double>& c ) {
  return ( lambda_bar( 1.68 ) - c[0] ) * shapeFunction( c[0], c );
}
double EvtVubNLO::subT( const std::vector<double>& c ) {
  return -3 * lambda2() * subS( c ) / mu_pi2( 1.68 );
}
double EvtVubNLO::subU( const std::vector<double>& c ) { return -2 * subS( c ); }
double EvtVubNLO::subV( const std::vector<double>& c ) { return -subT( c ); }
 
double EvtVubNLO::lambda_bar( double omega0 ) {
  if ( _lbar < 0 )
  {
    if ( _idSF == 1 )
    { // exponential SF
      double rat = omega0 * _b / lambda_SF();
      _lbar      = lambda_SF() / _b * ( Gamma( 1 + _b ) - Gamma( 1 + _b, rat ) ) /
              ( Gamma( _b ) - Gamma( _b, rat ) );
    }
    else if ( _idSF == 2 )
    { // Gaussian SF
      double c = cGaus( _b );
      _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 );
    }
  }
  return _lbar;
}
 
double EvtVubNLO::mu_pi2( double omega0 ) {
  if ( _mupi2 < 0 )
  {
    if ( _idSF == 1 )
    { // exponential SF
      double rat = omega0 * _b / lambda_SF();
      _mupi2 = 3 * ( pow( lambda_SF() / _b, 2 ) * ( Gamma( 2 + _b ) - Gamma( 2 + _b, rat ) ) /
                         ( Gamma( _b ) - Gamma( _b, rat ) ) -
                     pow( lambda_bar( omega0 ), 2 ) );
    }
    else if ( _idSF == 2 )
    { // Gaussian SF
      double c  = cGaus( _b );
      double m1 = Gamma( ( 3 + _b ) / 2 ) -
                  Gamma( ( 3 + _b ) / 2, pow( omega0 / lambda_SF(), 2 ) * c );
      double m2 =
          Gamma( 1 + _b / 2 ) - Gamma( 1 + _b / 2, pow( omega0 / lambda_SF(), 2 ) * c );
      double m3 = Gamma( ( 1 + _b ) / 2 ) -
                  Gamma( ( 1 + _b ) / 2, pow( omega0 / lambda_SF(), 2 ) * c );
      _mupi2 = 3 * pow( lambda_SF(), 2 ) * ( m1 / m3 - pow( m2 / m3, 2 ) ) / c;
    }
  }
  return _mupi2;
}
 
double EvtVubNLO::M0( double mui, double omega0 ) {
  double mf = omega0 - lambda_bar( omega0 );
  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 ) );
}
 
double EvtVubNLO::alphas( double mu ) {
  double Lambda4 = 0.302932;
  double lg      = 2 * log( mu / Lambda4 );
  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 ) ) );
}
 
double EvtVubNLO::gausShapeFunction( double omega, const std::vector<double>& sCoeffs ) {
  double b  = sCoeffs[3];
  double l  = sCoeffs[7];
  double wL = omega / l;
 
  return pow( wL, b ) * exp( -cGaus( b ) * wL * wL );
}
 
double EvtVubNLO::expShapeFunction( double omega, const std::vector<double>& sCoeffs ) {
  double b  = sCoeffs[3];
  double l  = sCoeffs[7];
  double wL = omega / l;
 
  return pow( wL, b - 1 ) * exp( -b * wL );
}
 
double EvtVubNLO::Gamma( double z ) {
 
  std::vector<double> gammaCoeffs( 6 );
  gammaCoeffs[0] = 76.18009172947146;
  gammaCoeffs[1] = -86.50532032941677;
  gammaCoeffs[2] = 24.01409824083091;
  gammaCoeffs[3] = -1.231739572450155;
  gammaCoeffs[4] = 0.1208650973866179e-2;
  gammaCoeffs[5] = -0.5395239384953e-5;
 
  // Lifted from Numerical Recipies in C
  double x, y, tmp, ser;
 
  int j;
  y = z;
  x = z;
 
  tmp = x + 5.5;
  tmp = tmp - ( x + 0.5 ) * log( tmp );
  ser = 1.000000000190015;
 
  for ( j = 0; j < 6; j++ )
  {
    y   = y + 1.0;
    ser = ser + gammaCoeffs[j] / y;
  }
 
  return exp( -tmp + log( 2.5066282746310005 * ser / x ) );
}
 
double EvtVubNLO::Gamma( double z, double tmin ) {
  std::vector<double> c( 1 );
  c[0]                       = z;
  EvtItgPtrFunction*   func  = new EvtItgPtrFunction( &dgamma, tmin, 100., c );
  EvtItgAbsIntegrator* jetSF = new EvtItgSimpsonIntegrator( *func, 0.001 );
  return jetSF->evaluate( tmin, 100. );
}