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.
EvtBcVHad.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 
21 #include "EvtGenModels/EvtBcVHad.hh"
22 
23 #include "EvtGenBase/EvtIdSet.hh"
24 #include "EvtGenBase/EvtPDL.hh"
25 #include "EvtGenBase/EvtParticle.hh"
26 #include "EvtGenBase/EvtReport.hh"
27 #include "EvtGenBase/EvtSpinType.hh"
28 #include "EvtGenBase/EvtTensor4C.hh"
29 #include "EvtGenBase/EvtVector4R.hh"
30 
31 #include "EvtGenModels/EvtBCVFF2.hh"
32 #include "EvtGenModels/EvtWHad.hh"
33 
34 #include <iostream>
35 
36 std::string EvtBcVHad::getName()
37 {
38  return "BC_VHAD";
39 }
40 
41 EvtDecayBase* EvtBcVHad::clone()
42 {
43  return new EvtBcVHad;
44 }
45 
46 //======================================================
47 
48 void EvtBcVHad::init()
49 {
50  checkNArg( 1 );
51  checkSpinParent( EvtSpinType::SCALAR );
52  checkSpinDaughter( 0, EvtSpinType::VECTOR );
53  for ( int i = 1; i <= ( getNDaug() - 1 ); i++ ) {
54  checkSpinDaughter( i, EvtSpinType::SCALAR );
55  }
56 
57  if ( getNDaug() < 2 || getNDaug() > 6 ) {
58  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
59  << "Have not yet implemented this final state in BC_VHAD model"
60  << std::endl;
61  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
62  << "Ndaug=" << getNDaug() << std::endl;
63  for ( int id = 0; id < ( getNDaug() - 1 ); id++ ) {
64  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
65  << "Daug " << id << " " << EvtPDL::name( getDaug( id ) ).c_str()
66  << std::endl;
67  }
68  return;
69  }
70 
71  idVector = getDaug( 0 ).getId();
72  whichfit = int( getArg( 0 ) + 0.1 );
73  ffmodel = std::make_unique<EvtBCVFF2>( idVector, whichfit );
74 
75  wcurr = std::make_unique<EvtWHad>();
76 
77  // determine the code of final hadronic state
78  EvtIdSet thePis( "pi+", "pi-", "pi0" );
79  EvtIdSet theK( "K+", "K-", "K_S0" );
80  if ( getNDaug() == 2 && thePis.contains( getDaug( 1 ) ) ) {
81  out_code = 1;
82  } else if ( getNDaug() == 3 && thePis.contains( getDaug( 1 ) ) &&
83  thePis.contains( getDaug( 2 ) ) ) {
84  out_code = 2;
85  } else if ( getNDaug() == 4 && thePis.contains( getDaug( 1 ) ) &&
86  thePis.contains( getDaug( 2 ) ) &&
87  thePis.contains( getDaug( 3 ) ) ) {
88  out_code = 3;
89  } else if ( getNDaug() == 5 && thePis.contains( getDaug( 1 ) ) &&
90  thePis.contains( getDaug( 2 ) ) &&
91  thePis.contains( getDaug( 3 ) ) &&
92  thePis.contains( getDaug( 4 ) ) ) {
93  out_code = 4;
94  } else if ( getNDaug() == 6 && thePis.contains( getDaug( 1 ) ) &&
95  thePis.contains( getDaug( 2 ) ) &&
96  thePis.contains( getDaug( 3 ) ) &&
97  thePis.contains( getDaug( 4 ) ) &&
98  thePis.contains( getDaug( 5 ) ) ) {
99  out_code = 5;
100  } else if ( getNDaug() == 4 && theK.contains( getDaug( 1 ) ) &&
101  theK.contains( getDaug( 2 ) ) &&
102  thePis.contains( getDaug( 3 ) ) ) {
103  out_code = 6;
104  } else if ( getNDaug() == 4 && theK.contains( getDaug( 1 ) ) &&
105  thePis.contains( getDaug( 2 ) ) &&
106  thePis.contains( getDaug( 3 ) ) ) {
107  out_code = 7;
108  } else if ( getNDaug() == 3 && theK.contains( getDaug( 1 ) ) &&
109  theK.contains( getDaug( 2 ) ) ) {
110  out_code = 8;
111  }
112 }
113 
114 //======================================================
115 
116 void EvtBcVHad::initProbMax()
117 {
118  if ( out_code == 3 ) {
119  if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
120  getNDaug() == 4 )
121  setProbMax( 42000. );
122  else if ( idVector == EvtPDL::getId( "J/psi" ).getId() &&
123  whichfit == 2 && getNDaug() == 4 )
124  setProbMax( 16000. );
125  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
126  whichfit == 1 && getNDaug() == 4 )
127  setProbMax( 1200. );
128  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
129  whichfit == 2 && getNDaug() == 4 )
130  setProbMax( 2600. );
131 
132  } else if ( out_code == 5 ) {
133  if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
134  getNDaug() == 6 )
135  setProbMax( 720000. );
136  else if ( idVector == EvtPDL::getId( "J/psi" ).getId() &&
137  whichfit == 2 && getNDaug() == 6 )
138  setProbMax( 471817. );
139  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
140  whichfit == 1 && getNDaug() == 6 )
141  setProbMax( 40000. );
142  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
143  whichfit == 2 && getNDaug() == 6 )
144  setProbMax( 30000. );
145 
146  } else if ( out_code == 6 ) {
147  if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 )
148  setProbMax( 50000. );
149  else if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 2 )
150  setProbMax( 20000.0 );
151  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 1 )
152  setProbMax( 2100.0 );
153  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 2 )
154  setProbMax( 1700.00 );
155 
156  } else if ( out_code == 7 ) {
157  if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 )
158  setProbMax( 2.2e+06 );
159  else if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 2 )
160  setProbMax( 830000 );
161  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 1 )
162  setProbMax( 92000.0 );
163  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 2 )
164  setProbMax( 93000.0 );
165 
166  } else if ( out_code == 8 ) {
167  if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 )
168  setProbMax( 2e2 );
169  else if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 2 )
170  setProbMax( 80 );
171  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 1 )
172  setProbMax( 10 );
173  else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() && whichfit == 2 )
174  setProbMax( 10 );
175 
176  } else {
177  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
178  << "Have not yet implemented this final state in BC_VHAD model"
179  << std::endl;
180  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
181  << "Ndaug=" << getNDaug() << std::endl;
182  for ( int id = 0; id < ( getNDaug() - 1 ); id++ ) {
183  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
184  << "Daug " << id << " " << EvtPDL::name( getDaug( id ) ).c_str()
185  << std::endl;
186  }
187  ::abort();
188  }
189 }
190 
191 //======================================================
192 
193 EvtVector4C EvtBcVHad::hardCurr( EvtParticle* root_particle ) const
194 {
195  EvtVector4C hardCur;
196 
197  if ( out_code == 1 ) {
198  // pi+
199  hardCur = wcurr->WCurrent( root_particle->getDaug( 1 )->getP4() );
200 
201  } else if ( out_code == 2 ) {
202  // pi+ pi0
203  hardCur = wcurr->WCurrent( root_particle->getDaug( 1 )->getP4(),
204  root_particle->getDaug( 2 )->getP4() );
205 
206  } else if ( out_code == 3 ) {
207  // pi+ pi+ pi-
208  hardCur = wcurr->WCurrent( root_particle->getDaug( 1 )->getP4(),
209  root_particle->getDaug( 2 )->getP4(),
210  root_particle->getDaug( 3 )->getP4() );
211 
212  } else if ( out_code == 5 ) {
213  // Bc -> psi pi+ pi+ pi- pi- pi+ from Kuhn, Was, hep-ph/0602162
214  hardCur = wcurr->WCurrent( root_particle->getDaug( 1 )->getP4(),
215  root_particle->getDaug( 2 )->getP4(),
216  root_particle->getDaug( 3 )->getP4(),
217  root_particle->getDaug( 4 )->getP4(),
218  root_particle->getDaug( 5 )->getP4() );
219 
220  } else if ( out_code == 6 ) {
221  // K+ K- pi+
222  hardCur = wcurr->WCurrent_KKP( root_particle->getDaug( 1 )->getP4(),
223  root_particle->getDaug( 2 )->getP4(),
224  root_particle->getDaug( 3 )->getP4() );
225 
226  } else if ( out_code == 7 ) {
227  // K+ pi+ pi-
228  hardCur = wcurr->WCurrent_KPP( root_particle->getDaug( 1 )->getP4(),
229  root_particle->getDaug( 2 )->getP4(),
230  root_particle->getDaug( 3 )->getP4() );
231 
232  } else if ( out_code == 8 ) {
233  // K_S0 K+
234  hardCur = wcurr->WCurrent_KSK( root_particle->getDaug( 1 )->getP4(),
235  root_particle->getDaug( 2 )->getP4() );
236 
237  } else {
238  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
239  << "Have not yet implemented this final state in BC_VHAD model"
240  << std::endl;
241  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
242  << "Ndaug=" << getNDaug() << std::endl;
243  for ( int id = 0; id < ( getNDaug() - 1 ); id++ ) {
244  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
245  << "Daug " << id << " " << EvtPDL::name( getDaug( id ) ).c_str()
246  << std::endl;
247  }
248  ::abort();
249  }
250 
251  return hardCur;
252 }
253 
254 //======================================================
255 
256 void EvtBcVHad::decay( EvtParticle* root_particle )
257 {
258  root_particle->initializePhaseSpace( getNDaug(), getDaugs() );
259 
260  // Calculate hadronic current
261  EvtVector4C hardCur = hardCurr( root_particle );
262 
263  EvtParticle* Jpsi = root_particle->getDaug( 0 );
264 
265  EvtVector4R p4b( root_particle->mass(), 0., 0., 0. ), // Bc momentum
266  p4meson = Jpsi->getP4(), // J/psi momenta
267  Q = p4b - p4meson, p4Sum = p4meson + p4b;
268  double Q2 = Q.mass2();
269 
270  // Calculate Bc -> V W form-factors
271  double a1f( 0.0 ), a2f( 0.0 ), vf( 0.0 ), a0f( 0.0 );
272 
273  double m_meson = Jpsi->mass();
274  double m_b = root_particle->mass();
275  double mVar = m_b + m_meson;
276 
277  ffmodel->getvectorff( root_particle->getId(), Jpsi->getId(), Q2, m_meson,
278  &a1f, &a2f, &vf, &a0f );
279 
280  double a3f = ( mVar / ( 2.0 * m_meson ) ) * a1f -
281  ( ( m_b - m_meson ) / ( 2.0 * m_meson ) ) * a2f;
282 
283  // Calculate Bc -> V W current
284  EvtTensor4C H = a1f * mVar * EvtTensor4C::g();
285  H.addDirProd( ( -a2f / mVar ) * p4b, p4Sum );
286  H += EvtComplex( 0.0, vf / mVar ) *
287  dual( EvtGenFunctions::directProd( p4Sum, Q ) );
288  H.addDirProd( ( a0f - a3f ) * 2.0 * ( m_meson / Q2 ) * p4b, Q );
289  EvtVector4C Heps = H.cont2( hardCur );
290 
291  for ( int i = 0; i < 4; i++ ) {
292  EvtVector4C eps =
293  Jpsi->epsParent( i ).conj(); // psi-meson polarization vector
294  EvtComplex amp = eps * Heps;
295  vertex( i, amp );
296  }
297 }
EvtGenFunctions::directProd
EvtTensor3C directProd(const EvtVector3C &c1, const EvtVector3C &c2)
Definition: EvtTensor3C.cpp:154
EvtBcVHad::init
void init() override
Definition: EvtBcVHad.cpp:48
EvtBcVHad::idVector
int idVector
Definition: EvtBcVHad.hh:56
EvtPDL::name
static std::string name(EvtId i)
Definition: EvtPDL.cpp:382
EvtBcVHad::initProbMax
void initProbMax() override
Definition: EvtBcVHad.cpp:116
EvtTensor4C::g
static const EvtTensor4C & g()
Definition: EvtTensor4C.cpp:44
EvtBcVHad::wcurr
std::unique_ptr< EvtWHad > wcurr
Definition: EvtBcVHad.hh:70
EvtDecayBase::getArg
double getArg(unsigned int j)
Definition: EvtDecayBase.cpp:604
EvtBcVHad::out_code
int out_code
Definition: EvtBcVHad.hh:67
dual
EvtTensor4C dual(const EvtTensor4C &t2)
Definition: EvtTensor4C.cpp:364
EvtGenReport
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
EvtBcVHad::whichfit
int whichfit
Definition: EvtBcVHad.hh:53
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cpp:130
EvtBcVHad::ffmodel
std::unique_ptr< EvtBCVFF2 > ffmodel
Definition: EvtBcVHad.hh:69
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:66
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition: EvtDecayBase.cpp:321
EvtVector4R::mass2
double mass2() const
Definition: EvtVector4R.hh:100
EvtDecayAmp::vertex
void vertex(const EvtComplex &amp)
Definition: EvtDecayAmp.hh:37
EvtBcVHad::clone
EvtDecayBase * clone() override
Definition: EvtBcVHad.cpp:41
EvtParticle::initializePhaseSpace
double initializePhaseSpace(unsigned int numdaughter, EvtId *daughters, bool forceResetMasses=false, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cpp:1102
EvtParticle::epsParent
virtual EvtVector4C epsParent(int i) const
Definition: EvtParticle.cpp:616
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition: EvtDecayBase.cpp:560
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition: EvtDecayBase.cpp:518
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cpp:150
EvtGenFunctions::eps
EvtTensor3C eps(const EvtVector3R &v)
Definition: EvtTensor3C.cpp:303
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition: EvtDecayBase.cpp:573
EvtParticle::mass
double mass() const
Definition: EvtParticle.cpp:165
EvtId::getId
int getId() const
Definition: EvtId.hh:42
EvtDecayBase::getNDaug
int getNDaug() const
Definition: EvtDecayBase.hh:65
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
EvtVector4C::conj
EvtVector4C conj() const
Definition: EvtVector4C.hh:202
EvtIdSet::contains
int contains(const EvtId id)
Definition: EvtIdSet.cpp:422
EvtBcVHad::decay
void decay(EvtParticle *p) override
Definition: EvtBcVHad.cpp:256
EvtBcVHad::getName
std::string getName() override
Definition: EvtBcVHad.cpp:36
EvtBcVHad::hardCurr
EvtVector4C hardCurr(EvtParticle *root_particle) const
Definition: EvtBcVHad.cpp:193
EvtDecayBase::getDaug
EvtId getDaug(int i) const
Definition: EvtDecayBase.hh:67
Generated by   doxygen 1.8.15