FFBLPRBase.cc

Go to the documentation of this file.

///
/// @file  FFBLPRBase.cc
/// @brief Hammer base class for BLPR form factors
///

//**** This file is a part of the HAMMER library
//**** Copyright (C) 2016 - 2020 The HAMMER Collaboration
//**** HAMMER is licensed under version 3 of the GPL; see COPYING for details
//**** Please note the MCnet academic guidelines; see GUIDELINES for details

// -*- C++ -*-
#include "Hammer/FormFactors/FFBLPRBase.hh"
#include "Hammer/Math/Utils.hh"
#include "Hammer/Math/Constants.hh"
#include "Hammer/Particle.hh"

using namespace std;

namespace Hammer {

    FFBLPRBase::FFBLPRBase() {
        setGroup("BLPR");
    }

    void FFBLPRBase::eval(const Particle& parent, const ParticleList& daughters,
                          const ParticleList&) {
        // Momenta
        const FourMomentum& pBmes = parent.momentum();
        const FourMomentum& pDstarmes = daughters[0].momentum();
        // const FourMomentum& pTau = daughters[2].momentum();


        // kinematic objects
        const double Mb = pBmes.mass();
        const double Mc = pDstarmes.mass();
        // const double Mt = pTau.mass();
        const double Sqq = pow(Mb, 2.) + pow(Mc, 2.) - 2. * (pBmes * pDstarmes);

        evalAtPSPoint({Sqq}, {Mb, Mc});
    }

    void FFBLPRBase::initVars(double w, double z) const {
        _wSq = w*w;
        _sqrt1wSq = sqrt(_wSq - 1);
        _zSq = z*z;
        _zCu = z*z*z;
        _zp1Sq = (z+1.)*(z+1.);
        _zm1Sq = (z-1.)*(z-1.);;
        _lnz = log(z);
        _wZ = (z + 1./z)/2.;
        _wP = w + _sqrt1wSq;
        _wM = w - _sqrt1wSq;
        _wmwZSq = (w - _wZ)*(w - _wZ);
        _lnwP = log(_wP);
        _rW = _lnwP/_sqrt1wSq;
        _Omega = w/(2.*_sqrt1wSq)*( 2*DiLog(1 - _wM*z) - 2*DiLog(1 - _wP*z) + DiLog(1 - _wP*_wP) - DiLog(1 - _wM*_wM) ) - w*_rW*_lnz + 1.;
    }

    double FFBLPRBase::CS(double w, double z) const {
        if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (2.*_Omega*(w - _wZ)*z + (-1. + _zSq)*_lnz - (-1. + w)*_zp1Sq*_rW)/(3.*(w - _wZ)*z);
    }

    double FFBLPRBase::CP(double w, double z) const {
        if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (2.*_Omega*(w - _wZ)*z + (-1. + _zSq)*_lnz - (1. + w)*_zm1Sq*_rW)/(3.*(w - _wZ)*z);
    }

    double FFBLPRBase::CV1(double w, double z) const {
        if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (4.*_Omega*(w - _wZ)*z + 12*(-w + _wZ)*z + _lnz - _zSq*_lnz + 2.*(1. + w)*(-1. + (-1. + 3.*w - z)*z)*_rW)/(6.*(w - _wZ)*z);
    }

    double FFBLPRBase::CV2(double w, double z) const {
        if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (-(z*(-2.*(w - _wZ)*(-1. + z) + (3. - 2.*w + (2. - 4.*w)*z + _zSq)*_lnz)) + (-2. + (-1. + 5.*w + 2.*_wSq)*z - 2.*w*(1. + 2.*w)*_zSq + (1. + w)*_zCu)*_rW)/(6.*_wmwZSq*_zSq);
    }

    double FFBLPRBase::CV3(double w, double z) const {
        if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (-2.*(w - _wZ)*(-1. + z)*z + (1. + (2. - 4.*w)*z + (3. - 2.*w)*_zSq)*_lnz + (1. + 2.*_wSq*(-2. + z)*z - _zSq - 2.*_zCu + w*(1. - 2.*z + 5.*_zSq))*_rW)/(6.*_wmwZSq*z);
    }

    double FFBLPRBase::CA1(double w, double z) const {
        if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (-12*w*z + 4.*_Omega*(w - _wZ)*z + 12*_wZ*z + _lnz - _zSq*_lnz + 2.*(-1. + w)*(-1. + z + 3.*w*z - _zSq)*_rW)/(6.*(w - _wZ)*z);
    }

    double FFBLPRBase::CA2(double w, double z) const {
       if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (-(z*(-2.*(w - _wZ)*(1. + z) + (3. + 2.*w + (-2. - 4.*w)*z + _zSq)*_lnz)) + (-2. + (1. + 5.*w - 2.*_wSq)*z + 2.*(1. - 2.*w)*w*_zSq + (-1. + w)*_zCu)*_rW)/(6.*_wmwZSq*_zSq);
    }

