Hammer  1.0.0
Helicity Amplitude Module for Matrix Element Reweighting
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RateBDLepNu.cc
Go to the documentation of this file.
1 ///
2 /// @file RateBDLepNu.cc
3 /// @brief \f$ B \rightarrow D \tau\nu \f$ total rate
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/Rates/RateBDLepNu.hh"
13 #include "Hammer/IndexLabels.hh"
14 #include "Hammer/Math/Constants.hh"
15 #include "Hammer/Particle.hh"
16 #include "Hammer/Tools/Pdg.hh"
17 #include "Hammer/Math/MultiDim/SparseContainer.hh"
18 #include <cmath>
19 
20 using namespace std;
21 
22 namespace Hammer {
23 
24  namespace MD = MultiDimensional;
25 
26  RateBDLepNu::RateBDLepNu() {
27  // Create tensor rank and dimensions
28  vector<IndexType> dims{{11, 4, 11, 4, _nPoints}};
29  string name{"RateBDLepNuQ2"};
30  auto& pdg = PID::instance();
31  addProcessSignature(PID::BPLUS, {-PID::D0, PID::NU_TAU, PID::ANTITAU});
32  addIntegrationBoundaries({PS::makeQ2Function(pdg.getMass(PID::ANTITAU), pdg.getMass(PID::BPLUS)-pdg.getMass(PID::D0))});
33  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCTAUNU, FF_BD, WILSON_BCTAUNU_HC, FF_BD_HC, INTEGRATION_INDEX})});
34 
35  addProcessSignature(PID::BZERO, {PID::DMINUS, PID::NU_TAU, PID::ANTITAU});
36  addIntegrationBoundaries({PS::makeQ2Function(pdg.getMass(PID::ANTITAU), pdg.getMass(PID::BZERO)-pdg.getMass(PID::DMINUS))});
37  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCTAUNU, FF_BD, WILSON_BCTAUNU_HC, FF_BD_HC, INTEGRATION_INDEX})});
38 
39  addProcessSignature(PID::BPLUS, {-PID::D0, PID::NU_MU, PID::ANTIMUON});
40  addIntegrationBoundaries({PS::makeQ2Function(pdg.getMass(PID::ANTIMUON), pdg.getMass(PID::BPLUS)-pdg.getMass(PID::D0))});
41  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCMUNU, FF_BD, WILSON_BCMUNU_HC, FF_BD_HC, INTEGRATION_INDEX})});
42 
43  addProcessSignature(PID::BZERO, {PID::DMINUS, PID::NU_MU, PID::ANTIMUON});
44  addIntegrationBoundaries({PS::makeQ2Function(pdg.getMass(PID::ANTIMUON), pdg.getMass(PID::BZERO)-pdg.getMass(PID::DMINUS))});
45  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCMUNU, FF_BD, WILSON_BCMUNU_HC, FF_BD_HC, INTEGRATION_INDEX})});
46 
47  addProcessSignature(PID::BPLUS, {-PID::D0, PID::NU_E, PID::POSITRON});
48  addIntegrationBoundaries({PS::makeQ2Function(pdg.getMass(PID::POSITRON), pdg.getMass(PID::BPLUS)-pdg.getMass(PID::D0))});
49  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCENU, FF_BD, WILSON_BCENU_HC, FF_BD_HC, INTEGRATION_INDEX})});
50 
51  addProcessSignature(PID::BZERO, {PID::DMINUS, PID::NU_E, PID::POSITRON});
52  addIntegrationBoundaries({PS::makeQ2Function(pdg.getMass(PID::POSITRON), pdg.getMass(PID::BZERO)-pdg.getMass(PID::DMINUS))});
53  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCENU, FF_BD, WILSON_BCENU_HC, FF_BD_HC, INTEGRATION_INDEX})});
54 
55  //bs -> cs
56  addProcessSignature(PID::BS, {PID::DSMINUS, PID::NU_TAU, PID::ANTITAU});
57  addIntegrationBoundaries({PS::makeQ2Function(pdg.getMass(PID::ANTITAU), pdg.getMass(PID::BS)-pdg.getMass(PID::DSMINUS))});
58  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCTAUNU, FF_BSDS, WILSON_BCTAUNU_HC, FF_BSDS_HC, INTEGRATION_INDEX})});
59 
60  addProcessSignature(PID::BS, {PID::DSMINUS, PID::NU_MU, PID::ANTIMUON});
61  addIntegrationBoundaries({PS::makeQ2Function(pdg.getMass(PID::ANTIMUON), pdg.getMass(PID::BS)-pdg.getMass(PID::DSMINUS))});
62  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCMUNU, FF_BSDS, WILSON_BCMUNU_HC, FF_BSDS_HC, INTEGRATION_INDEX})});
63 
64  addProcessSignature(PID::BS, {PID::DSMINUS, PID::NU_E, PID::POSITRON});
65  addIntegrationBoundaries({PS::makeQ2Function(pdg.getMass(PID::POSITRON), pdg.getMass(PID::BS)-pdg.getMass(PID::DSMINUS))});
66  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCENU, FF_BSDS, WILSON_BCENU_HC, FF_BSDS_HC, INTEGRATION_INDEX})});
67 
68  setSignatureIndex();
69  }
70 
71  Tensor RateBDLepNu::evalAtPSPoint(const vector<double>& point) {
72  auto labs = getTensor().labels();
73  labs.pop_back();
74  auto dimensions = getTensor().dims();
75  dimensions.pop_back();
76  Tensor result{"RateBDLepNu", MD::makeEmptySparse(dimensions, labs)};
77 
78  const double Mb = masses()[0];
79  const double Mc = masses()[1];
80  const double Mt = masses()[3];
81 
82  const double Mb2 = Mb*Mb;
83  const double Mc2 = Mc*Mc;
84  const double Mt2 = Mt*Mt;
85 
86  // kinematic objects
87  const double Sqq = point[0];
88  const double Sqq2 = pow(Sqq, 2.);
89  const double Ew = (Mb2 - Mc2 + Sqq) / (2 * Mb);
90  const double Ew2 = Ew*Ew;
91  const double Pw = sqrt(Ew2 - Sqq);
92  const double Pw2 = Pw*Pw;
93  const double Mb2Mc2Sq = pow(Mb2 - Mc2, 2.);
94 
95 
96  double RateNorm = (pow(GFermi,2.)*Pw*pow(-Mt2 + Sqq,2.))/(32.*Mb2*pi3*Sqq);
97 
98  // set non-zero tensor elements
99  result.element({0, 1, 0, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
100  result.element({0, 1, 1, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
101  result.element({0, 1, 3, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
102  result.element({0, 1, 5, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
103  result.element({0, 1, 7, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
104  result.element({0, 2, 0, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
105  result.element({0, 2, 5, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
106  result.element({0, 2, 7, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
107  result.element({0, 2, 9, 3}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
108  result.element({1, 0, 0, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
109  result.element({1, 0, 1, 0}) = 0.5;
110  result.element({1, 0, 3, 0}) = 0.5;
111  result.element({1, 0, 5, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
112  result.element({1, 0, 7, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
113  result.element({2, 0, 2, 0}) = 0.5;
114  result.element({2, 0, 4, 0}) = 0.5;
115  result.element({2, 0, 6, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
116  result.element({2, 0, 8, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
117  result.element({3, 0, 0, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
118  result.element({3, 0, 1, 0}) = 0.5;
119  result.element({3, 0, 3, 0}) = 0.5;
120  result.element({3, 0, 5, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
121  result.element({3, 0, 7, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
122  result.element({4, 0, 2, 0}) = 0.5;
123  result.element({4, 0, 4, 0}) = 0.5;
124  result.element({4, 0, 6, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
125  result.element({4, 0, 8, 1}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
126  result.element({5, 1, 0, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
127  result.element({5, 1, 1, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
128  result.element({5, 1, 3, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
129  result.element({5, 1, 5, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
130  result.element({5, 1, 7, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
131  result.element({5, 2, 0, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
132  result.element({5, 2, 5, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
133  result.element({5, 2, 7, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
134  result.element({5, 2, 9, 3}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
135  result.element({6, 1, 2, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
136  result.element({6, 1, 4, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
137  result.element({6, 1, 6, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
138  result.element({6, 1, 8, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
139  result.element({6, 2, 6, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
140  result.element({6, 2, 8, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
141  result.element({6, 2, 10, 3}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
142  result.element({7, 1, 0, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
143  result.element({7, 1, 1, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
144  result.element({7, 1, 3, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
145  result.element({7, 1, 5, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
146  result.element({7, 1, 7, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
147  result.element({7, 2, 0, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
148  result.element({7, 2, 5, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
149  result.element({7, 2, 7, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
150  result.element({7, 2, 9, 3}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
151  result.element({8, 1, 2, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
152  result.element({8, 1, 4, 0}) = ((-Mb2 + Mc2) * Mt) / (2. * Sqq);
153  result.element({8, 1, 6, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
154  result.element({8, 1, 8, 1}) = (Mb2Mc2Sq * Mt2) / (2. * Sqq2);
155  result.element({8, 2, 6, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
156  result.element({8, 2, 8, 2}) = (2. * Mb2 * Pw2 * (Mt2 + 2. * Sqq)) / (3. * Sqq2);
157  result.element({8, 2, 10, 3}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
158  result.element({9, 3, 0, 2}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
159  result.element({9, 3, 5, 2}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
160  result.element({9, 3, 7, 2}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
161  result.element({9, 3, 9, 3}) = (32. * Mb2 * Pw2 * (2. * Mt2 + Sqq)) / (3. * Sqq);
162  result.element({10, 3, 6, 2}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
163  result.element({10, 3, 8, 2}) = (8. * Mb2 * Mt * (-Ew2 + Sqq)) / Sqq;
164  result.element({10, 3, 10, 3}) = (32. * Mb2 * Pw2 * (2. * Mt2 + Sqq)) / (3. * Sqq);
165 
166  result *= (RateNorm);
167 
168  return result;
169  }
170 
171 } // namespace Hammer
Hammer::MultiDimensional::makeEmptySparse
TensorData makeEmptySparse(const IndexList &dimensions, const LabelsList &labels)
Definition: SparseContainer.cc:287
Hammer::Tensor::labels
LabelsList labels() const
get the labels of all the indices at once
Definition: Tensor.cc:83
Hammer::GFermi
static constexpr double GFermi
Definition: Constants.hh:29
IndexLabels.hh
Tensor indices label definitions.
Hammer::pi3
static constexpr double pi3
Definition: Constants.hh:25
SparseContainer.hh
Sparse tensor data container.
Hammer::Tensor
Multidimensional tensor class with complex numbers as elements.
Definition: Tensor.hh:33
Constants.hh
Various numerical constants.
Pdg.hh
Hammer particle data class.
RateBDLepNu.hh
total rate
Hammer::PS::makeQ2Function
BoundaryFunction makeQ2Function(double qmin, double qmax)
Definition: RateBase.hh:122
Particle.hh
Hammer particle class.