FFBtoDstarBLPR.cc

Go to the documentation of this file.

///
/// @file  FFBtoDstarBLPR.cc
/// @brief \f$ B \rightarrow D^* \f$ 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/FFBtoDstarBLPR.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;

    FFBtoDstarBLPR::FFBtoDstarBLPR() {
        // Create tensor rank and dimensions
        vector<IndexType> dims = {8};
        string name{"FFBtoDstarBLPR"};
        
        setPrefix("BtoD*");
        addProcessSignature(PID::BPLUS, {-PID::DSTAR});
        addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BDSTAR})});

        addProcessSignature(PID::BZERO, {PID::DSTARMINUS});
        addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BDSTAR})});

        setPrefix("BstoDs*");
        addProcessSignature(PID::BS, {PID::DSSTARMINUS});
        addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BSDSSTAR})});
        
        setSignatureIndex();
    }

    void FFBtoDstarBLPR::defineSettings() {
        //_FFErrNames = ;
        setPath(getFFErrPrefixGroup().get());
        setUnits("GeV");

        addSetting<double>("RhoSq",1.24);
        addSetting<double>("a",1.509/sqrt2);
        //1S scheme: mb1S = 4710, mbBar = 5313, lambda1 = -3 10^5 MeV^2, delta mb-mc = 3400, alpha_s = 26/100
        addSetting<double>("as",0.26);
        addSetting<double>("mb",4.710);
        addSetting<double>("mc",4.710 - 3.400);
        addSetting<double>("la",0.57115);
        //Renormalon improvement (1S scheme)
        addSetting<double>("ebR/eb",0.861);
        addSetting<double>("ecR/ec",0.822);

        addSetting<double>("chi1p",0.);
        addSetting<double>("chi21",-0.06);
        addSetting<double>("chi2p",-0.00);
        addSetting<double>("chi3p",0.05);
        addSetting<double>("eta1",0.30);
        addSetting<double>("etap",-0.05);
        addSetting<double>("dV20",0.);

        initialized = true;

    }

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

        if(!initialized){
            MSG_WARNING("Warning, Setting 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 Mb3 = Mb*Mb*Mb;
        const double Mc2 = Mc*Mc;

        double Sqq = point[0];

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

        //Conformal parameters
        const double a = *getSetting<double>("a");
        const double zCon = (sqrt(w+1) - sqrt2*a)/(sqrt(w+1) + sqrt2*a);

        //BLPR expansion parameters
        const double zBC = (*getSetting<double>("mc"))/(*getSetting<double>("mb"));
        const double eb = (*getSetting<double>("la"))/(*getSetting<double>("mb")*2.);
        const double ebReb = (*getSetting<double>("ebR/eb"));
        const double ec = (*getSetting<double>("la"))/(*getSetting<double>("mc")*2.);
        const double ecRec = (*getSetting<double>("ecR/ec"));
        const double as = (*getSetting<double>("as"))/pi;

        //(Sub)leading IW function parameters
        const double RhoSq = *getSetting<double>("RhoSq");
        const double chi21 = (*getSetting<double>("chi21"));
        const double chi2p = (*getSetting<double>("chi2p"));
        const double chi3p = (*getSetting<double>("chi3p"));
        const double eta1 = (*getSetting<double>("eta1"));
        const double etap = (*getSetting<double>("etap"));

        //L and helper variables
        const double L1 = -4.*(w-1)*(chi21 + (w-1.)*chi2p)+12.*chi3p*(w-1.);
        const double L2 = -4.*chi3p*(w-1.);
        const double L3 = 4.*(chi21 + (w-1.)*chi2p);
        // const double L4 = 2.*(eta1 + etap*(w-1.))-1.;
        const double L4b = (2.*(eta1 + etap*(w-1.))-1.*ebReb);
        // const double L4c = (2.*(eta1 + etap*(w-1.))-1.*ecRec);
        // const double L5 = -1.;
        const double L5c = -ecRec;
        // const double L6 = -2.*(1+(eta1 + etap*(w-1.)))/(w+1.);
        const double L6c = -2.*(ecRec + (eta1 + etap*(w-1.)))/(w+1.);

        // Variables for leading IW function (derived from G(1))
        const double rD = 1867./5280.;
        // const double V11 = 8*pow(a,2);
        const double V21 = 57.0;
        const double V20 = 7.5 + (*getSetting<double>("dV20"));

        //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) at w_0 = 2a^2-1. Calculated once and stored.
        if(!initCs){
            const double w0 = 2*a*a - 1.;
            Cv1z =  CV1(w0, zBC);
            Cv2z =  CV2(w0, zBC);
            Cv3z =  CV3(w0, zBC);

            Cv1zp =  (CV1(w0 + 1e-5, zBC) - Cv1z)/1e-5;
            Cv2zp =  (CV2(w0 + 1e-5, zBC) - Cv2z)/1e-5;
            Cv3zp =  (CV3(w0 + 1e-5, zBC) - Cv3z)/1e-5;

            Cv1zpp =  (Cv1zp - (Cv1z - CV1(w0 - 1e-5, zBC))/1e-5)/1e-5;
            Cv2zpp =  (Cv2zp - (Cv2z - CV2(w0 - 1e-5, zBC))/1e-5)/1e-5;
            Cv3zpp =  (Cv3zp - (Cv3z - CV3(w0 - 1e-5, zBC))/1e-5)/1e-5;

            initCs = true;
        }
//        //alpha_s derivatives at w_0 = 2a^2-1
//        // const double Cv1z = 0.97543966986296275;
//        const double Cv2z = -0.416360981425492067;
//        const double Cv3z = -0.186896365253791276;
//
//        const double Cv1zp = 0.11470494247949626;
//        const double Cv2zp = 0.16684760479981954;
//        const double Cv3zp = 0.038382492791834666;
//
//        const double Cv1zpp = -0.2557645997620378;
//        const double Cv2zpp = -0.1297628996251381;
//        const double Cv3zpp = -0.02238541998227300;


        // Leading and sub-leading IW function
        //-------------------------------------------------------------------------------------------
        const double a2 = a*a;
        const double a4 = a2*a2;
        const double a6 = a4*a2;
        const double Xi = 64*a4*RhoSq - 16*a2 - V21;
        double xiIW = 1 - 8*a2*RhoSq*zCon + zCon*zCon*(
                    V21*RhoSq - V20 + (eb - ec)*(2*Xi*etap * (1-rD)/(1+rD))
                    + (eb + ec)*(Xi* (12*chi3p - 4*chi21) - 16*((a2-1)*Xi - 16* a4)*chi2p)
                    + as*(Xi*(Cv1zp +(Cv3z + rD*Cv2z)/(1 + rD)) + 2*a2*(Xi - 32*a2)*(Cv3zp + rD*Cv2zp)/(1+rD) - 64*a6*(Cv3zpp + rD*Cv2zpp)/(1+rD) -32*a4*Cv1zpp ));
        xiIW /= 1 - 8*a2*RhoSq*(1-a)/(1+a) + pow((1-a)/(1+a),2.)*(
                    V21*RhoSq - V20 + (eb - ec)*(2*Xi*etap * (1-rD)/(1+rD))
                    + (eb + ec)*(Xi* (12*chi3p - 4*chi21) - 16*((a2-1)*Xi - 16* a4)*chi2p)
                    + as*(Xi*(Cv1zp +(Cv3z + rD*Cv2z)/(1 + rD)) + 2*a2*(Xi - 32*a2)*(Cv3zp + rD*Cv2zp)/(1+rD) - 64*a6*(Cv3zpp + rD*Cv2zpp)/(1+rD) -32*a4*Cv1zpp ));

        // LO IW function linear only
        //double xiIW = 1. - 8 * pow(a, 2.) * RhoSq * zCon;

        const double chi1=0;
        const double h=xiIW + 2.*(eb+ec)*chi1;

        //Create hatted h FFs
        const double Hps = 1.+as*Cps+ec*(L2+L3*(w-1)+L5c-L6c*(w+1))+eb*(L1-L4b);
        const double Hv = 1.+as*Cv1+ec*(L2-L5c)+eb*(L1-L4b);
        const double Ha1 = 1.+as*Ca1+ec*(L2-L5c*(w-1)/(w+1))+eb*(L1-L4b*(w-1)/(w+1));
        const double Ha2 = as*Ca2+ec*(L3+L6c);
        const double Ha3 = 1.+as*(Ca1+Ca3)+ec*(L2-L3-L5c+L6c)+eb*(L1-L4b);
        const double Ht1 = 1.+as*(Ct1+0.5*(w-1)*(Ct2-Ct3))+ec*(L2)+eb*(L1);
        const double Ht2 = 0.5*(w+1)*as*(Ct2+Ct3)+ec*(L5c)-eb*(L4b);
        const double Ht3 = as*(Ct2)+ec*(L6c-L3);


        // set elements, mapping to amplitude FF basis
        // Fs
        result.element({0}) = -((Hps * Mc) / sqMbMc);
        // Ff
        result.element({1}) = (Ha1 * ((Mb + Mc)*(Mb + Mc) - Sqq)) / (2. * sqMbMc);
        // Fg
        result.element({2}) = Hv / (2. * sqMbMc);
        // Fm
        result.element({3}) = -(Ha2 * sqrt(Mc / Mb3)) / 2. + Ha3 / (2. * sqMbMc);
        // Fp
        result.element({4}) = -(Ha2 * sqrt(Mc / Mb3)) / 2. - Ha3 / (2. * sqMbMc);
        // Fzt
        result.element({5}) = (Ht3 * Mc) / (2. * pow(Mb * Mc, 1.5));
        // Fmt
        result.element({6}) = (Ht1 * (Mb - Mc)) / (2. * sqMbMc) - (Ht2 * (Mb + Mc)) / (2. * sqMbMc);
        // Fpt
        result.element({7}) = (Ht2 * (Mb - Mc)) / (2. * sqMbMc) - (Ht1 * (Mb + Mc)) / (2. * sqMbMc);

        result *= h;

    }

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

} // namespace Hammer