Hammer  1.0.0
Helicity Amplitude Module for Matrix Element Reweighting
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FFBtoDBLPRVar.cc
Go to the documentation of this file.
1 ///
2 /// @file FFBtoDBLPRVar.cc
3 /// @brief \f$ B \rightarrow D \f$ BLPRVar form factors with variations
4 ///
5 
6 //**** This file is a part of the HAMMER library
7 //**** Copyright (C) 2016 - 2020 The HAMMER Collaboration
8 //**** HAMMER is licensed under version 3 of the GPL; see COPYING for details
9 //**** Please note the MCnet academic guidelines; see GUIDELINES for details
10 
11 // -*- C++ -*-
12 #include "Hammer/FormFactors/FFBtoDBLPRVar.hh"
13 #include "Hammer/IndexLabels.hh"
14 #include "Hammer/Math/Constants.hh"
15 #include "Hammer/Math/Utils.hh"
16 #include "Hammer/Particle.hh"
17 #include "Hammer/Tools/Pdg.hh"
18 #include "Hammer/Math/MultiDim/SparseContainer.hh"
19 #include <cmath>
20 
21 #include <iostream>
22 
23 using namespace std;
24 
25 namespace Hammer {
26 
27  namespace MD = MultiDimensional;
28 
29  FFBtoDBLPRVar::FFBtoDBLPRVar() {
30  // Create tensor rank and dimensions
31  vector<IndexType> dims = {4, 8}; // size is _FFErrNames + 1 to include central values in zeroth component
32  setGroup("BLPRVar"); //override BLPR base class FF group setting
33  string name{"FFBtoDBLPRVar"};
34 
35  setPrefix("BtoD");
36  _FFErrLabel = FF_BD_VAR;
37  addProcessSignature(PID::BPLUS, {-PID::D0});
38  addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BD, FF_BD_VAR})});
39 
40  addProcessSignature(PID::BZERO, {PID::DMINUS});
41  addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BD, FF_BD_VAR})});
42 
43  setPrefix("BstoDs");
44  _FFErrLabel = FF_BSDS_VAR;
45  addProcessSignature(PID::BS, {PID::DSMINUS});
46  addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BSDS, FF_BSDS_VAR})});
47 
48  setSignatureIndex();
49  }
50 
51  void FFBtoDBLPRVar::defineSettings() {
52  _FFErrNames = {"delta_RhoSq","delta_chi21","delta_chi2p","delta_chi3p","delta_eta1","delta_etap","delta_dV20"};
53  setPath(getFFErrPrefixGroup().get());
54  setUnits("GeV");
55 
56  addSetting<double>("RhoSq",1.24);
57  addSetting<double>("a",1.509/sqrt2);
58  //1S scheme: mb1S = 4710, mbBar = 5313, lambda1 = -3 10^5 MeV^2, delta mb-mc = 3400, alpha_s = 26/100
59  addSetting<double>("as",0.26);
60  addSetting<double>("mb",4.710);
61  addSetting<double>("mc",4.710 - 3.400);
62  addSetting<double>("la",0.57115);
63  //Renormalon improvement (1S scheme)
64  addSetting<double>("ebR/eb",0.861);
65  addSetting<double>("ecR/ec",0.822);
66 
67  addSetting<double>("chi21",-0.06);
68  addSetting<double>("chi2p",-0.00);
69  addSetting<double>("chi3p",0.05);
70  addSetting<double>("eta1",0.30);
71  addSetting<double>("etap",-0.05);
72  addSetting<double>("dV20",0.);
73 
74  //correlation matrix and set zero error eigenvectors
75  //Row basis is RhoSq, chi2(1), chi2(1)prime, chi3(1)prime, eta(1), eta(1)prime, V20
76  vector<vector<double>> sliwmat{ {1., 0., 0., 0., 0., 0., 0.},
77  {0., 1., 0., 0., 0., 0., 0.},
78  {0., 0., 1., 0., 0., 0., 0.},
79  {0., 0., 0., 1., 0., 0., 0.},
80  {0., 0., 0., 0., 1., 0., 0.},
81  {0., 0., 0., 0., 0., 1., 0.},
82  {0., 0., 0., 0., 0., 0., 1.}};
83  addSetting<vector<vector<double>>>("sliwmatrix",sliwmat);
84 
85  for (auto elem : _FFErrNames) {
86  addSetting<double>(elem, 0.);
87  }
88 
89  initialized = true;
90 
91  }
92 
93  void FFBtoDBLPRVar::evalAtPSPoint(const vector<double>& point, const vector<double>& masses) {
94  Tensor& result = getTensor();
95  result.clearData();
96 
97  if(!initialized){
98  MSG_WARNING("Warning, Settings have not been defined!");
99  }
100 
101  double Mb = 0.;
102  double Mc = 0.;
103  // double unitres = 1.;
104  if(masses.size() >= 2) {
105  Mb = masses[0];
106  Mc = masses[1];
107  // unitres = _units;
108  }
109  else {
110  Mb = this->masses()[0];
111  Mc = this->masses()[1];
112  }
113  const double Mb2 = Mb*Mb;
114  const double Mc2 = Mc*Mc;
115 
116  double Sqq = point[0];
117 
118  // const double sqSqq = sqrt(Sqq);
119  const double sqMbMc = sqrt(Mb*Mc);
120  double w = (Mb2 + Mc2 - Sqq)/(2.*Mb*Mc);
121  //safety measure if w==1.0
122  if(isZero(w - 1.0)) w += 1e-6;
123 
124  //Conformal parameters
125  const double a = *getSetting<double>("a");
126  const double zCon = (sqrt(w+1) - sqrt2*a)/(sqrt(w+1) + sqrt2*a);
127 
128  //BLPRVar expansion parameters
129  const double zBC = (*getSetting<double>("mc"))/(*getSetting<double>("mb"));
130  const double eb = (*getSetting<double>("la"))/(*getSetting<double>("mb")*2.);
131  const double ebReb = (*getSetting<double>("ebR/eb"));
132  const double ec = (*getSetting<double>("la"))/(*getSetting<double>("mc")*2.);
133  const double ecRec = (*getSetting<double>("ecR/ec"));
134  const double as = (*getSetting<double>("as"))/pi;
135 
136  //(Sub)leading IW function parameters
137  const double RhoSq = *getSetting<double>("RhoSq");
138  const double chi21 = (*getSetting<double>("chi21"));
139  const double chi2p = (*getSetting<double>("chi2p"));
140  const double chi3p = (*getSetting<double>("chi3p"));
141  const double eta1 = (*getSetting<double>("eta1"));
142  const double etap = (*getSetting<double>("etap"));
143 
144  const vector<vector<double>>& sliwmat = (*getSetting<vector<vector<double>>>("sliwmatrix"));
145 
146  //L and helper variables
147  const double L1 = -4.*(w-1)*(chi21 + (w-1.)*chi2p)+12.*chi3p*(w-1.);
148  // const double L2 = -4.*chi3p*(w-1.);
149  // const double L3 = 4.*(chi21 + (w-1.)*chi2p);
150  // const double L4 = 2.*(eta1 + etap*(w-1.))-1.;
151  const double L4b = (2.*(eta1 + etap*(w-1.))-1.*ebReb);
152  const double L4c = (2.*(eta1 + etap*(w-1.))-1.*ecRec);
153  // const double L5 = -1.;
154  // const double L5c = -ecRec;
155  // const double L6 = -2.*(1+(eta1 + etap*(w-1.)))/(w+1.);
156  // const double L6c = -2.*(ecRec + (eta1 + etap*(w-1.)))/(w+1.);
157 
158  // Variables for leading IW function (derived from G(1))
159  const double rD = 1867./5280.;
160  // const double V11 = 8*pow(a,2);
161  const double V21 = 57.0;
162  const double V20 = 7.5 + (*getSetting<double>("dV20"));
163 
164  //QCD correction functions
165  const double Cs = CS(w, zBC);
166  // const double Cps = CP(w, zBC);
167  const double Cv1 = CV1(w, zBC);
168  const double Cv2 = CV2(w, zBC);
169  const double Cv3 = CV3(w, zBC);
170  // const double Ca1 = CA1(w, zBC);
171  // const double Ca2 = CA2(w, zBC);
172  // const double Ca3 = CA3(w, zBC);
173  const double Ct1 = CT1(w, zBC);
174  const double Ct2 = CT2(w, zBC);
175  const double Ct3 = CT3(w, zBC);
176  //Derivatives (approx) at w_0 = 2a^2-1. Calculated once and stored.
177  if(!initCs){
178  const double w0 = 2*a*a - 1.;
179  Cv1z = CV1(w0, zBC);
180  Cv2z = CV2(w0, zBC);
181  Cv3z = CV3(w0, zBC);
182 
183  Cv1zp = (CV1(w0 + 1e-5, zBC) - Cv1z)/1e-5;
184  Cv2zp = (CV2(w0 + 1e-5, zBC) - Cv2z)/1e-5;
185  Cv3zp = (CV3(w0 + 1e-5, zBC) - Cv3z)/1e-5;
186 
187  Cv1zpp = (Cv1zp - (Cv1z - CV1(w0 - 1e-5, zBC))/1e-5)/1e-5;
188  Cv2zpp = (Cv2zp - (Cv2z - CV2(w0 - 1e-5, zBC))/1e-5)/1e-5;
189  Cv3zpp = (Cv3zp - (Cv3z - CV3(w0 - 1e-5, zBC))/1e-5)/1e-5;
190 
191  initCs = true;
192  }
193 
194 // //alpha_s derivatives at w_0 = 2a^2-1
195 // // const double Cv1z = 0.97543966986296275;
196 // const double Cv2z = -0.416360981425492067;
197 // const double Cv3z = -0.186896365253791276;
198 //
199 // const double Cv1zp = 0.11470494247949626;
200 // const double Cv2zp = 0.16684760479981954;
201 // const double Cv3zp = 0.038382492791834666;
202 //
203 // const double Cv1zpp = -0.2557645997620378;
204 // const double Cv2zpp = -0.1297628996251381;
205 // const double Cv3zpp = -0.02238541998227300;
206 
207 
208  // Leading and sub-leading IW function
209  //-------------------------------------------------------------------------------------------
210  const double a2 = a*a;
211  const double a4 = a2*a2;
212  const double a6 = a4*a2;
213  const double Xi = 64*a4*RhoSq - 16*a2 - V21;
214  const double xiIW = 1 - 8*a2*RhoSq*zCon + zCon*zCon*(
215  V21*RhoSq - V20 + (eb - ec)*(2*Xi*etap * (1-rD)/(1+rD))
216  + (eb + ec)*(Xi* (12*chi3p - 4*chi21) - 16*((a2-1)*Xi - 16* a4)*chi2p)
217  + 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 ));
218  const double xiIWN = 1 - 8*a2*RhoSq*(1-a)/(1+a) + pow((1-a)/(1+a),2.)*(
219  V21*RhoSq - V20 + (eb - ec)*(2*Xi*etap * (1-rD)/(1+rD))
220  + (eb + ec)*(Xi* (12*chi3p - 4*chi21) - 16*((a2-1)*Xi - 16* a4)*chi2p)
221  + 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 ));
222 
223  const double chi1=0;
224  const double h=xiIW/xiIWN + 2.*(eb+ec)*chi1;
225 
226  // Create hatted h FFs
227  const double Hs = 1.+as*Cs + ec*(L1-(L4c*(w-1)/(w+1))) + eb*(L1-(L4b*(w-1)/(w+1)));
228  const double Hp = 1.+as*(Cv1+0.5*(w+1)*(Cv2+Cv3))+(ec+eb)*L1;
229  const double Hm = as*0.5*(w+1)*(Cv2-Cv3)+(ec*L4c-eb*L4b);
230  const double Ht = 1.+as*(Ct1-Ct2+Ct3)+(ec+eb)*(L1) - (ec*L4c+eb*L4b);
231 
232 
233  // set elements, mapping to amplitude FF basis
234  // Fs
235  result.element({0, 0}) = (Hs * ((Mb + Mc)*(Mb + Mc) - Sqq)) / (2. * sqMbMc);
236  // Fz
237  result.element({1, 0}) = (Hp * ((Mb + Mc)*(Mb + Mc) - Sqq)) / (2. * sqMbMc * (Mb + Mc)) +
238  (Hm * (-(Mb - Mc)*(Mb - Mc) + Sqq)) / (2. * (Mb - Mc) * sqMbMc);
239  // Fp
240  result.element({2, 0}) = (Hm * (-Mb + Mc)) / (2. * sqMbMc) + (Hp * (Mb + Mc)) / (2. * sqMbMc);
241  // Ft
242  result.element({3, 0}) = Ht / (2. * sqMbMc);
243 
244  //Now create remaining FF tensor entries
245  //n indexes rows ie FFs; idx indexes the columns that contract with delta.
246  const double rC = Mc/Mb;
247  const double sqrC = sqrt(rC);
248  const double wSq = w*w;
249  const double zConSq = zCon*zCon;
250  const double wm1Sq = (w-1)*(w-1);
251  const double am1Sq = (a-1)*(a-1);
252  const double ap1Sq = (a+1)*(a+1);
253 
254  //Linearization of hatted FFs
255  vector<vector<double>> FFmat{
256  {0, -4*(eb + ec)*Mb*(-1 + wSq)*sqrC, -4*(eb + ec)*Mb*wm1Sq*(1 + w)*sqrC, 12*(eb + ec)*Mb*(-1 + wSq)*sqrC, -2*(eb + ec)*Mb*(-1 + w)*sqrC, -2*(eb + ec)*Mb*wm1Sq*sqrC, 0},
257  {0, (-4*(eb + ec)*(-1 + wSq)*sqrC)/(1 + rC), (-4*(eb + ec)*wm1Sq*(1 + w)*sqrC)/(1 + rC), (12*(eb + ec)*(-1 + wSq)*sqrC)/(1 + rC), (-2*(eb - ec)*(-1 + w)*sqrC)/(-1 + rC), (-2*(eb - ec)*wm1Sq*sqrC)/(-1 + rC), 0},
258  {0, (-2*(eb + ec)*(-1 + w)*(1 + rC))/sqrC, (-2*(eb + ec)*wm1Sq*(1 + rC))/sqrC, (6*(eb + ec)*(-1 + w)*(1 + rC))/sqrC, -(((eb - ec)*(-1 + rC))/sqrC), -(((eb - ec)*(-1 + w)*(-1 + rC))/sqrC), 0},
259  {0, (-2*(eb + ec)*(-1 + w))/(Mb*sqrC), (-2*(eb + ec)*wm1Sq)/(Mb*sqrC), (6*(eb + ec)*(-1 + w))/(Mb*sqrC), -((eb + ec)/(Mb*sqrC)), -(((eb + ec)*(-1 + w))/(Mb*sqrC)), 0}};
260 
261  //Linearization of LO IW functions
262  const vector<double> xiIWL = {
263  -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)),
264  -4*zConSq*(eb + ec)*Xi,
265  -16*zConSq*(eb + ec)*(-(a2*(-16 + V21)) + V21 + 64*a6*RhoSq - 32*a4*(1 + 2*RhoSq)),
266  12*zConSq*(eb + ec)*Xi,
267  0,
268  (-2*(-1 + rD)*zConSq*(eb - ec)*Xi)/(1 + rD),
269  -zConSq};
270 
271  const vector<double> xiIWLN = {
272  (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)),
273  (-4*am1Sq*(eb + ec)*Xi)/ap1Sq,
274  (-16*am1Sq*(eb + ec)*(-(a2*(-16 + V21)) + V21 + 64*a6*RhoSq - 32*a4*(1 + 2*RhoSq)))/ap1Sq,
275  (12*am1Sq*(eb + ec)*Xi)/ap1Sq,
276  0,
277  (-2*am1Sq*(-1 + rD)*(eb - ec)*Xi)/(ap1Sq*(1 + rD)),
278  (-1 + 2*a - a2)/ap1Sq};
279 
280  for(size_t n = 0; n < 4; ++n){
281  auto ni = static_cast<IndexType>(n);
282  for(size_t idx = 0; idx < 7; ++idx){
283  complex<double> entry = 0;
284  auto idxi = static_cast<IndexType>(idx+1u);
285  for(size_t idx1 = 0; idx1 < 7; ++idx1){
286  entry += sliwmat[idx1][idx]*(FFmat[n][idx1] + result.element({ni, 0})*(xiIWL[idx1]/xiIW - xiIWLN[idx1]/xiIWN));
287  }
288  result.element({ni, idxi}) = entry;
289  }
290  }
291 
292  result *= h;
293  result.toVector();
294  }
295 
296  std::unique_ptr<FormFactorBase> FFBtoDBLPRVar::clone(const std::string& label) {
297  MAKE_CLONE(FFBtoDBLPRVar, label);
298  }
299 
300 } // 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
IndexLabels.hh
Tensor indices label definitions.
Hammer::IndexType
uint16_t IndexType
Definition: MultiDimensional.fhh:23
Hammer::sqrt2
static constexpr double sqrt2
Definition: Constants.hh:26
SparseContainer.hh
Sparse tensor data container.
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.
FFBtoDBLPRVar.hh
BLPRVar form factors with variations
Pdg.hh
Hammer particle data class.
Particle.hh
Hammer particle class.
MAKE_CLONE
#define MAKE_CLONE(OBJ, LABEL)
Definition: FormFactorBase.hh:149
Hammer::pi
static constexpr double pi
Definition: Constants.hh:21