Hammer  1.0.0
Helicity Amplitude Module for Matrix Element Reweighting
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FFBtoD1BLRVar.cc
Go to the documentation of this file.
1 ///
2 /// @file FFBtoD1BLRVar.cc
3 /// @brief \f$ B \rightarrow D_1 \f$ BLR 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/FFBtoD1BLRVar.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 #include <iostream>
21 
22 using namespace std;
23 
24 namespace Hammer {
25 
26  namespace MD = MultiDimensional;
27 
28  FFBtoD1BLRVar::FFBtoD1BLRVar() {
29  // Create tensor rank and dimensions
30  vector<IndexType> dims = {8, 8}; // size is _FFErrNames + 1 to include central values in zeroth component
31  setGroup("BLRVar"); // override BLPR base class FF group setting
32  string name{"FFBtoD1BLRVar"};
33 
34  setPrefix("BtoD**1");
35  _FFErrLabel = FF_BDSSD1_VAR;
36  addProcessSignature(PID::BPLUS, {-PID::DSSD1});
37  addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BDSSD1, FF_BDSSD1_VAR})});
38 
39  addProcessSignature(PID::BZERO, {PID::DSSD1MINUS});
40  addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BDSSD1, FF_BDSSD1_VAR})});
41 
42  setPrefix("BstoDs**1");
43  _FFErrLabel = FF_BSDSSDS1_VAR;
44  addProcessSignature(PID::BS, {PID::DSSDS1MINUS});
45  addTensor(Tensor{name, MD::makeEmptySparse(dims, {FF_BSDSSDS1, FF_BSDSSDS1_VAR})});
46 
47  setSignatureIndex();
48  }
49 
50  void FFBtoD1BLRVar::defineSettings() {
51  _FFErrNames = {"delta_tau1", "delta_tau2", "delta_t1", "delta_tp", "delta_eta1", "delta_eta2", "delta_eta3"};
52  //_FFErrNames = ;
53  setPath(getFFErrPrefixGroup().get());
54  setUnits("GeV");
55 
56  //1S scheme: mb1S = 4710, mbBar = 5313, lambda1 = -3 10^5 MeV^2, delta mb-mc = 3400, alpha_s = 26/100
57  addSetting<double>("as",0.26);
58  addSetting<double>("mb",4.710);
59  addSetting<double>("mc",4.710 - 3.400);
60  addSetting<double>("t1", 0.7);
61  addSetting<double>("tp", -1.6);
62  addSetting<double>("tau1", -0.5);
63  addSetting<double>("tau2", 2.9);
64  addSetting<double>("eta1", 0.);
65  addSetting<double>("eta2", 0.);
66  addSetting<double>("eta3", 0.);
67  addSetting<double>("laB", 0.4);
68  addSetting<double>("laP", 0.8);
69 
70  // correlation matrix and set zero error eigenvectors
71  // Row basis is t1, tp, tau1, tau2, eta1(1), eta2(1), eta3(1)
72  vector<vector<double>> sliwmat{ {1., 0., 0., 0., 0., 0., 0.},
73  {0., 1., 0., 0., 0., 0., 0.},
74  {0., 0., 1., 0., 0., 0., 0.},
75  {0., 0., 0., 1., 0., 0., 0.},
76  {0., 0., 0., 0., 1., 0., 0.},
77  {0., 0., 0., 0., 0., 1., 0.},
78  {0., 0., 0., 0., 0., 0., 1.}};
79  addSetting<vector<vector<double>>>("sliwmatrix",sliwmat);
80 
81  for (auto elem : _FFErrNames) {
82  addSetting<double>(elem, 0.);
83  }
84 
85  initialized = true;
86  }
87 
88  void FFBtoD1BLRVar::evalAtPSPoint(const vector<double>& point, const vector<double>& masses) {
89  Tensor& result = getTensor();
90  result.clearData();
91 
92  if(!initialized){
93  MSG_WARNING("Warning, Setting have not been defined!");
94  }
95 
96  double Mb = 0.;
97  double Mc = 0.;
98  // double unitres = 1.;
99  if(masses.size() >= 2) {
100  Mb = masses[0];
101  Mc = masses[1];
102  //unitres = _units;
103  }
104  else {
105  Mb = this->masses()[0];
106  Mc = this->masses()[1];
107  }
108  const double Mb2 = Mb*Mb;
109  const double Mc2 = Mc*Mc;
110 
111  double Sqq = point[0];
112 
113  // const double sqSqq = sqrt(Sqq);
114  double w = (Mb2 + Mc2 - Sqq)/(2.*Mb*Mc);
115  //safety measure if w==1.0
116  if(isZero(w - 1.0)) w += 1e-6;
117 
118  //BLR parameters
119  const double zBC = (*getSetting<double>("mc"))/(*getSetting<double>("mb"));
120  const double eB = 1./(*getSetting<double>("mb")*2.);
121  const double eC = 1./(*getSetting<double>("mc")*2.);
122  const double as = (*getSetting<double>("as"))/pi;
123  const double t1 = (*getSetting<double>("t1"));
124  const double tp = (*getSetting<double>("tp"));
125  const double tau1 = (*getSetting<double>("tau1"));
126  const double tau2 = (*getSetting<double>("tau2"));
127  const double eta1 = (*getSetting<double>("eta1"));
128  const double eta2 = (*getSetting<double>("eta2"));
129  const double eta3 = (*getSetting<double>("eta3"));
130  const double laB = (*getSetting<double>("laB"));
131  const double laP = (*getSetting<double>("laP"));
132 
133  const vector<vector<double>>& sliwmat = (*getSetting<vector<vector<double>>>("sliwmatrix"));
134 
135  const double LambdaD32 = -laB + laP * w;
136  const double Fb = laB + laP - tau2 - tau1*(1 + 2*w);
137  const double LOIWtau = t1 + (w-1)*t1*tp;
138 
139  //QCD correction functions
140  const double Cs = CS(w, zBC);
141  // const double Cps = CP(w, zBC);
142  const double Cv1 = CV1(w, zBC);
143  const double Cv2 = CV2(w, zBC);
144  const double Cv3 = CV3(w, zBC);
145  const double Ca1 = CA1(w, zBC);
146  // const double Ca2 = CA2(w, zBC);
147  // const double Ca3 = CA3(w, zBC);
148  const double Ct1 = CT1(w, zBC);
149  const double Ct2 = CT2(w, zBC);
150  const double Ct3 = CT3(w, zBC);
151 
152  //Form factors
153  const double Fsc = (-(eC*(4*LambdaD32 + 2*(1 + w)*(6*eta1 + 2*(w-1)*eta2 - eta3) - 2*(w-1)*((1 + 2*w)*tau1 + tau2))) - 2*(1 + w)*(1 + as*Cs) - 2*(w-1)*eB*Fb)/sqrt(6.);
154  const double Fv1 = (-(eC*(4*(1 + w)*LambdaD32 - (-1 + w*w)*(2*eta1 + 3*eta3 + 3*tau1 - 3*tau2))) + (1 - w*w)*(1 + as*Cv1) - (-1 + w*w)*eB*Fb)/sqrt(6.);
155  const double Fv2 = (-3 - eC*(10*eta1 + 4*(w-1)*eta2 - 5*eta3 + (-1 + 4*w)*tau1 + 5*tau2) - as*(3*Cv1 + 2*(1 + w)*Cv2) - 3*eB*Fb)/sqrt(6.);
156  const double Fv3 = (-2 + w + eC*(4*LambdaD32 - 2*(6 + w)*eta1 - 4*(w-1)*eta2 - (-2 + 3*w)*eta3 + (2 + w)*tau1 + (2 + 3*w)*tau2) - as*((2 - w)*Cv1 + 2*(1 + w)*Cv3) + (2 + w)*eB*Fb)/sqrt(6.);
157  const double Fa = (-(eC*(4*LambdaD32 - (1 + w)*(2*eta1 + 3*eta3) - 3*(w-1)*(tau1 - tau2))) + (-1 - w)*(1 + as*Ca1) - (w-1)*eB*Fb)/sqrt(6.);
158  const double Ft1 = (-(eC*(4*LambdaD32 + (1 + w)*(2*eta1 + 3*eta3) - 3*(w-1)*(tau1 - tau2))) + (1 + w)*(1 + as*(Ct1 + (w-1)*Ct2)) + (w-1)*eB*Fb)/sqrt(6.);
159  const double Ft2 = (-(eC*(4*LambdaD32 - (1 + w)*(2*eta1 + 3*eta3) - 3*(w-1)*(tau1 - tau2))) + (-1 - w)*(1 + as*(Ct1 - (w-1)*Ct3)) + (w-1)*eB*Fb)/sqrt(6.);
160  const double Ft3 = (3 - eC*(-10*eta1 - 4*(w-1)*eta2 + 5*eta3 + (-1 + 4*w)*tau1 + 5*tau2) + as*(3*Ct1 - (2 - w)*Ct2 + 3*Ct3) + 3*eB*Fb)/sqrt(6.);
161 
162 
163  // defined as 1/LOIWtau * D[LOIWtau * f, var]
164  const array<double, 7> FscDer = {
165  Fsc / t1, Fsc / LOIWtau * t1 * (w - 1), sqrt(2./3.)*(eB + eC)*(-1 + w)*(1 + 2*w), sqrt(2./3.)*(eB + eC)*(-1 + w),
166  -2*sqrt(6.)*eC*(1 + w), -2*sqrt(2./3.)* eC*(-1 + w*w), sqrt(2./3.)* eC*(1 + w)};
167  const array<double, 7> Fv1Der = {Fv1 / t1, Fv1/LOIWtau * t1 * (w-1), ((eB + 3*eC + 2*eB*w)*(-1 + w*w))/sqrt(6.),
168  ((eB - 3*eC)*(-1 + w*w))/sqrt(6.), sqrt(2./3.)*eC*(-1 + w*w), 0., sqrt(3./2.)*eC*(-1 + w*w)};
169  const array<double, 7> Fv2Der = {Fv2 / t1, Fv2 / LOIWtau * t1 * (w - 1), (eC - 4*eC*w + eB*(3 + 6*w)) / sqrt(6.),
170  (3*eB - 5*eC) / sqrt(6.), -5*sqrt(2./3.)* eC, -2*sqrt(2./3.)* eC*(-1 + w), (5*eC) / sqrt(6.)};
171  const array<double, 7> Fv3Der = {Fv3 / t1, Fv3 / LOIWtau * t1 * (w - 1), -(((2 + w)*(eB - eC + 2*eB*w)) / sqrt(6.)),
172  (-eB*(2 + w) + eC*(2 + 3*w)) / sqrt(6.), -sqrt(2. / 3.)*eC*(6 + w), -2*sqrt(2./3.)* eC*(-1 + w), (eC*(2 - 3*w)) / sqrt(6.)};
173  const array<double, 7> FaDer = {Fa / t1, Fa / LOIWtau * t1 * (w - 1), ((-1 + w)*(eB + 3*eC + 2*eB*w)) / sqrt(6.),
174  ((eB - 3*eC)*(-1 + w)) / sqrt(6.), sqrt(2./3.)* eC*(1 + w), 0., sqrt(3. / 2.)*eC*(1 + w)};
175  const array<double, 7> Ft1Der = {Ft1 / t1, Ft1 / LOIWtau * t1 * (w - 1), -(((-1 + w)*(eB - 3*eC + 2*eB*w)) / sqrt(6.)),
176  -(((eB + 3*eC)*(-1 + w)) / sqrt(6.)), -sqrt(2. / 3.)*eC*(1 + w), 0., -sqrt(3. / 2.)*eC*(1 + w)};
177  const array<double, 7> Ft2Der = {Ft2 / t1, Ft2 / LOIWtau * t1 * (w - 1), -(((-1 + w)*(eB - 3*eC + 2*eB*w)) / sqrt(6.)),
178  -(((eB + 3*eC)*(-1 + w)) / sqrt(6.)), sqrt(2./3.)* eC*(1 + w), 0., sqrt(3. / 2.)*eC*(1 + w)};
179  const array<double, 7> Ft3Der = {Ft3 / t1, Ft3 / LOIWtau * t1 * (w - 1), (eC - 4*eC*w - 3*eB*(1 + 2*w)) / sqrt(6.),
180  -((3*eB + 5*eC) / sqrt(6.)), 5*sqrt(2./3.)* eC, 2*sqrt(2./3.)* eC*(-1 + w), -((5*eC) / sqrt(6.))};
181 
182  //Set elements
183  result.element({0, 0}) = Fsc;
184  result.element({1, 0}) = Fv1;
185  result.element({2, 0}) = Fv2;
186  result.element({3, 0}) = Fv3;
187  result.element({4, 0}) = Fa;
188  result.element({5, 0}) = Ft1;
189  result.element({6, 0}) = Ft2;
190  result.element({7, 0}) = Ft3;
191 
192  for (IndexType i1 = 1; i1 <= 7; ++i1) {
193  for (size_t i2 = 0; i2 < 7; ++i2) {
194  result.element({0, i1}) += sliwmat[i2][i1 - 1] * FscDer[i2];
195  result.element({1, i1}) += sliwmat[i2][i1 - 1] * Fv1Der[i2];
196  result.element({2, i1}) += sliwmat[i2][i1 - 1] * Fv2Der[i2];
197  result.element({3, i1}) += sliwmat[i2][i1 - 1] * Fv3Der[i2];
198  result.element({4, i1}) += sliwmat[i2][i1 - 1] * FaDer[i2];
199  result.element({5, i1}) += sliwmat[i2][i1 - 1] * Ft1Der[i2];
200  result.element({6, i1}) += sliwmat[i2][i1 - 1] * Ft2Der[i2];
201  result.element({7, i1}) += sliwmat[i2][i1 - 1] * Ft3Der[i2];
202  }
203  }
204 
205  result *= LOIWtau;
206  result.toVector();
207  }
208 
209  std::unique_ptr<FormFactorBase> FFBtoD1BLRVar::clone(const std::string& label) {
210  MAKE_CLONE(FFBtoD1BLRVar, label);
211  }
212 
213 } // 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
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.
Pdg.hh
Hammer particle data class.
FFBtoD1BLRVar.hh
BLR form factors with variations
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