Hammer  1.0.0
Helicity Amplitude Module for Matrix Element Reweighting
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AmplBD1starLepNu.cc
Go to the documentation of this file.
1 ///
2 /// @file AmplBD1starLepNu.cc
3 /// @brief \f$ B \rightarrow D_1^* \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/AmplBD1starLepNu.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  AmplBD1starLepNu::AmplBD1starLepNu() {
28  // Create tensor rank and dimensions
29  vector<IndexType> dims{{11, 8, 3, 2, 2, 2}};
30  string name{"AmplBD1starLepNu"};
31  addProcessSignature(PID::BPLUS, {-PID::DSSD1STAR, PID::NU_TAU, PID::ANTITAU});
32  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCTAUNU, FF_BDSSD1STAR, SPIN_DSSD1STAR, SPIN_NUTAU, SPIN_NUTAU_REF, SPIN_TAUP})});
33 
34  addProcessSignature(PID::BZERO, {PID::DSSD1STARMINUS, PID::NU_TAU, PID::ANTITAU});
35  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCTAUNU, FF_BDSSD1STAR, SPIN_DSSD1STAR, SPIN_NUTAU, SPIN_NUTAU_REF, SPIN_TAUP})});
36 
37  addProcessSignature(PID::BPLUS, {-PID::DSSD1STAR, PID::NU_MU, PID::ANTIMUON});
38  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCMUNU, FF_BDSSD1STAR, SPIN_DSSD1STAR, SPIN_NUMU, SPIN_NUMU_REF, SPIN_MUP})});
39 
40  addProcessSignature(PID::BZERO, {PID::DSSD1STARMINUS, PID::NU_MU, PID::ANTIMUON});
41  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCMUNU, FF_BDSSD1STAR, SPIN_DSSD1STAR, SPIN_NUMU, SPIN_NUMU_REF, SPIN_MUP})});
42 
43  addProcessSignature(PID::BPLUS, {-PID::DSSD1STAR, PID::NU_E, PID::POSITRON});
44  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCENU, FF_BDSSD1STAR, SPIN_DSSD1STAR, SPIN_NUE, SPIN_NUE_REF, SPIN_EP})});
45 
46  addProcessSignature(PID::BZERO, {PID::DSSD1STARMINUS, PID::NU_E, PID::POSITRON});
47  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCENU, FF_BDSSD1STAR, SPIN_DSSD1STAR, SPIN_NUE, SPIN_NUE_REF, SPIN_EP})});
48 
49  //bs -> cs
50  addProcessSignature(PID::BS, {PID::DSSDS1STARMINUS, PID::NU_TAU, PID::ANTITAU});
51  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCTAUNU, FF_BSDSSDS1STAR, SPIN_DSSDS1STAR, SPIN_NUTAU, SPIN_NUTAU_REF, SPIN_TAUP})});
52 
53  addProcessSignature(PID::BS, {PID::DSSDS1STARMINUS, PID::NU_MU, PID::ANTIMUON});
54  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCMUNU, FF_BSDSSDS1STAR, SPIN_DSSDS1STAR, SPIN_NUMU, SPIN_NUMU_REF, SPIN_MUP})});
55 
56  addProcessSignature(PID::BS, {PID::DSSDS1STARMINUS, PID::NU_E, PID::POSITRON});
57  addTensor(Tensor{name, MD::makeEmptySparse(dims, {WILSON_BCENU, FF_BSDSSDS1STAR, SPIN_DSSDS1STAR, SPIN_NUE, SPIN_NUE_REF, SPIN_EP})});
58 
59  setSignatureIndex();
60  }
61 
62  void AmplBD1starLepNu::eval(const Particle& parent, const ParticleList& daughters,
63  const ParticleList&) {
64  // Momenta
65  const FourMomentum& pBmes = parent.momentum();
66  const FourMomentum& pD1starmes = 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 Mc = pD1starmes.mass();
74  const double Mc2 = Mc*Mc;
75  const double Mt = pTau.mass();
76  const double Mt2 = Mt*Mt;
77  const double Sqq = Mb2 + Mc2 - 2. * (pBmes * pD1starmes);
78  // const double sqSqq = sqrt(Sqq);
79  const double Ew = (Mb2 - Mc2 + Sqq) / (2 * Mb);
80  const double BNuTau = (pBmes * kNuTau);
81  const double NuTauQ = (pBmes * kNuTau) - (pD1starmes * kNuTau);
82  const double BQ = Mb2 - (pBmes * pD1starmes);
83 
84  const double w = (Mb2 + Mc2 - Sqq)/(2 * Mb * Mc);
85  const double rC = Mc/Mb;
86  const double rt = Mt/Mb;
87 
88  const double mSqq = Sqq/Mb2;
89  const double w2m1 = w*w - 1;
90  const double Sqw2m1 = sqrt(w2m1);
91  const double SqmSqq = sqrt(mSqq);
92  const double Sqw2m1OnmSqq = sqrt(w2m1/mSqq);
93  // const double SqmSqqw2m1 = sqrt(mSqq*w2m1);
94  // const double w2m132 = pow(w2m1,1.5);
95 
96  // Helicity Angles
97  const double CosTt = -((Ew * (Sqq * BNuTau - NuTauQ * BQ)) / (sqrt(Ew*Ew - Sqq) * NuTauQ * BQ));
98  const double SinTt = sqrt(1. - CosTt*CosTt);
99  const double CosTtHalfSq = (1. + CosTt) / 2.;
100  const double SinTtHalfSq = (1. - CosTt) / 2.;
101 
102  const double prefactor = 2*GFermi*sqrt(Mb*Mc)*sqrt(Sqq - Mt2);
103 
104  // initialize tensor elements to zero
105  Tensor& t = getTensor();
106  t.clearData();
107 
108  // set non-zero tensor elements
109  // NB: Amplitudes include additional D1star phase, defined wrt D1star spin!
110  t.element({0,1,0,0,0,0}) = (rt*SinTt)/(2.*SqmSqq);
111  t.element({0,4,0,0,0,0}) = (rt*SinTt*Sqw2m1OnmSqq)/2.;
112  t.element({5,1,0,0,0,0}) = (rt*SinTt)/(2.*SqmSqq);
113  t.element({5,4,0,0,0,0}) = -(rt*SinTt*Sqw2m1OnmSqq)/2.;
114  t.element({7,1,0,0,0,0}) = (rt*SinTt)/(2.*SqmSqq);
115  t.element({7,4,0,0,0,0}) = (rt*SinTt*Sqw2m1OnmSqq)/2.;
116  t.element({9,5,0,0,0,0}) = (2.*SinTt*(-1 + rC*(w + Sqw2m1)))/SqmSqq;
117  t.element({9,6,0,0,0,0}) = (2.*SinTt*(rC - w + Sqw2m1))/SqmSqq;
118  t.element({0,1,0,0,0,1}) = CosTtHalfSq;
119  t.element({0,4,0,0,0,1}) = CosTtHalfSq*Sqw2m1;
120  t.element({5,1,0,0,0,1}) = CosTtHalfSq;
121  t.element({5,4,0,0,0,1}) = -(CosTtHalfSq*Sqw2m1);
122  t.element({7,1,0,0,0,1}) = CosTtHalfSq;
123  t.element({7,4,0,0,0,1}) = CosTtHalfSq*Sqw2m1;
124  t.element({9,5,0,0,0,1}) = (4.*CosTtHalfSq*rt*(-1 + rC*(w + Sqw2m1)))/mSqq;
125  t.element({9,6,0,0,0,1}) = (4.*CosTtHalfSq*rt*(rC - w + Sqw2m1))/mSqq;
126  t.element({6,1,0,1,1,0}) = -SinTtHalfSq;
127  t.element({6,4,0,1,1,0}) = SinTtHalfSq*Sqw2m1;
128  t.element({8,1,0,1,1,0}) = -SinTtHalfSq;
129  t.element({8,4,0,1,1,0}) = -(SinTtHalfSq*Sqw2m1);
130  t.element({10,5,0,1,1,0}) = (4.*rt*SinTtHalfSq*(1 + rC*(-w + Sqw2m1)))/mSqq;
131  t.element({10,6,0,1,1,0}) = (4.*rt*SinTtHalfSq*(-rC + w + Sqw2m1))/mSqq;
132  t.element({6,1,0,1,1,1}) = -(rt*SinTt)/(2.*SqmSqq);
133  t.element({6,4,0,1,1,1}) = (rt*SinTt*Sqw2m1OnmSqq)/2.;
134  t.element({8,1,0,1,1,1}) = -(rt*SinTt)/(2.*SqmSqq);
135  t.element({8,4,0,1,1,1}) = -(rt*SinTt*Sqw2m1OnmSqq)/2.;
136  t.element({10,5,0,1,1,1}) = (2.*SinTt*(1 + rC*(-w + Sqw2m1)))/SqmSqq;
137  t.element({10,6,0,1,1,1}) = (2.*SinTt*(-rC + w + Sqw2m1))/SqmSqq;
138  t.element({0,1,1,0,0,0}) = (rt*(CosTt*(rC - w) + Sqw2m1))/(sqrt2*mSqq);
139  t.element({0,2,1,0,0,0}) = -((rt*((-1 + rC*w)*Sqw2m1 + CosTt*rC*w2m1))/(sqrt2*mSqq));
140  t.element({0,3,1,0,0,0}) = (rt*((-rC + w)*Sqw2m1 - CosTt*w2m1))/(sqrt2*mSqq);
141  t.element({1,0,1,0,0,0}) = Sqw2m1/sqrt2;
142  t.element({3,0,1,0,0,0}) = Sqw2m1/sqrt2;
143  t.element({5,1,1,0,0,0}) = (rt*(CosTt*(rC - w) + Sqw2m1))/(sqrt2*mSqq);
144  t.element({5,2,1,0,0,0}) = -((rt*((-1 + rC*w)*Sqw2m1 + CosTt*rC*w2m1))/(sqrt2*mSqq));
145  t.element({5,3,1,0,0,0}) = (rt*((-rC + w)*Sqw2m1 - CosTt*w2m1))/(sqrt2*mSqq);
146  t.element({7,1,1,0,0,0}) = (rt*(CosTt*(rC - w) + Sqw2m1))/(sqrt2*mSqq);
147  t.element({7,2,1,0,0,0}) = -((rt*((-1 + rC*w)*Sqw2m1 + CosTt*rC*w2m1))/(sqrt2*mSqq));
148  t.element({7,3,1,0,0,0}) = (rt*((-rC + w)*Sqw2m1 - CosTt*w2m1))/(sqrt2*mSqq);
149  t.element({9,5,1,0,0,0}) = 2.*sqrt2*CosTt*w;
150  t.element({9,6,1,0,0,0}) = 2.*sqrt2*CosTt;
151  t.element({9,7,1,0,0,0}) = -2.*sqrt2*CosTt*w2m1;
152  t.element({0,1,1,0,0,1}) = (SinTt*(-rC + w))/(sqrt2*SqmSqq);
153  t.element({0,2,1,0,0,1}) = (rC*SinTt*w2m1)/(sqrt2*SqmSqq);
154  t.element({0,3,1,0,0,1}) = (SinTt*w2m1)/(sqrt2*SqmSqq);
155  t.element({5,1,1,0,0,1}) = (SinTt*(-rC + w))/(sqrt2*SqmSqq);
156  t.element({5,2,1,0,0,1}) = (rC*SinTt*w2m1)/(sqrt2*SqmSqq);
157  t.element({5,3,1,0,0,1}) = (SinTt*w2m1)/(sqrt2*SqmSqq);
158  t.element({7,1,1,0,0,1}) = (SinTt*(-rC + w))/(sqrt2*SqmSqq);
159  t.element({7,2,1,0,0,1}) = (rC*SinTt*w2m1)/(sqrt2*SqmSqq);
160  t.element({7,3,1,0,0,1}) = (SinTt*w2m1)/(sqrt2*SqmSqq);
161  t.element({9,5,1,0,0,1}) = (-2.*sqrt2*rt*SinTt*w)/SqmSqq;
162  t.element({9,6,1,0,0,1}) = (-2.*sqrt2*rt*SinTt)/SqmSqq;
163  t.element({9,7,1,0,0,1}) = (2.*sqrt2*rt*SinTt*w2m1)/SqmSqq;
164  t.element({6,1,1,1,1,0}) = (SinTt*(-rC + w))/(sqrt2*SqmSqq);
165  t.element({6,2,1,1,1,0}) = (rC*SinTt*w2m1)/(sqrt2*SqmSqq);
166  t.element({6,3,1,1,1,0}) = (SinTt*w2m1)/(sqrt2*SqmSqq);
167  t.element({8,1,1,1,1,0}) = (SinTt*(-rC + w))/(sqrt2*SqmSqq);
168  t.element({8,2,1,1,1,0}) = (rC*SinTt*w2m1)/(sqrt2*SqmSqq);
169  t.element({8,3,1,1,1,0}) = (SinTt*w2m1)/(sqrt2*SqmSqq);
170  t.element({10,5,1,1,1,0}) = (-2.*sqrt2*rt*SinTt*w)/SqmSqq;
171  t.element({10,6,1,1,1,0}) = (-2.*sqrt2*rt*SinTt)/SqmSqq;
172  t.element({10,7,1,1,1,0}) = (2.*sqrt2*rt*SinTt*w2m1)/SqmSqq;
173  t.element({2,0,1,1,1,1}) = -(Sqw2m1/sqrt2);
174  t.element({4,0,1,1,1,1}) = -(Sqw2m1/sqrt2);
175  t.element({6,1,1,1,1,1}) = -((rt*(CosTt*(rC - w) + Sqw2m1))/(sqrt2*mSqq));
176  t.element({6,2,1,1,1,1}) = (rt*((-1 + rC*w)*Sqw2m1 + CosTt*rC*w2m1))/(sqrt2*mSqq);
177  t.element({6,3,1,1,1,1}) = (rt*((rC - w)*Sqw2m1 + CosTt*w2m1))/(sqrt2*mSqq);
178  t.element({8,1,1,1,1,1}) = -((rt*(CosTt*(rC - w) + Sqw2m1))/(sqrt2*mSqq));
179  t.element({8,2,1,1,1,1}) = (rt*((-1 + rC*w)*Sqw2m1 + CosTt*rC*w2m1))/(sqrt2*mSqq);
180  t.element({8,3,1,1,1,1}) = (rt*((rC - w)*Sqw2m1 + CosTt*w2m1))/(sqrt2*mSqq);
181  t.element({10,5,1,1,1,1}) = -2.*sqrt2*CosTt*w;
182  t.element({10,6,1,1,1,1}) = -2.*sqrt2*CosTt;
183  t.element({10,7,1,1,1,1}) = 2.*sqrt2*CosTt*w2m1;
184  t.element({0,1,2,0,0,0}) = (rt*SinTt)/(2.*SqmSqq);
185  t.element({0,4,2,0,0,0}) = -(rt*SinTt*Sqw2m1OnmSqq)/2.;
186  t.element({5,1,2,0,0,0}) = (rt*SinTt)/(2.*SqmSqq);
187  t.element({5,4,2,0,0,0}) = (rt*SinTt*Sqw2m1OnmSqq)/2.;
188  t.element({7,1,2,0,0,0}) = (rt*SinTt)/(2.*SqmSqq);
189  t.element({7,4,2,0,0,0}) = -(rt*SinTt*Sqw2m1OnmSqq)/2.;
190  t.element({9,5,2,0,0,0}) = (-2.*SinTt*(1 + rC*(-w + Sqw2m1)))/SqmSqq;
191  t.element({9,6,2,0,0,0}) = (-2.*SinTt*(-rC + w + Sqw2m1))/SqmSqq;
192  t.element({0,1,2,0,0,1}) = -SinTtHalfSq;
193  t.element({0,4,2,0,0,1}) = SinTtHalfSq*Sqw2m1;
194  t.element({5,1,2,0,0,1}) = -SinTtHalfSq;
195  t.element({5,4,2,0,0,1}) = -(SinTtHalfSq*Sqw2m1);
196  t.element({7,1,2,0,0,1}) = -SinTtHalfSq;
197  t.element({7,4,2,0,0,1}) = SinTtHalfSq*Sqw2m1;
198  t.element({9,5,2,0,0,1}) = (4.*rt*SinTtHalfSq*(1 + rC*(-w + Sqw2m1)))/mSqq;
199  t.element({9,6,2,0,0,1}) = (4.*rt*SinTtHalfSq*(-rC + w + Sqw2m1))/mSqq;
200  t.element({6,1,2,1,1,0}) = CosTtHalfSq;
201  t.element({6,4,2,1,1,0}) = CosTtHalfSq*Sqw2m1;
202  t.element({8,1,2,1,1,0}) = CosTtHalfSq;
203  t.element({8,4,2,1,1,0}) = -(CosTtHalfSq*Sqw2m1);
204  t.element({10,5,2,1,1,0}) = (4.*CosTtHalfSq*rt*(-1 + rC*(w + Sqw2m1)))/mSqq;
205  t.element({10,6,2,1,1,0}) = (4.*CosTtHalfSq*rt*(rC - w + Sqw2m1))/mSqq;
206  t.element({6,1,2,1,1,1}) = -(rt*SinTt)/(2.*SqmSqq);
207  t.element({6,4,2,1,1,1}) = -(rt*SinTt*Sqw2m1OnmSqq)/2.;
208  t.element({8,1,2,1,1,1}) = -(rt*SinTt)/(2.*SqmSqq);
209  t.element({8,4,2,1,1,1}) = (rt*SinTt*Sqw2m1OnmSqq)/2.;
210  t.element({10,5,2,1,1,1}) = (-2.*SinTt*(-1 + rC*(w + Sqw2m1)))/SqmSqq;
211  t.element({10,6,2,1,1,1}) = (-2.*SinTt*(rC - w + Sqw2m1))/SqmSqq;
212 
213  t *= prefactor;
214  }
215 
216  void AmplBD1starLepNu::addRefs() const {
217  if(!getSettingsHandler()->checkReference("Bernlochner:2017jxt")){
218  string ref =
219  "@article{Bernlochner:2017jxt,\n"
220  " author = \"Bernlochner, Florian U. and Ligeti, Zoltan and Robinson, Dean J.\",\n"
221  " title = \"{Model independent analysis of semileptonic $B$ decays to $D^{**}$ for arbitrary new physics}\",\n"
222  " journal = \"Phys. Rev.\",\n"
223  " volume = \"D97\",\n"
224  " year = \"2018\",\n"
225  " number = \"7\",\n"
226  " pages = \"075011\",\n"
227  " doi = \"10.1103/PhysRevD.97.075011\",\n"
228  " eprint = \"1711.03110\",\n"
229  " archivePrefix = \"arXiv\",\n"
230  " primaryClass = \"hep-ph\",\n"
231  " SLACcitation = \"%%CITATION = ARXIV:1711.03110;%%\"\n"
232  "}\n";
233  getSettingsHandler()->addReference("Bernlochner:2017jxt", ref);
234  }
235  }
236 
237 } // 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