include/EvtGenModels/EvtVubNLO.hh

0001
0002 /***********************************************************************
0003 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
0004 *                                                                      *
0005 * This file is part of EvtGen.                                         *
0006 *                                                                      *
0007 * EvtGen is free software: you can redistribute it and/or modify       *
0008 * it under the terms of the GNU General Public License as published by *
0009 * the Free Software Foundation, either version 3 of the License, or    *
0010 * (at your option) any later version.                                  *
0011 *                                                                      *
0012 * EvtGen is distributed in the hope that it will be useful,            *
0013 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
0014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
0015 * GNU General Public License for more details.                         *
0016 *                                                                      *
0017 * You should have received a copy of the GNU General Public License    *
0018 * along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
0019 ***********************************************************************/
0020
0021 #ifndef EVTVUBNLO_HH
0022 #define EVTVUBNLO_HH
0023
0024 #include "EvtGenBase/EvtDecayIncoherent.hh"
0025
0026 #include <vector>
0027
0028 class EvtParticle;
0029 class RandGeneral;
0030
0031 // Description:
0032 // Class to generate inclusive B to X_u l nu decays according to various
0033 // decay models. Implemtented are ACCM, parton-model and a QCD model.
0034 // Description: Routine to decay B->Xulnu according to Bosch, Lange, Neubert, and Paz hep-ph/0402094
0035 //              Equation numbers refer to this paper
0036
0037 class EvtVubNLO : public EvtDecayIncoherent {
0038   public:
0039     EvtVubNLO() = default;
0040     ~EvtVubNLO();
0041
0042     std::string getName() override;
0043
0044     EvtDecayBase* clone() override;
0045
0046     void initProbMax() override;
0047
0048     void init() override;
0049
0050     void decay( EvtParticle* p ) override;
0051
0052   private:
0053     // cache
0054     double _lbar;
0055     double _mupi2;
0056
0057     double _mb;    // the b-quark pole mass in GeV
0058     double _mB;
0059     double _lambdaSF;
0060     double _b;    // Parameter for the Fermi Motion
0061     double _kpar;
0062     double _mui;       // renormalization scale (preferred value=1.5 GeV)
0063     double _SFNorm;    // SF normalization
0064     double _dGMax;     // max dGamma*p2 value;
0065     int _idSF;         // which shape function?
0066     std::vector<double> _masses;
0067     std::vector<double> _weights;
0068
0069     double _gmax;
0070     int _ngood, _ntot;
0071
0072     double tripleDiff( double pp, double pl, double pm );
0073     double SFNorm( const std::vector<double>& coeffs );
0074     static double integrand( double omega, const std::vector<double>& coeffs );
0075     double F10( const std::vector<double>& coeffs );
0076     static double F1Int( double omega, const std::vector<double>& coeffs );
0077     double F20( const std::vector<double>& coeffs );
0078     static double F2Int( double omega, const std::vector<double>& coeffs );
0079     double F30( const std::vector<double>& coeffs );
0080     static double F3Int( double omega, const std::vector<double>& coeffs );
0081     static double g1( double y, double z );
0082     static double g2( double y, double z );
0083     static double g3( double y, double z );
0084
0085     static double Gamma( double z );    // Euler Gamma Function
0086     static double dgamma( double t, const std::vector<double>& c )
0087     {
0088         return pow( t, c[0] - 1 ) * exp( -t );
0089     }
0090     static double Gamma( double z, double tmax );
0091
0092     // theory parameters
0093     inline double mu_i() { return _mui; }    // intermediate scale
0094     inline double mu_bar() { return _mui; }
0095     inline double mu_h() { return _mb / sqrt( 2.0 ); }    // high scale
0096     inline double lambda1() { return -_mupi2; }
0097
0098     // expansion coefficients for RGE
0099     static double beta0( int nf = 4 ) { return 11. - 2. / 3. * nf; }
0100     static double beta1( int nf = 4 ) { return 34. * 3. - 38. / 3. * nf; }
0101     static double beta2( int nf = 4 )
0102     {
0103         return 1428.5 - 5033. / 18. * nf + 325. / 54. * nf * nf;
0104     }
0105     static double gamma0() { return 16. / 3.; }
0106     static double gamma1( int nf = 4 )
0107     {
0108         return 4. / 3. * ( 49.85498 - 40. / 9. * nf );
0109     }
0110     static double gamma2( int nf = 4 )
0111     {
0112         return 64. / 3. * ( 55.07242 - 8.58691 * nf - nf * nf / 27. );
0113     } /*  zeta3=1.20206 */
0114     static double gammap0() { return -20. / 3.; }
0115     static double gammap1( int nf = 4 )
0116     {
0117         return -32. / 3. * ( 6.92653 - 0.9899 * nf );
0118     } /* ??  zeta3=1.202 */
0119
0120     // running constants
0121
0122     static double alphas( double mu );
0123     static double C_F( double mu )
0124     {
0125         return ( 4.0 / 3.0 ) * alphas( mu ) / 4. / EvtConst::pi;
0126     }
0127
0128     // Shape Functions
0129
0130     inline double lambda_SF() { return _lambdaSF; }
0131     double lambda_bar( double omega0 );
0132     inline double lambda2() { return 0.12; }
0133     double mu_pi2( double omega0 );
0134     inline double lambda( double ) { return _mB - _mb; }
0135
0136     // specail for gaussian SF
0137     static double cGaus( double b )
0138     {
0139         return pow( Gamma( 1 + b / 2. ) / Gamma( ( 1 + b ) / 2. ), 2 );
0140     }
0141
0142     double M0( double mui, double omega0 );
0143     static double shapeFunction( double omega, const std::vector<double>& coeffs );
0144     static double expShapeFunction( double omega,
0145                                     const std::vector<double>& coeffs );
0146     static double gausShapeFunction( double omega,
0147                                      const std::vector<double>& coeffs );
0148     // SSF (not yet implemented)
0149     double subS( const std::vector<double>& coeffs );
0150     double subT( const std::vector<double>& coeffs );
0151     double subU( const std::vector<double>& coeffs );
0152     double subV( const std::vector<double>& coeffs );
0153
0154     // Sudakov
0155
0156     inline double S0( double a, double r )
0157     {
0158         return -gamma0() / 4 / a / pow( beta0(), 2 ) * ( 1 / r - 1 + log( r ) );
0159     }
0160     inline double S1( double /*a*/, double r )
0161     {
0162         return gamma0() / 4. / pow( beta0(), 2 ) *
0163                ( pow( log( r ), 2 ) * beta1() / 2. / beta0() +
0164                  ( gamma1() / gamma0() - beta1() / beta0() ) *
0165                      ( 1. - r + log( r ) ) );
0166     }
0167     inline double S2( double a, double r )
0168     {
0169         return gamma0() * a / 4. / pow( beta0(), 2 ) *
0170                ( -0.5 * pow( ( 1 - r ), 2 ) *
0171                      ( pow( beta1() / beta0(), 2 ) - beta2() / beta0() -
0172                        beta1() / beta0() * gamma1() / gamma0() +
0173                        gamma2() / gamma0() ) +
0174                  ( pow( beta1() / beta0(), 2 ) - beta2() / beta0() ) *
0175                      ( 1 - r ) * log( r ) +
0176                  ( beta1() / beta0() * gamma1() / gamma0() - beta2() / beta0() ) *
0177                      ( 1 - r + r * log( r ) ) );
0178     }
0179     inline double dSudakovdepsi( double mu1, double mu2 )
0180     {
0181         return S2( alphas( mu1 ) / ( 4 * EvtConst::pi ),
0182                    alphas( mu2 ) / alphas( mu1 ) );
0183     }
0184     inline double Sudakov( double mu1, double mu2, double epsi = 0 )
0185     {
0186         double fp( 4 * EvtConst::pi );
0187         return S0( alphas( mu1 ) / fp, alphas( mu2 ) / alphas( mu1 ) ) +
0188                S1( alphas( mu1 ) / fp, alphas( mu2 ) / alphas( mu1 ) ) +
0189                epsi * dSudakovdepsi( mu1, mu2 );
0190     }
0191
0192     // RG
0193     inline double dGdepsi( double mu1, double mu2 )
0194     {
0195         return 1. / 8. / EvtConst::pi * ( alphas( mu2 ) - alphas( mu1 ) ) *
0196                ( gamma1() / beta0() - beta1() * gamma0() / pow( beta0(), 2 ) );
0197     }
0198     inline double aGamma( double mu1, double mu2, double epsi = 0 )
0199     {
0200         return gamma0() / 2 / beta0() * log( alphas( mu2 ) / alphas( mu1 ) ) +
0201                epsi * dGdepsi( mu1, mu2 );
0202     }
0203     inline double dgpdepsi( double mu1, double mu2 )
0204     {
0205         return 1. / 8. / EvtConst::pi * ( alphas( mu2 ) - alphas( mu1 ) ) *
0206                ( gammap1() / beta0() - beta1() * gammap0() / pow( beta0(), 2 ) );
0207     }
0208     inline double agammap( double mu1, double mu2, double epsi = 0 )
0209     {
0210         return gammap0() / 2 / beta0() * log( alphas( mu2 ) / alphas( mu1 ) ) +
0211                epsi * dgpdepsi( mu1, mu2 );
0212     }
0213     inline double U1( double mu1, double mu2, double epsi = 0 )
0214     {
0215         return exp( 2 * ( Sudakov( mu1, mu2, epsi ) - agammap( mu1, mu2, epsi ) -
0216                           aGamma( mu1, mu2, epsi ) * log( _mb / mu1 ) ) );
0217     }
0218     inline double U1lo( double mu1, double mu2 ) { return U1( mu1, mu2 ); }
0219     inline double U1nlo( double mu1, double mu2 )
0220     {
0221         return U1( mu1, mu2 ) *
0222                ( 1 + 2 * ( dSudakovdepsi( mu1, mu2 ) - dgpdepsi( mu1, mu2 ) -
0223                            log( _mb / mu1 ) * dGdepsi( mu1, mu2 ) ) );
0224     }
0225     inline double alo( double mu1, double mu2 )
0226     {
0227         return -2 * aGamma( mu1, mu2 );
0228     }
0229     inline double anlo( double mu1, double mu2 )
0230     {
0231         return -2 * dGdepsi( mu1, mu2 );
0232     }
0233 };
0234
0235 #endif