FFLbtoLcBLRSVar.cc

Go to the documentation of this file.

///
/// @file  FFLbtoLcBLRSVar.cc
/// @brief \f$ \Lambda_b \rightarrow Lambda_c \f$ BLRS 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/FFLbtoLcBLRSVar.hh"
#include "Hammer/IndexLabels.hh"
#include "Hammer/Math/Constants.hh"
#include "Hammer/Math/Utils.hh"
#include "Hammer/Particle.hh"
#include "Hammer/Tools/Pdg.hh"
#include "Hammer/Math/MultiDim/SparseContainer.hh"
#include <cmath>

#include <iostream>

using namespace std;

namespace Hammer {

    namespace MD = MultiDimensional;

    FFLbtoLcBLRSVar::FFLbtoLcBLRSVar() {
        // Create tensor rank and dimensions
        vector<IndexType> dims = {12,5}; // size is _FFErrNames + 1 to include central values in zeroth component
        setGroup("BLRSVar"); //override BLRS base class FF group setting
        string name{"FFLbtoLcBLRSVar"};
        
        setPrefix("LbtoLc");
        _FFErrLabel = FF_LBLC_VAR;
        addProcessSignature(-PID::LAMBDAB, {PID::LAMBDACMINUS});
        addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_LBLC, FF_LBLC_VAR})});

        setSignatureIndex();
    }

    void FFLbtoLcBLRSVar::defineSettings() {
        _FFErrNames = {"delta_z1", "delta_z2", "delta_b1", "delta_b2"};
        setPath(getFFErrPrefixGroup().get());
        setUnits("GeV");

        //1S scheme: mb1S = 4721.5, delta mb-mc = 3400., alpha_s = 26/100
        addSetting<double>("as",0.26);
        addSetting<double>("mb1S",4.7215); //GeV
        addSetting<double>("dmbmc",3.400);  //GeV
        addSetting<double>("mLb",5.620); //GeV
        addSetting<double>("mLc",2.286);  //GeV

        addSetting<double>("z1",-2.037);
        addSetting<double>("z2",3.159);
        addSetting<double>("b1",-0.459); //GeV^2
        addSetting<double>("b2",-0.390); //GeV^2

        //correlation matrix and set zero error eigenvectors
        //Row basis is z1, z2, b1 ,b2
        vector<vector<double>> sliwmat{ {1., 0., 0., 0.},
                                        {0., 1., 0., 0.},
                                        {0., 0., 1., 0.},
                                        {0., 0., 0., 1.}};
        addSetting<vector<vector<double>>>("sliwmatrix",sliwmat);
        
        for (auto elem : _FFErrNames) {
            addSetting<double>(elem, 0.);
        }
        
        initialized = true;
    }

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

        // kinematic objects
        const double Mb = pLbmes.mass();
        const double Mc = pLcmes.mass();
        // const double Mt = pTau.mass();
        const double Sqq = Mb*Mb + Mc*Mc - 2. * (pLbmes * pLcmes);

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

    void FFLbtoLcBLRSVar::evalAtPSPoint(const vector<double>& point, const vector<double>& masses) {
        Tensor& result = getTensor();
        result.clearData();

        if(!initialized){
            MSG_WARNING("Warning, Settings have not been defined!");
        }

        double Mb = 0.;
        double Mc = 0.;
        double unitres = 1.;
        if(masses.size() >= 2) {
            Mb = masses[0];
            Mc = masses[1];
            unitres = _units;
        }
        else {
            Mb = this->masses()[0];
            Mc = this->masses()[1];
        }
        const double Mb2 = Mb*Mb;
        const double Mc2 = Mc*Mc;

        double Sqq = point[0];

        // const double sqSqq = sqrt(Sqq);
        double w = (Mb2 + Mc2 - Sqq)/(2.*Mb*Mc);
        //safety measure if w==1.0
        if(isZero(w - 1.0)) w += 1e-6;
        const double wSq = w*w;
        const double wm1sq = (w-1.)*(w-1.);

        //Renormalon improvement scheme
        const double aS = (*getSetting<double>("as"))/pi;
        const double mb1S = unitres*(*getSetting<double>("mb1S"));
        const double dmbmc = unitres*(*getSetting<double>("dmbmc"));
        const double mLb = unitres*(*getSetting<double>("mLb"));
        const double mLc = unitres*(*getSetting<double>("mLc"));
        const double mc1S = mb1S - dmbmc;
        const double mc2 = pow(mc1S,2.);

//        const double mb = mb1S*(1. + 2.*pow(aS*pi,2.)/9.);
//        const double mc = mb - dmbmc;
//        const double la = (mb1S*(mLb - mb) - mc1S*(mLc - mc))/(mb1S - mc1S);
        const double la1S = (mb1S*(mLb - mb1S) - mc1S*(mLc - mc1S))/(mb1S - mc1S);

        const double eB0 = la1S/(2.*mb1S);
        const double eC0 = la1S/(2.*mc1S);
//        const double eB = la/(2.*mb);
//        const double eC = la/(2.*mc);

        //Explicit expansion to O(eps)
        const double eps = 2.*pow(aS*pi,2.)/9.;
        const double eB = eB0 + eps*(dmbmc*mc1S - mb1S*mLb + mc1S*mLc)/(mb1S - mc1S)/(2.*mb1S);
        const double eC = eC0 + eps*((mb1S*(dmbmc*mb1S - mb1S*mLb + mc1S*mLc))/(mc1S*(mb1S - mc1S)))/(2.*mc1S);

        const double zBC = mc1S/mb1S;

        //(Sub)leading IW function parameters
        const double z1 = (*getSetting<double>("z1"));
        const double z2 = (*getSetting<double>("z2"));
        const double b1 = unitres*unitres*(*getSetting<double>("b1"));
        const double b2 = unitres*unitres*(*getSetting<double>("b2"));

        const vector<vector<double>>& sliwmat = (*getSetting<vector<vector<double>>>("sliwmatrix"));
        
        // Leading and sub-leading IW function
        const double zetaIW = 1. + z1*(w-1) + 0.5*z2*pow(w-1,2.);

        //QCD correction functions
        const double Cs  = CS(w, zBC);
        const double Cps = CP(w, zBC);
        const double Cv1 = CV1(w, zBC);
        const double Cv2 = CV2(w, zBC);
        const double Cv3 = CV3(w, zBC);
        const double Ca1 = CA1(w, zBC);
        const double Ca2 = CA2(w, zBC);
        const double Ca3 = CA3(w, zBC);
        const double Ct1 = CT1(w, zBC);
        const double Ct2 = CT2(w, zBC);
        const double Ct3 = CT3(w, zBC);
        //Derivatives (approx)
        const double Csp  =  (CS(w + 1e-6, zBC) - Cs)/1e-6;
        const double Cpsp =  (CP(w + 1e-6, zBC) - Cps)/1e-6;
        const double Cv1p =  (CV1(w + 1e-6, zBC) - Cv1)/1e-6;
        const double Cv2p =  (CV2(w + 1e-6, zBC) - Cv2)/1e-6;
        const double Cv3p =  (CV3(w + 1e-6, zBC) - Cv3)/1e-6;
        const double Ca1p =  (CA1(w + 1e-6, zBC) - Ca1)/1e-6;
        const double Ca2p =  (CA2(w + 1e-6, zBC) - Ca2)/1e-6;
        const double Ca3p =  (CA3(w + 1e-6, zBC) - Ca3)/1e-6;
        const double Ct1p =  (CT1(w + 1e-6, zBC) - Ct1)/1e-6;
        const double Ct2p =  (CT2(w + 1e-6, zBC) - Ct2)/1e-6;
        const double Ct3p =  (CT3(w + 1e-6, zBC) - Ct3)/1e-6;

        // Create hatted h FFs
        const double hs = 1 + b1/(4.*mc2) + (eB*(w - 1))/(1 + w) + (eC*(w - 1))/(1 + w) +
                    aS*(Cs + (eB0*(w - 1)*(Cs + 2.*Csp*(1 + w)))/(1 + w) + (eC0*(w - 1)*(Cs + 2.*Csp*(1 + w)))/(1 + w)) ;
        const double hp = 1 + eB + eC + b1/(4.*mc2) - b2/(4.*mc2) + aS*(Cps + eB0*(Cps + 2.*Cpsp*(w - 1)) + eC0*(Cps + 2.*Cpsp*(w - 1))) ;
        const double f1 = 1 + eB + eC + b1/(4.*mc2) - b2/(4.*mc2) + aS*(Cv1 + eB0*(Cv1 + 2.*Cv1p*(w - 1)) + eC0*(Cv1 + 2.*Cv1p*(w - 1))) ;
        const double f2 = b2/(4.*mc2) - (2.*eC)/(1 + w) + aS*(Cv2 + (eC0*(-2.*Cv1 - Cv2 - 2.*Cv3 + Cv2*w
                                                                            + 2.*Cv2p*(-1 + wSq)))/(1 + w) +
                    (eB0*(Cv2*(-1 + 3.*w) + 2.*Cv2p*(-1 + wSq)))/(1 + w)) ;
        const double f3 = (-2.*eB)/(1 + w) + aS*(Cv3 + (eB0*(-2.*Cv1 - 2.*Cv2 - Cv3 + Cv3*w
                                                                            + 2.*Cv3p*(-1 + wSq)))/(1 + w) +
                    (eC0*(Cv3*(-1 + 3.*w) + 2.*Cv3p*(-1 + wSq)))/(1 + w)) ;
        const double g1 = 1 + b1/(4.*mc2) + (eB*(w - 1))/(1 + w) + (eC*(w - 1))/(1 + w) +
                    aS*(Ca1 + (eB0*(w - 1)*(Ca1 + 2.*Ca1p*(1 + w)))/(1 + w) + (eC0*(w - 1)*(Ca1 + 2.*Ca1p*(1 + w)))/(1 + w)) ;
        const double g2 = b2/(4.*mc2) - (2.*eC)/(1 + w) + aS*(Ca2 + (eC0*(-2.*Ca1 + Ca2 - 2.*Ca3 + Ca2*w
                                                                            + 2.*Ca2p*(-1 + wSq)))/(1 + w) +
                    (eB0*(Ca2 + 3.*Ca2*w + 2.*Ca2p*(-1 + wSq)))/(1 + w)) ;
        const double g3 = (2.*eB)/(1 + w) + aS*(Ca3 + (eB0*(2.*Ca1 - 2.*Ca2 + Ca3 + Ca3*w
                                                                            + 2.*Ca3p*(-1 + wSq)))/(1 + w) +
                    (eC0*(Ca3 + 3.*Ca3*w + 2.*Ca3p*(-1 + wSq)))/(1 + w)) ;
        const double h1 = 1 + b1/(4.*mc2) + (eB*(w - 1))/(1 + w) + (eC*(w - 1))/(1 + w) +
                    aS*(Ct1 + (eB0*(w - 1)*(Ct1 + 2.*Ct1p*(1 + w)))/(1 + w) + (eC0*(w - 1)*(Ct1 + 2.*Ct1p*(1 + w)))/(1 + w)) ;
        const double h2 = b2/(4.*mc2) - (2.*eC)/(1 + w) + aS*(Ct2 + (eC0*(-2.*Ct1 + Ct2 - 2.*Ct3 + Ct2*w
                                                                            + 2.*Ct2p*(-1 + wSq)))/(1 + w) +
                    (eB0*(Ct2 + 3.*Ct2*w + 2.*Ct2p*(-1 + wSq)))/(1 + w)) ;
        const double h3 = (2.*eB)/(1 + w) + aS*(Ct3 + (eB0*(2.*Ct1 - 2.*Ct2 + Ct3 + Ct3*w
                                                                            + 2.*Ct3p*(-1 + wSq)))/(1 + w) +
                    (eC0*(Ct3 + 3.*Ct3*w + 2.*Ct3p*(-1 + wSq)))/(1 + w)) ;
        const double h4 = aS*((-2.*Ct2*eB0)/(1 + w) + (2.*Ct3*eC0)/(1 + w)) ;

        // set elements
        result.element({0, 0}) = hs;
        result.element({1, 0}) = hp;
        result.element({2, 0}) = f1;
        result.element({3, 0}) = f2;
        result.element({4, 0}) = f3;
        result.element({5, 0}) = g1;
        result.element({6, 0}) = g2;
        result.element({7, 0}) = g3;
        result.element({8, 0}) = h1;
        result.element({9, 0}) = h2;
        result.element({10, 0}) = h3;
        result.element({11, 0}) = h4;

        //Now create remaining FF tensor entries
        //n indexes rows ie FFs; idx indexes the columns that contract with delta
        
        //Linearization of hatted FFs
        vector<vector<double>> FFmat{{(hs*(-1 + w))/zetaIW, (hs*wm1sq)/(2.*zetaIW), 1./(4.*mc2), 0.},
                                     {(hp*(-1 + w))/zetaIW, (hp*wm1sq)/(2.*zetaIW), 1./(4.*mc2), -1./(4.*mc2)},
                                     {(f1*(-1 + w))/zetaIW, (f1*wm1sq)/(2.*zetaIW), 1./(4.*mc2), -1./(4.*mc2)},
                                     {(f2*(-1 + w))/zetaIW, (f2*wm1sq)/(2.*zetaIW), 0., 1./(4.*mc2)},
                                     {(f3*(-1 + w))/zetaIW, (f3*wm1sq)/(2.*zetaIW), 0., 0.},
                                     {(g1*(-1 + w))/zetaIW, (g1*wm1sq)/(2.*zetaIW), 1./(4.*mc2), 0.},
                                     {(g2*(-1 + w))/zetaIW, (g2*wm1sq)/(2.*zetaIW), 0., 1./(4.*mc2)},
                                     {(g3*(-1 + w))/zetaIW, (g3*wm1sq)/(2.*zetaIW), 0., 0.},
                                     {(h1*(-1 + w))/zetaIW, (h1*wm1sq)/(2.*zetaIW), 1./(4.*mc2), 0.},
                                     {(h2*(-1 + w))/zetaIW, (h2*wm1sq)/(2.*zetaIW), 0., 1./(4.*mc2)},
                                     {(h3*(-1 + w))/zetaIW, (h3*wm1sq)/(2.*zetaIW), 0., 0.},
                                     {(h4*(-1 + w))/zetaIW, (h4*wm1sq)/(2.*zetaIW), 0., 0.}};
        
        for(size_t n = 0; n < 12; ++n){
            auto ni = static_cast<IndexType>(n);
            for(size_t idx = 0; idx < 4; ++idx){
                complex<double> entry = 0;
                auto idxi = static_cast<IndexType>(idx+1u);
                for(size_t idx1 = 0; idx1 < 4; ++idx1){
                    entry += sliwmat[idx1][idx]*FFmat[n][idx1];
                }
                result.element({ni, idxi}) = entry;
            }
        }
        
        
        result *= zetaIW;

    }

    std::unique_ptr<FormFactorBase> FFLbtoLcBLRSVar::clone(const std::string& label) {
        MAKE_CLONE(FFLbtoLcBLRSVar, label);
    }

} // namespace Hammer