evtgen is hosted by Hepforge, IPPP Durham
EvtGen  2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
EvtbTosllVectorAmp.cpp
Go to the documentation of this file.
1 
2 /***********************************************************************
3 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4 * *
5 * This file is part of EvtGen. *
6 * *
7 * EvtGen is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * EvtGen is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19 ***********************************************************************/
20 
22 
23 #include "EvtGenBase/EvtAmp.hh"
25 #include "EvtGenBase/EvtGenKine.hh"
26 #include "EvtGenBase/EvtId.hh"
27 #include "EvtGenBase/EvtIdSet.hh"
28 #include "EvtGenBase/EvtPDL.hh"
30 #include "EvtGenBase/EvtPatches.hh"
31 #include "EvtGenBase/EvtReport.hh"
34 
37 
39  EvtbTosllFF* formFactors )
40 {
41  //Add the lepton and neutrino 4 momenta to find q2
42 
43  EvtVector4R q = parent->getDaug( 1 )->getP4() + parent->getDaug( 2 )->getP4();
44  double q2 = ( q.mass2() );
45 
46  double a1, a2, a0, v, t1, t2, t3;
47  double mesonmass = parent->getDaug( 0 )->mass();
48  double parentmass = parent->mass();
49 
50  formFactors->getVectorFF( parent->getId(), parent->getDaug( 0 )->getId(),
51  q2, mesonmass, a1, a2, a0, v, t1, t2, t3 );
52 
53  EvtId daught = parent->getDaug( 0 )->getId();
54  EvtId parentId = parent->getId();
55  bool btod = false;
56  bool nnlo = true;
57  if ( ( parentId == EvtPDL::getId( "B0" ) ||
58  parentId == EvtPDL::getId( "anti-B0" ) ||
59  parentId == EvtPDL::getId( "B+" ) ||
60  parentId == EvtPDL::getId( "B-" ) ) &&
61  ( daught == EvtPDL::getId( std::string( "rho+" ) ) ||
62  daught == EvtPDL::getId( std::string( "rho-" ) ) ||
63  daught == EvtPDL::getId( std::string( "rho0" ) ) ||
64  daught == EvtPDL::getId( std::string( "omega" ) ) ) ) {
65  btod = true;
66  }
67  if ( ( parentId == EvtPDL::getId( "B_s0" ) ||
68  parentId == EvtPDL::getId( "anti-B_s0" ) ) &&
69  ( daught == EvtPDL::getId( std::string( "K*0" ) ) ||
70  daught == EvtPDL::getId( std::string( "anti-K*0" ) ) ||
71  daught == EvtPDL::getId( std::string( "K*+" ) ) ||
72  daught == EvtPDL::getId( std::string( "K*-" ) ) ) ) {
73  btod = true;
74  }
75 
76  EvtVector4R p4b;
77  p4b.set( parent->mass(), 0.0, 0.0, 0.0 );
78  EvtVector4R p4meson = parent->getDaug( 0 )->getP4();
79 
80  EvtVector4C l11, l12;
81  EvtVector4C l21, l22;
82 
83  EvtVector4C a11, a12;
84  EvtVector4C a21, a22;
85 
86  EvtId parentID = parent->getId();
87 
88  //EvtId l_num = parent->getDaug(1)->getId();
89 
90  EvtVector4R pbhat = p4b / parentmass;
91  EvtVector4R qhat = q / parentmass;
92  EvtVector4R pkstarhat = p4meson / parentmass;
93  EvtVector4R phat = pbhat + pkstarhat;
94 
95  EvtComplex c7eff = EvtbTosllAmp::GetC7Eff( q2, nnlo );
96  EvtComplex c9eff = EvtbTosllAmp::GetC9Eff( q2, nnlo, btod );
97  EvtComplex c10eff = EvtbTosllAmp::GetC10Eff( q2, nnlo );
98  EvtComplex uniti( 0.0, 1.0 );
99 
100  double mhatb = 4.4 / ( parentmass );
101  double mhatkstar = mesonmass / ( parentmass );
102  double shat = q2 / ( parentmass * parentmass );
103 
104  EvtComplex a;
105  a = c9eff * v * 2 / ( 1 + mhatkstar ) + 4 * mhatb * c7eff * t1 / shat;
106  EvtComplex b;
107  b = ( 1 + mhatkstar ) *
108  ( c9eff * a1 + 2 * mhatb * ( 1 - mhatkstar ) * c7eff * t2 / shat );
109  EvtComplex c;
110  c = ( ( 1 - mhatkstar ) * c9eff * a2 +
111  2 * mhatb * c7eff * ( t3 + ( 1 - mhatkstar * mhatkstar ) * t2 / shat ) ) /
112  ( 1 - mhatkstar * mhatkstar );
113  EvtComplex d;
114  d = ( c9eff * ( ( 1 + mhatkstar ) * a1 - ( 1 - mhatkstar ) * a2 -
115  2 * mhatkstar * a0 ) -
116  2 * mhatb * c7eff * t3 ) /
117  shat;
118  EvtComplex e;
119  e = 2 * c10eff * v / ( 1 + mhatkstar );
120  EvtComplex f;
121  f = ( 1 + mhatkstar ) * c10eff * a1;
122  EvtComplex g;
123  g = c10eff * a2 / ( 1 + mhatkstar );
124  EvtComplex h;
125  h = c10eff *
126  ( ( 1 + mhatkstar ) * a1 - ( 1 - mhatkstar ) * a2 - 2 * mhatkstar * a0 ) /
127  shat;
128 
129  EvtTensor4C T1, T2;
130 
131  static EvtIdSet bmesons( "B-", "anti-B0", "anti-B_s0" );
132  static EvtIdSet bbarmesons( "B+", "B0", "B_s0" );
133 
134  EvtParticle* lepPlus = 0;
135  EvtParticle* lepMinus = 0;
136 
137  int charge1 = EvtPDL::chg3( parent->getDaug( 1 )->getId() );
138  int charge2 = EvtPDL::chg3( parent->getDaug( 2 )->getId() );
139 
140  lepPlus = ( charge1 > charge2 ) ? parent->getDaug( 1 ) : parent->getDaug( 2 );
141  lepMinus = ( charge1 < charge2 ) ? parent->getDaug( 1 )
142  : parent->getDaug( 2 );
143 
144  if ( bmesons.contains( parentID ) ) {
145  T1 = a * dual( EvtGenFunctions::directProd( pbhat, pkstarhat ) ) -
146  b * uniti * EvtTensor4C::g() +
147  c * uniti * EvtGenFunctions::directProd( pbhat, phat ) +
148  d * uniti * EvtGenFunctions::directProd( pbhat, qhat );
149 
150  T2 = e * dual( EvtGenFunctions::directProd( pbhat, pkstarhat ) ) -
151  f * uniti * EvtTensor4C::g() +
152  g * uniti * EvtGenFunctions::directProd( pbhat, phat ) +
153  h * uniti * EvtGenFunctions::directProd( pbhat, qhat );
154 
155  l11 = EvtLeptonVCurrent( lepPlus->spParent( 0 ), lepMinus->spParent( 0 ) );
156  l21 = EvtLeptonVCurrent( lepPlus->spParent( 1 ), lepMinus->spParent( 0 ) );
157  l12 = EvtLeptonVCurrent( lepPlus->spParent( 0 ), lepMinus->spParent( 1 ) );
158  l22 = EvtLeptonVCurrent( lepPlus->spParent( 1 ), lepMinus->spParent( 1 ) );
159 
160  a11 = EvtLeptonACurrent( lepPlus->spParent( 0 ), lepMinus->spParent( 0 ) );
161  a21 = EvtLeptonACurrent( lepPlus->spParent( 1 ), lepMinus->spParent( 0 ) );
162  a12 = EvtLeptonACurrent( lepPlus->spParent( 0 ), lepMinus->spParent( 1 ) );
163  a22 = EvtLeptonACurrent( lepPlus->spParent( 1 ), lepMinus->spParent( 1 ) );
164 
165  } else {
166  if ( bbarmesons.contains( parentID ) ) {
167  T1 = -a * dual( EvtGenFunctions::directProd( pbhat, pkstarhat ) ) -
168  b * uniti * EvtTensor4C::g() +
169  c * uniti * EvtGenFunctions::directProd( pbhat, phat ) +
170  d * uniti * EvtGenFunctions::directProd( pbhat, qhat );
171 
172  T2 = -e * dual( EvtGenFunctions::directProd( pbhat, pkstarhat ) ) -
173  f * uniti * EvtTensor4C::g() +
174  g * uniti * EvtGenFunctions::directProd( pbhat, phat ) +
175  h * uniti * EvtGenFunctions::directProd( pbhat, qhat );
176 
177  l11 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
178  lepMinus->spParent( 1 ) );
179  l21 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
180  lepMinus->spParent( 1 ) );
181  l12 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
182  lepMinus->spParent( 0 ) );
183  l22 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
184  lepMinus->spParent( 0 ) );
185 
186  a11 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
187  lepMinus->spParent( 1 ) );
188  a21 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
189  lepMinus->spParent( 1 ) );
190  a12 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
191  lepMinus->spParent( 0 ) );
192  a22 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
193  lepMinus->spParent( 0 ) );
194 
195  } else {
196  EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Wrong lepton number\n";
197  T1.zero();
198  T2.zero(); // Set all tensor terms to zero.
199  }
200  }
201 
202  int i;
203 
204  for ( i = 0; i < 3; i++ ) {
205  EvtVector4C eps = parent->getDaug( 0 )->epsParent( i ).conj();
206 
207  EvtVector4C E1 = T1.cont1( eps );
208  EvtVector4C E2 = T2.cont1( eps );
209 
210  amp.vertex( i, 0, 0, l11 * E1 + a11 * E2 );
211  amp.vertex( i, 0, 1, l12 * E1 + a12 * E2 );
212  amp.vertex( i, 1, 0, l21 * E1 + a21 * E2 );
213  amp.vertex( i, 1, 1, l22 * E1 + a22 * E2 );
214  }
215 }
EvtTensor3C directProd(const EvtVector3C &c1, const EvtVector3C &c2)
EvtComplex GetC7Eff(double q2, bool nnlo=true)
static const EvtTensor4C & g()
Definition: EvtTensor4C.cpp:44
EvtTensor4C dual(const EvtTensor4C &t2)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
const double a2
EvtId getId() const
void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtbTosllFF *formFactors) override
const double a1
void set(int i, double d)
Definition: EvtVector4R.hh:167
static int chg3(EvtId i)
Definition: EvtPDL.cpp:372
virtual EvtDiracSpinor spParent(int) const
EvtComplex GetC10Eff(double q2, bool nnlo=true)
Definition: EvtId.hh:27
EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void vertex(const EvtComplex &amp)
Definition: EvtAmp.cpp:454
virtual void getVectorFF(EvtId, EvtId, double, double, double &, double &, double &, double &, double &, double &, double &)
Definition: EvtbTosllFF.hh:38
double mass2() const
Definition: EvtVector4R.hh:100
virtual EvtVector4C epsParent(int i) const
Definition: EvtAmp.hh:30
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
const EvtVector4R & getP4() const
EvtTensor3C eps(const EvtVector3R &v)
double mass() const
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
EvtVector4C conj() const
Definition: EvtVector4C.hh:202
EvtComplex GetC9Eff(double q2, bool nnlo=true, bool btod=false)
int contains(const EvtId id)
Definition: EvtIdSet.cpp:422