Hammer  1.0.0
Helicity Amplitude Module for Matrix Element Reweighting
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AmplBDstarLepNu.cc
Go to the documentation of this file.
1 ///
2 /// @file AmplBDstarLepNu.cc
3 /// @brief \f$ B \rightarrow D^* \tau\nu \f$ amplitude
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/Amplitudes/AmplBDstarLepNu.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 using namespace complex_literals;
22 
23 namespace Hammer {
24 
25  namespace MD = MultiDimensional;
26 
27  AmplBDstarLepNu::AmplBDstarLepNu() {
28  // Create tensor rank and dimensions
29  vector<IndexType> dims{{11, 8, 3, 2, 2, 2}};
30  string name{"AmplBDstarLepNu"};
31  addProcessSignature(PID::BPLUS, {-PID::DSTAR, PID::NU_TAU, PID::ANTITAU});
32  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCTAUNU, FF_BDSTAR, SPIN_DSTAR, SPIN_NUTAU, SPIN_NUTAU_REF, SPIN_TAUP})});
33 
34  addProcessSignature(PID::BZERO, {PID::DSTARMINUS, PID::NU_TAU, PID::ANTITAU});
35  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCTAUNU, FF_BDSTAR, SPIN_DSTAR, SPIN_NUTAU, SPIN_NUTAU_REF, SPIN_TAUP})});
36 
37  addProcessSignature(PID::BPLUS, {-PID::DSTAR, PID::NU_MU, PID::ANTIMUON});
38  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCMUNU, FF_BDSTAR, SPIN_DSTAR, SPIN_NUMU, SPIN_NUMU_REF, SPIN_MUP})});
39 
40  addProcessSignature(PID::BZERO, {PID::DSTARMINUS, PID::NU_MU, PID::ANTIMUON});
41  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCMUNU, FF_BDSTAR, SPIN_DSTAR, SPIN_NUMU, SPIN_NUMU_REF, SPIN_MUP})});
42 
43  addProcessSignature(PID::BPLUS, {-PID::DSTAR, PID::NU_E, PID::POSITRON});
44  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCENU, FF_BDSTAR, SPIN_DSTAR, SPIN_NUE, SPIN_NUE_REF, SPIN_EP})});
45 
46  addProcessSignature(PID::BZERO, {PID::DSTARMINUS, PID::NU_E, PID::POSITRON});
47  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCENU, FF_BDSTAR, SPIN_DSTAR, SPIN_NUE, SPIN_NUE_REF, SPIN_EP})});
48 
49  //bs -> cs
50  addProcessSignature(PID::BS, {PID::DSSTARMINUS, PID::NU_TAU, PID::ANTITAU});
51  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCTAUNU, FF_BSDSSTAR, SPIN_DSSTAR, SPIN_NUTAU, SPIN_NUTAU_REF, SPIN_TAUP})});
52 
53  addProcessSignature(PID::BS, {PID::DSSTARMINUS, PID::NU_MU, PID::ANTIMUON});
54  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCMUNU, FF_BSDSSTAR, SPIN_DSSTAR, SPIN_NUMU, SPIN_NUMU_REF, SPIN_MUP})});
55 
56  addProcessSignature(PID::BS, {PID::DSSTARMINUS, PID::NU_E, PID::POSITRON});
57  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCENU, FF_BSDSSTAR, SPIN_DSSTAR, SPIN_NUE, SPIN_NUE_REF, SPIN_EP})});
58 
59  setSignatureIndex();
60  }
61 
62  void AmplBDstarLepNu::eval(const Particle& parent, const ParticleList& daughters,
63  const ParticleList&) {
64  // Momenta
65  const FourMomentum& pBmes = parent.momentum();
66  const FourMomentum& pDstarmes = daughters[0].momentum();
67  const FourMomentum& kNuTau = daughters[1].momentum();
68  const FourMomentum& pTau = daughters[2].momentum();
69 
70  // kinematic objects
71  const double Mb = pBmes.mass();
72  const double Mb2 = Mb*Mb;
73  const double Mds = pDstarmes.mass();
74  const double Mds2 = Mds*Mds;
75  const double Mt = pTau.mass();
76  const double Sqq = Mb2 + Mds2 - 2. * (pBmes * pDstarmes);
77  const double sqSqq = sqrt(Sqq);
78  const double Ew = (Mb2 - Mds2 + Sqq) / (2 * Mb);
79  // const double Eds = (Mb2 + Mds2 - Sqq) / (2 * Mb);
80  const double Pw = sqrt(Ew*Ew - Sqq);
81  const double Pw2 = Pw*Pw;
82  const double BNuTau = (pBmes * kNuTau);
83  const double NuTauQ = (pBmes * kNuTau) - (pDstarmes * kNuTau);
84  const double BQ = Mb2 - (pBmes * pDstarmes);
85 
86  // Helicity Angles
87  const double CosTt = -((Ew * (Sqq * BNuTau - NuTauQ * BQ)) / (sqrt(Ew*Ew - Sqq) * NuTauQ * BQ));
88  const double SinTt = sqrt(1. - CosTt*CosTt);
89  const double CosTtHalfSq = (1. + CosTt) / 2.;
90  const double SinTtHalfSq = (1. - CosTt) / 2.;
91 
92  const double TwoGfSq = 2. * GFermi * sqrt(Sqq - Mt*Mt);
93 
94  // initialize tensor elements to zero
95  Tensor& t = getTensor();
96  t.clearData();
97 
98  // set non-zero tensor elements
99  // NB: Amplitudes include additional Dstar phase, defined wrt Dstar spin, that removes tau phase. Must be include in D* decay amplitudes.
100  t.element({0, 1, 0, 0, 0, 0}) = -(Mt * SinTt)/(2. * sqSqq);
101  t.element({0, 2, 0, 0, 0, 0}) = -((Mb * Mt * Pw * SinTt)/sqSqq);
102  t.element({5, 1, 0, 0, 0, 0}) = (Mt * SinTt)/(2. * sqSqq);
103  t.element({5, 2, 0, 0, 0, 0}) = -((Mb * Mt * Pw * SinTt)/sqSqq);
104  t.element({7, 1, 0, 0, 0, 0}) = -(Mt * SinTt)/(2. * sqSqq);
105  t.element({7, 2, 0, 0, 0, 0}) = -((Mb * Mt * Pw * SinTt)/sqSqq);
106  t.element({9, 6, 0, 0, 0, 0}) = 2. * sqSqq * SinTt;
107  t.element({9, 7, 0, 0, 0, 0}) = (2. * (Mb2 - Mds2 - 2. * Mb * Pw) * SinTt)/sqSqq;
108  t.element({0, 1, 0, 0, 0, 1}) = -CosTtHalfSq;
109  t.element({0, 2, 0, 0, 0, 1}) = -2. * Mb * Pw * CosTtHalfSq;
110  t.element({5, 1, 0, 0, 0, 1}) = CosTtHalfSq;
111  t.element({5, 2, 0, 0, 0, 1}) = -2. * Mb * Pw * CosTtHalfSq;
112  t.element({7, 1, 0, 0, 0, 1}) = -CosTtHalfSq;
113  t.element({7, 2, 0, 0, 0, 1}) = -2. * Mb * Pw * CosTtHalfSq;
114  t.element({9, 6, 0, 0, 0, 1}) = 4. * Mt * CosTtHalfSq;
115  t.element({9, 7, 0, 0, 0, 1}) = (4. * Mt * (Mb2 - Mds2 - 2. * Mb * Pw) * CosTtHalfSq)/Sqq;
116  t.element({6, 1, 0, 1, 1, 0}) = -SinTtHalfSq;
117  t.element({6, 2, 0, 1, 1, 0}) = 2. * Mb * Pw * SinTtHalfSq;
118  t.element({8, 1, 0, 1, 1, 0}) = SinTtHalfSq;
119  t.element({8, 2, 0, 1, 1, 0}) = 2. * Mb * Pw * SinTtHalfSq;
120  t.element({10, 6, 0, 1, 1, 0}) = 4. * Mt * SinTtHalfSq;
121  t.element({10, 7, 0, 1, 1, 0}) = (4. * Mt * (Mb2 - Mds2 + 2. * Mb * Pw) * SinTtHalfSq)/Sqq;
122  t.element({6, 1, 0, 1, 1, 1}) = -(Mt * SinTt)/(2. * sqSqq);
123  t.element({6, 2, 0, 1, 1, 1}) = (Mb * Mt * Pw * SinTt)/sqSqq;
124  t.element({8, 1, 0, 1, 1, 1}) = (Mt * SinTt)/(2. * sqSqq);
125  t.element({8, 2, 0, 1, 1, 1}) = (Mb * Mt * Pw * SinTt)/sqSqq;
126  t.element({10, 6, 0, 1, 1, 1}) = 2. * sqSqq * SinTt;
127  t.element({10, 7, 0, 1, 1, 1}) = (2. * (Mb2 - Mds2 + 2. * Mb * Pw) * SinTt)/sqSqq;
128  t.element({0, 1, 1, 0, 0, 0}) = -(Mt * (2. * Mb * Pw + (-Mb2 + Mds2 + Sqq) * CosTt))/(2. * sqrt2 * Mds * Sqq);
129  t.element({0, 3, 1, 0, 0, 0}) = -((Mb * Mt * Pw)/(sqrt2 * Mds));
130  t.element({0, 4, 1, 0, 0, 0}) = (Mb * Mt * Pw * (-Mb2 + Mds2 + 2. * Mb * Pw * CosTt))/(sqrt2 * Mds * Sqq);
131  t.element({1, 0, 1, 0, 0, 0}) = (Mb * Pw)/(sqrt2 * Mds);
132  t.element({3, 0, 1, 0, 0, 0}) = -((Mb * Pw)/(sqrt2 * Mds));
133  t.element({5, 1, 1, 0, 0, 0}) = (Mt * (2. * Mb * Pw + (-Mb2 + Mds2 + Sqq) * CosTt))/(2. * sqrt2 * Mds * Sqq);
134  t.element({5, 3, 1, 0, 0, 0}) = (Mb * Mt * Pw)/(sqrt2 * Mds);
135  t.element({5, 4, 1, 0, 0, 0}) = -((Mb * Mt * Pw * (-Mb2 + Mds2 + 2. * Mb * Pw * CosTt))/(sqrt2 * Mds * Sqq));
136  t.element({7, 1, 1, 0, 0, 0}) = -(Mt * (2. * Mb * Pw + (-Mb2 + Mds2 + Sqq) * CosTt))/(2. * sqrt2 * Mds * Sqq);
137  t.element({7, 3, 1, 0, 0, 0}) = -((Mb * Mt * Pw)/(sqrt2 * Mds));
138  t.element({7, 4, 1, 0, 0, 0}) = (Mb * Mt * Pw * (-Mb2 + Mds2 + 2. * Mb * Pw * CosTt))/(sqrt2 * Mds * Sqq);
139  t.element({9, 5, 1, 0, 0, 0}) = (-4. * sqrt2 * Mb2 * Pw2 * CosTt)/Mds;
140  t.element({9, 6, 1, 0, 0, 0}) = (sqrt2 * (-Mb2 + Mds2 + Sqq) * CosTt)/Mds;
141  t.element({9, 7, 1, 0, 0, 0}) = (sqrt2 * (-Mb2 - 3. * Mds2 + Sqq) * CosTt)/Mds;
142  t.element({0, 1, 1, 0, 0, 1}) = ((-Mb2 + Mds2 + Sqq) * SinTt)/(2. * sqrt2 * Mds * sqSqq);
143  t.element({0, 4, 1, 0, 0, 1}) = -((sqrt2 * Mb2 * Pw2 * SinTt)/(Mds * sqSqq));
144  t.element({5, 1, 1, 0, 0, 1}) = -((-Mb2 + Mds2 + Sqq) * SinTt)/(2. * sqrt2 * Mds * sqSqq);
145  t.element({5, 4, 1, 0, 0, 1}) = (sqrt2 * Mb2 * Pw2 * SinTt)/(Mds * sqSqq);
146  t.element({7, 1, 1, 0, 0, 1}) = ((-Mb2 + Mds2 + Sqq) * SinTt)/(2. * sqrt2 * Mds * sqSqq);
147  t.element({7, 4, 1, 0, 0, 1}) = -((sqrt2 * Mb2 * Pw2 * SinTt)/(Mds * sqSqq));
148  t.element({9, 5, 1, 0, 0, 1}) = (4. * sqrt2 * Mb2 * Mt * Pw2 * SinTt)/(Mds * sqSqq);
149  t.element({9, 6, 1, 0, 0, 1}) = -((sqrt2 * Mt * (-Mb2 + Mds2 + Sqq) * SinTt)/(Mds * sqSqq));
150  t.element({9, 7, 1, 0, 0, 1}) = (sqrt2 * Mt * (Mb2 + 3. * Mds2 - Sqq) * SinTt)/(Mds * sqSqq);
151  t.element({6, 1, 1, 1, 1, 0}) = -((-Mb2 + Mds2 + Sqq) * SinTt)/(2. * sqrt2 * Mds * sqSqq);
152  t.element({6, 4, 1, 1, 1, 0}) = (sqrt2 * Mb2 * Pw2 * SinTt)/(Mds * sqSqq);
153  t.element({8, 1, 1, 1, 1, 0}) = ((-Mb2 + Mds2 + Sqq) * SinTt)/(2. * sqrt2 * Mds * sqSqq);
154  t.element({8, 4, 1, 1, 1, 0}) = -((sqrt2 * Mb2 * Pw2 * SinTt)/(Mds * sqSqq));
155  t.element({10, 5, 1, 1, 1, 0}) = (-4. * sqrt2 * Mb2 * Mt * Pw2 * SinTt)/(Mds * sqSqq);
156  t.element({10, 6, 1, 1, 1, 0}) = (sqrt2 * Mt * (-Mb2 + Mds2 + Sqq) * SinTt)/(Mds * sqSqq);
157  t.element({10, 7, 1, 1, 1, 0}) = (sqrt2 * Mt * (-Mb2 - 3. * Mds2 + Sqq) * SinTt)/(Mds * sqSqq);
158  t.element({2, 0, 1, 1, 1, 1}) = -((Mb * Pw)/(sqrt2 * Mds));
159  t.element({4, 0, 1, 1, 1, 1}) = (Mb * Pw)/(sqrt2 * Mds);
160  t.element({6, 1, 1, 1, 1, 1}) = -(Mt * (2. * Mb * Pw + (-Mb2 + Mds2 + Sqq) * CosTt))/(2. * sqrt2 * Mds * Sqq);
161  t.element({6, 3, 1, 1, 1, 1}) = -((Mb * Mt * Pw)/(sqrt2 * Mds));
162  t.element({6, 4, 1, 1, 1, 1}) = (Mb * Mt * Pw * (-Mb2 + Mds2 + 2. * Mb * Pw * CosTt))/(sqrt2 * Mds * Sqq);
163  t.element({8, 1, 1, 1, 1, 1}) = (Mt * (2. * Mb * Pw + (-Mb2 + Mds2 + Sqq) * CosTt))/(2. * sqrt2 * Mds * Sqq);
164  t.element({8, 3, 1, 1, 1, 1}) = (Mb * Mt * Pw)/(sqrt2 * Mds);
165  t.element({8, 4, 1, 1, 1, 1}) = -((Mb * Mt * Pw * (-Mb2 + Mds2 + 2. * Mb * Pw * CosTt))/(sqrt2 * Mds * Sqq));
166  t.element({10, 5, 1, 1, 1, 1}) = (-4. * sqrt2 * Mb2 * Pw2 * CosTt)/Mds;
167  t.element({10, 6, 1, 1, 1, 1}) = (sqrt2 * (-Mb2 + Mds2 + Sqq) * CosTt)/Mds;
168  t.element({10, 7, 1, 1, 1, 1}) = (sqrt2 * (-Mb2 - 3. * Mds2 + Sqq) * CosTt)/Mds;
169  t.element({0, 1, 2, 0, 0, 0}) = -(Mt * SinTt)/(2. * sqSqq);
170  t.element({0, 2, 2, 0, 0, 0}) = (Mb * Mt * Pw * SinTt)/sqSqq;
171  t.element({5, 1, 2, 0, 0, 0}) = (Mt * SinTt)/(2. * sqSqq);
172  t.element({5, 2, 2, 0, 0, 0}) = (Mb * Mt * Pw * SinTt)/sqSqq;
173  t.element({7, 1, 2, 0, 0, 0}) = -(Mt * SinTt)/(2. * sqSqq);
174  t.element({7, 2, 2, 0, 0, 0}) = (Mb * Mt * Pw * SinTt)/sqSqq;
175  t.element({9, 6, 2, 0, 0, 0}) = 2. * sqSqq * SinTt;
176  t.element({9, 7, 2, 0, 0, 0}) = (2. * (Mb2 - Mds2 + 2. * Mb * Pw) * SinTt)/sqSqq;
177  t.element({0, 1, 2, 0, 0, 1}) = SinTtHalfSq;
178  t.element({0, 2, 2, 0, 0, 1}) = -2. * Mb * Pw * SinTtHalfSq;
179  t.element({5, 1, 2, 0, 0, 1}) = -SinTtHalfSq;
180  t.element({5, 2, 2, 0, 0, 1}) = -2. * Mb * Pw * SinTtHalfSq;
181  t.element({7, 1, 2, 0, 0, 1}) = SinTtHalfSq;
182  t.element({7, 2, 2, 0, 0, 1}) = -2. * Mb * Pw * SinTtHalfSq;
183  t.element({9, 6, 2, 0, 0, 1}) = -4. * Mt * SinTtHalfSq;
184  t.element({9, 7, 2, 0, 0, 1}) = (4. * Mt * (Mds2 - Mb * (Mb + 2. * Pw)) * SinTtHalfSq)/Sqq;
185  t.element({6, 1, 2, 1, 1, 0}) = CosTtHalfSq;
186  t.element({6, 2, 2, 1, 1, 0}) = 2. * Mb * Pw * CosTtHalfSq;
187  t.element({8, 1, 2, 1, 1, 0}) = -CosTtHalfSq;
188  t.element({8, 2, 2, 1, 1, 0}) = 2. * Mb * Pw * CosTtHalfSq;
189  t.element({10, 6, 2, 1, 1, 0}) = -4. * Mt * CosTtHalfSq;
190  t.element({10, 7, 2, 1, 1, 0}) = (4. * Mt * (-Mb2 + Mds2 + 2. * Mb * Pw) * CosTtHalfSq)/Sqq;
191  t.element({6, 1, 2, 1, 1, 1}) = -(Mt * SinTt)/(2. * sqSqq);
192  t.element({6, 2, 2, 1, 1, 1}) = -((Mb * Mt * Pw * SinTt)/sqSqq);
193  t.element({8, 1, 2, 1, 1, 1}) = (Mt * SinTt)/(2. * sqSqq);
194  t.element({8, 2, 2, 1, 1, 1}) = -((Mb * Mt * Pw * SinTt)/sqSqq);
195  t.element({10, 6, 2, 1, 1, 1}) = 2. * sqSqq * SinTt;
196  t.element({10, 7, 2, 1, 1, 1}) = (2. * (Mb2 - Mds2 - 2. * Mb * Pw) * SinTt)/sqSqq;
197 
198  t *= TwoGfSq;
199  }
200 
201  void AmplBDstarLepNu::addRefs() const {
202  if(!getSettingsHandler()->checkReference("Ligeti:2016npd")){
203  string ref =
204  "@article{Ligeti:2016npd,\n"
205  " author = \"Ligeti, Zoltan and Papucci, Michele and Robinson, Dean J.\",\n"
206  " title = \"{New Physics in the Visible Final States of $B\\to D^{(*)}\\tau\\nu$}\",\n"
207  " journal = \"JHEP\",\n"
208  " volume = \"01\",\n"
209  " year = \"2017\",\n"
210  " pages = \"083\",\n"
211  " doi = \"10.1007/JHEP01(2017)083\",\n"
212  " eprint = \"1610.02045\",\n"
213  " archivePrefix = \"arXiv\",\n"
214  " primaryClass = \"hep-ph\",\n"
215  " SLACcitation = \"%%CITATION = ARXIV:1610.02045;%%\"\n"
216  "}\n";
217  getSettingsHandler()->addReference("Ligeti:2016npd", ref);
218  }
219  }
220 
221 
222 } // namespace Hammer
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::GFermi
static constexpr double GFermi
Definition: Constants.hh:29
IndexLabels.hh
Tensor indices label definitions.
Hammer::sqrt2
static constexpr double sqrt2
Definition: Constants.hh:26
Hammer::FourMomentum::mass
double mass() const
returns the invariant mass if the invariant mass squared is negative returns
Hammer::ParticleList
std::vector< Particle > ParticleList
Definition: Particle.fhh:20
SparseContainer.hh
Sparse tensor data container.
Hammer::Particle
Particle class.
Definition: Particle.hh:30
Hammer::Tensor
Multidimensional tensor class with complex numbers as elements.
Definition: Tensor.hh:33
Hammer::Tensor::clearData
void clearData()
sets all the elements to 0
Definition: Tensor.cc:229
Hammer::Particle::momentum
const FourMomentum & momentum() const
Definition: Particle.cc:56
Constants.hh
Various numerical constants.
Pdg.hh
Hammer particle data class.
Particle.hh
Hammer particle class.
Hammer::FourMomentum
4-momentum class
Definition: FourMomentum.hh:30