    double FFBLPRBase::CA3(double w, double z) const {
        if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (-2.*(w - _wZ)*z*(1. + z) + (-1. + (2. + 4.*w)*z + (-3. - 2.*w)*_zSq)*_lnz + (1. - _zSq + 2.*_zCu + 2.*_wSq*z*(2. + z) - w*(1. + 2.*z + 5.*_zSq))*_rW)/(6.*_wmwZSq*z);
    }

    double FFBLPRBase::CT1(double w, double z) const {
        if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (-6.*w*z + 2.*_Omega*(w - _wZ)*z + 6.*_wZ*z + _lnz - _zSq*_lnz + (-1. + w)*(-1. + (2. + 4.*w)*z - _zSq)*_rW)/(3.*(w - _wZ)*z);
    }

    double FFBLPRBase::CT2(double w, double z) const {
        if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (-2.*((1. - w*z)*_rW + z*_lnz))/(3.*(-w + _wZ)*z);
    }

    double FFBLPRBase::CT3(double w, double z) const {
       if(isZero(w - 1.)) w += 1e-6;
        initVars(w , z);
        return (-2.*((w - z)*_rW + _lnz))/(3.*(-w + _wZ));
}

    double FFBLPRBase::DiLog(double z) const {
// ====================================================================
// This file is part of Dilogarithm.
//
// Dilogarithm is licenced under the GNU Lesser General Public License
// (GNU LGPL) version 3.
// Code translated by R.Brun from CERNLIB DILOG function C332
// ====================================================================
   constexpr double HF  = 0.5;
   constexpr double PI3 = pi2/3.;
   constexpr double PI6 = pi2/6.;
   constexpr double PI12 = pi2/12.;
   constexpr double C[20] = {0.42996693560813697, 0.40975987533077105,
     -0.01858843665014592, 0.00145751084062268,-0.00014304184442340,
      0.00001588415541880,-0.00000190784959387, 0.00000024195180854,
     -0.00000003193341274, 0.00000000434545063,-0.00000000060578480,
      0.00000000008612098,-0.00000000001244332, 0.00000000000182256,
     -0.00000000000027007, 0.00000000000004042,-0.00000000000000610,
      0.00000000000000093,-0.00000000000000014, 0.00000000000000002};

   double T,H,Y,S,A,ALFA,B1,B2,B0;

   if (isZero(z - 1.)) {
       H = PI6;
   } else if (isZero(z + 1.)) {
       H = -PI12;
   } else {
       T = -z;
       if (T <= -2) {
           Y = -1/(1+T);
           S = 1;
           B1= log(-T);
           B2= log(1+1/T);
           A = -PI3+HF*(B1*B1-B2*B2);
       } else if (T < -1) {
           Y = -1-T;
           S = -1;
           A = log(-T);
           A = -PI6+A*(A+log(1+1/T));
       } else if (T <= -0.5) {
           Y = -(1+T)/T;
           S = 1;
           A = log(-T);
           A = -PI6+A*(-HF*A+log(1+T));
       } else if (T < 0) {
           Y = -T/(1+T);
           S = -1;
           B1= log(1+T);
           A = HF*B1*B1;
       } else if (T <= 1) {
           Y = T;
           S = 1;
           A = 0;
       } else {
           Y = 1/T;
           S = -1;
           B1= log(T);
           A = PI6+HF*B1*B1;
       }
       H    = Y+Y-1;
       ALFA = H+H;
       B1   = 0;
       B2   = 0;
       B0   = 0;
       for (int i=19;i>=0;i--){
          B0 = C[i] + ALFA*B1-B2;
          B2 = B1;
          B1 = B0;
       }
       H = -(S*(B0-H*B2)+A);
    }
        return H;
    }

    void FFBLPRBase::addRefs() const {
        if(!getSettingsHandler()->checkReference("Bernlochner:2017jka")){
            string ref =
                "@article{Bernlochner:2017jka,\n"
                "      author         = \"Bernlochner, Florian U. and Ligeti, Zoltan and Papucci, Michele and Robinson, Dean J.\",\n"
                "      title          = \"{Combined analysis of semileptonic $B$ decays to $D$ and $D^*$: $R(D^{(*)})$, $|V_{cb}|$, and new physics}\",\n"
                "      journal        = \"Phys. Rev.\",\n"
                "      volume         = \"D95\",\n"
                "      year           = \"2017\",\n"
                "      number         = \"11\",\n"
                "      pages          = \"115008\",\n"
                "      doi            = \"10.1103/PhysRevD.95.115008, 10.1103/PhysRevD.97.059902\",\n"
                "      note           = \"[erratum: Phys. Rev.D97,no.5,059902(2018)]\",\n"
                "      eprint         = \"1703.05330\",\n"
                "      archivePrefix  = \"arXiv\",\n"
                "      primaryClass   = \"hep-ph\",\n"
                "      SLACcitation   = \"%%CITATION = ARXIV:1703.05330;%%\"\n"
                "}\n";
            getSettingsHandler()->addReference("Bernlochner:2017jka", ref);
        }
    }
}