code_browser/FFBtoDstarBLPRVar_8cc_source.html

 ///

 /// @file  FFBtoDstarBLPRVar.cc

 /// @brief \f$ B \rightarrow D^* \f$ BLPRVar form factors with variations

 ///


 //**** 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/FFBtoDstarBLPRVar.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;


     FFBtoDstarBLPRVar::FFBtoDstarBLPRVar() {

         // Create tensor rank and dimensions

         vector<IndexType> dims = {8, 8}; // size is _FFErrNames + 1 to include central values in zeroth component

         setGroup("BLPRVar"); //override BLPR base class FF group setting

         string name{"FFBtoDstarBLPRVar"};


         setPrefix("BtoD*");

         _FFErrLabel = FF_BDSTAR_VAR;

         addProcessSignature(PID::BPLUS, {-PID::DSTAR});

         addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BDSTAR, FF_BDSTAR_VAR})});


         addProcessSignature(PID::BZERO, {PID::DSTARMINUS});

         addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BDSTAR, FF_BDSTAR_VAR})});


         setPrefix("BstoDs*");

         _FFErrLabel = FF_BSDSSTAR_VAR;

         addProcessSignature(PID::BS, {PID::DSSTARMINUS});

         addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BSDSSTAR, FF_BSDSSTAR_VAR})});


         setSignatureIndex();

     }


     void FFBtoDstarBLPRVar::defineSettings() {

         _FFErrNames = {"delta_RhoSq","delta_chi21","delta_chi2p","delta_chi3p","delta_eta1","delta_etap","delta_dV20"};

         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>("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.);


         //correlation matrix and set zero error eigenvectors

         //Row basis is RhoSq, chi2(1), chi2(1)prime, chi3(1)prime, eta(1), eta(1)prime, V20

         vector<vector<double>> sliwmat{ {1., 0., 0., 0., 0., 0., 0.},

                                         {0., 1., 0., 0., 0., 0., 0.},

                                         {0., 0., 1., 0., 0., 0., 0.},

                                         {0., 0., 0., 1., 0., 0., 0.},

                                         {0., 0., 0., 0., 1., 0., 0.},

                                         {0., 0., 0., 0., 0., 1., 0.},

                                         {0., 0., 0., 0., 0., 0., 1.}};

         addSetting<vector<vector<double>>>("sliwmatrix",sliwmat);


         for (auto elem : _FFErrNames) {

             addSetting<double>(elem, 0.);

         }


         initialized = true;


     }


     void FFBtoDstarBLPRVar::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 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);


         //BLPRVar 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"));


         const vector<vector<double>>& sliwmat = (*getSetting<vector<vector<double>>>("sliwmatrix"));


         //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;

         const 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 ));

         const double xiIWN = 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 ));


         const double chi1=0;

         const double h=xiIW/xiIWN + 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, 0}) = -((Hps * Mc) / sqMbMc);

         // Ff

         result.element({1, 0}) = (Ha1 * ((Mb + Mc)*(Mb + Mc) - Sqq)) / (2. * sqMbMc);

         // Fg

         result.element({2, 0}) = Hv / (2. * sqMbMc);

         // Fm

         result.element({3, 0}) = -(Ha2 * sqrt(Mc / Mb3)) / 2. + Ha3 / (2. * sqMbMc);

         // Fp

         result.element({4, 0}) = -(Ha2 * sqrt(Mc / Mb3)) / 2. - Ha3 / (2. * sqMbMc);

         // Fzt

         result.element({5, 0}) = (Ht3 * Mc) / (2. * pow(Mb * Mc, 1.5));

         // Fmt

         result.element({6, 0}) = (Ht1 * (Mb - Mc)) / (2. * sqMbMc) - (Ht2 * (Mb + Mc)) / (2. * sqMbMc);

         // Fpt

         result.element({7, 0}) = (Ht2 * (Mb - Mc)) / (2. * sqMbMc) - (Ht1 * (Mb + Mc)) / (2. * sqMbMc);


         //Now create remaining FF tensor entries

         //n indexes rows ie FFs; idx indexes the columns that contract with delta.

         const double rC = Mc/Mb;

         const double sqrC = sqrt(rC);

         const double wSq = w*w;

         const double zConSq = zCon*zCon;

         const double wm1Sq = (w-1)*(w-1);

         const double am1Sq = (a-1)*(a-1);

         const double ap1Sq = (a+1)*(a+1);


         //Linearization of hatted FFs

         vector<vector<double>> FFmat{

             {0, 4*(eb - ec)*(-1 + w)*sqrC, 4*(eb - ec)*wm1Sq*sqrC, -4*(3*eb - ec)*(-1 + w)*sqrC, 2*(eb - ec)*sqrC, 2*(eb - ec)*(-1 + w)*sqrC, 0},

             {0, -4*eb*Mb*(-1 + wSq)*sqrC, -4*eb*Mb*wm1Sq*(1 + w)*sqrC, 4*(3*eb - ec)*Mb*(-1 + wSq)*sqrC, -2*eb*Mb*(-1 + w)*sqrC, -2*eb*Mb*wm1Sq*sqrC, 0},

             {0, (-2*eb*(-1 + w))/(Mb*sqrC), (-2*eb*wm1Sq)/(Mb*sqrC), (2*(3*eb - ec)*(-1 + w))/(Mb*sqrC), -(eb/(Mb*sqrC)), -((eb*(-1 + w))/(Mb*sqrC)), 0},

             {0, (-2*(eb*(-1 + w) + ec*(1 + rC)))/(Mb*sqrC), (-2*(-1 + w)*(eb*(-1 + w) + ec*(1 + rC)))/(Mb*sqrC), (2*(3*eb - ec)*(-1 + w))/(Mb*sqrC), -((eb + ec + eb*w - ec*rC)/(Mb*(1 + w)*sqrC)), -(((-1 + w)*(eb + ec + eb*w - ec*rC))/(Mb*(1 + w)*sqrC)), 0},

             {0, (2*(ec + eb*(-1 + w) - ec*rC))/(Mb*sqrC), (2*(-1 + w)*(ec + eb*(-1 + w) - ec*rC))/(Mb*sqrC), (-2*(3*eb - ec)*(-1 + w))/(Mb*sqrC), (eb + ec + eb*w + ec*rC)/(Mb*(1 + w)*sqrC), ((-1 + w)*(eb + ec + eb*w + ec*rC))/(Mb*(1 + w)*sqrC), 0},

             {0, (-2*ec)/(Mb2*sqrC), (-2*ec*(-1 + w))/(Mb2*sqrC), 0, -(ec/(Mb2*(1 + w)*sqrC)), -((ec*(-1 + w))/(Mb2*(1 + w)*sqrC)), 0},

             {0, (2*eb*(-1 + w)*(-1 + rC))/sqrC, (2*eb*wm1Sq*(-1 + rC))/sqrC, (-2*(3*eb - ec)*(-1 + w)*(-1 + rC))/sqrC, (eb*(1 + rC))/sqrC, (eb*(-1 + w)*(1 + rC))/sqrC, 0},

             {0, (2*eb*(-1 + w)*(1 + rC))/sqrC, (2*eb*wm1Sq*(1 + rC))/sqrC, (-2*(3*eb - ec)*(-1 + w)*(1 + rC))/sqrC, (eb*(-1 + rC))/sqrC, (eb*(-1 + w)*(-1 + rC))/sqrC, 0}};


         //Linearization of LO IW functions

         const vector<double> xiIWL = {

                 -8*a2*zCon + zConSq*(V21 + (64*a4*Cv1zp + (64*a4*Cv3z)/(1 + rD) + (128*a6*Cv3zp)/(1 + rD) + (64*a4*Cv2z*rD)/(1 + rD) + (128*a6*Cv2zp*rD)/(1 + rD))*as - 256*a4*(eb + ec)*chi21 - 16*(-64*a4 + 64*a6)*(eb + ec)*chi2p + 768*a4*(eb + ec)*chi3p - (128*a4*(-1 + rD)*(eb - ec)*etap)/(1 + rD)),

                 -4*zConSq*(eb + ec)*Xi,

                 -16*zConSq*(eb + ec)*(-(a2*(-16 + V21)) + V21 + 64*a6*RhoSq - 32*a4*(1 + 2*RhoSq)),

                 12*zConSq*(eb + ec)*Xi,

                 0,

                 (-2*(-1 + rD)*zConSq*(eb - ec)*Xi)/(1 + rD),

                 -zConSq};


         const vector<double> xiIWLN = {

                 (8*a4 + a2*(-8 + V21) + V21 - 2*a*V21)/ap1Sq + ((64*am1Sq*a4*Cv1zp)/ap1Sq + (64*am1Sq*a4*Cv3z)/(ap1Sq*(1 + rD)) + (128*am1Sq*a6*Cv3zp)/(ap1Sq*(1 + rD)) + (64*am1Sq*a4*Cv2z*rD)/(ap1Sq*(1 + rD)) + (128*am1Sq*a6*Cv2zp*rD)/(ap1Sq*(1 + rD)))*as - (256*am1Sq*a4*(eb + ec)*chi21)/ap1Sq - (16*am1Sq*(-64*a4 + 64*a6)*(eb + ec)*chi2p)/ap1Sq + (768*am1Sq*a4*(eb + ec)*chi3p)/ap1Sq - (128*am1Sq*a4*(-1 + rD)*(eb - ec)*etap)/(ap1Sq*(1 + rD)),

                 (-4*am1Sq*(eb + ec)*Xi)/ap1Sq,

                 (-16*am1Sq*(eb + ec)*(-(a2*(-16 + V21)) + V21 + 64*a6*RhoSq - 32*a4*(1 + 2*RhoSq)))/ap1Sq,

                 (12*am1Sq*(eb + ec)*Xi)/ap1Sq,

                 0,

                 (-2*am1Sq*(-1 + rD)*(eb - ec)*Xi)/(ap1Sq*(1 + rD)),

                 (-1 + 2*a - a2)/ap1Sq};


         for(size_t n = 0; n < 8; ++n){

             auto ni = static_cast<IndexType>(n);

             for(size_t idx = 0; idx < 7; ++idx){

                 complex<double> entry = 0;

                 auto idxi = static_cast<IndexType>(idx+1u);

                 for(size_t idx1 = 0; idx1 < 7; ++idx1){

                     entry += sliwmat[idx1][idx]*(FFmat[n][idx1] + result.element({ni, 0})*(xiIWL[idx1]/xiIW -  xiIWLN[idx1]/xiIWN));

                 }

                 result.element({ni, idxi}) = entry;

             }

         }


         result *= h;

         result.toVector();

     }


     std::unique_ptr<FormFactorBase> FFBtoDstarBLPRVar::clone(const std::string& label) {

         MAKE_CLONE(FFBtoDstarBLPRVar, label);

     }


 } // namespace Hammer

MSG_WARNING
#define MSG_WARNING(x)
Definition: Logging.hh:366

Hammer::MultiDimensional::makeEmptySparse
TensorData makeEmptySparse(const IndexList &dimensions, const LabelsList &labels)
Definition: SparseContainer.cc:287

Hammer::Tensor::element
std::complex< double > & element(const IndexList &indices={})
access an element given its indices
Definition: Tensor.cc:67

Hammer::Tensor::toVector
Tensor & toVector()
forces conversion of a tensor to vector type
Definition: Tensor.cc:171

FFBtoDstarBLPRVar.hh
 BLPRVar form factors with variations

IndexLabels.hh
Tensor indices label definitions.

Utils.hh

Hammer::FF_BSDSSTAR
Definition: IndexLabels.hh:158

Hammer::IndexType
uint16_t IndexType
Definition: MultiDimensional.fhh:23

Hammer::sqrt2
static constexpr double sqrt2
Definition: Constants.hh:26

Hammer::FF_BDSTAR
Definition: IndexLabels.hh:140

SparseContainer.hh
Sparse tensor data container.

Hammer::FF_BSDSSTAR_VAR
Definition: IndexLabels.hh:207

Hammer::Tensor
Multidimensional tensor class with complex numbers as elements.
Definition: Tensor.hh:33

Hammer::isZero
bool isZero(const std::complex< double > val)
Definition: Math/Utils.hh:25

Hammer::Tensor::clearData
void clearData()
sets all the elements to 0
Definition: Tensor.cc:229

Constants.hh
Various numerical constants.

Pdg.hh
Hammer particle data class.

Hammer::FFBtoDstarBLPRVar
Definition: FFBtoDstarBLPRVar.hh:19

Particle.hh
Hammer particle class.

MAKE_CLONE
#define MAKE_CLONE(OBJ, LABEL)
Definition: FormFactorBase.hh:149

Hammer::FF_BDSTAR_VAR
Definition: IndexLabels.hh:189

Hammer::pi
static constexpr double pi
Definition: Constants.hh:21