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.
EvtGoityRoberts.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/EvtGoityRoberts.hh"
22 
23 #include "EvtGenBase/EvtDiracSpinor.hh"
24 #include "EvtGenBase/EvtGenKine.hh"
25 #include "EvtGenBase/EvtPDL.hh"
26 #include "EvtGenBase/EvtParticle.hh"
27 #include "EvtGenBase/EvtPatches.hh"
28 #include "EvtGenBase/EvtReport.hh"
29 #include "EvtGenBase/EvtTensor4C.hh"
30 #include "EvtGenBase/EvtVector4C.hh"
31 
32 #include <stdlib.h>
33 #include <string>
34 
35 std::string EvtGoityRoberts::getName()
36 {
37  return "GOITY_ROBERTS";
38 }
39 
40 EvtDecayBase* EvtGoityRoberts::clone()
41 {
42  return new EvtGoityRoberts;
43 }
44 
45 void EvtGoityRoberts::init()
46 {
47  // check that there are 0 arguments
48  checkNArg( 0 );
49  checkNDaug( 4 );
50 
51  checkSpinParent( EvtSpinType::SCALAR );
52  checkSpinDaughter( 1, EvtSpinType::SCALAR );
53  checkSpinDaughter( 2, EvtSpinType::DIRAC );
54  checkSpinDaughter( 3, EvtSpinType::NEUTRINO );
55 }
56 
57 void EvtGoityRoberts::initProbMax()
58 {
59  setProbMax( 3000.0 );
60 }
61 
62 void EvtGoityRoberts::decay( EvtParticle* p )
63 {
64  //added by Lange Jan4,2000
65  static EvtId DST0 = EvtPDL::getId( "D*0" );
66  static EvtId DSTB = EvtPDL::getId( "anti-D*0" );
67  static EvtId DSTP = EvtPDL::getId( "D*+" );
68  static EvtId DSTM = EvtPDL::getId( "D*-" );
69  static EvtId D0 = EvtPDL::getId( "D0" );
70  static EvtId D0B = EvtPDL::getId( "anti-D0" );
71  static EvtId DP = EvtPDL::getId( "D+" );
72  static EvtId DM = EvtPDL::getId( "D-" );
73 
74  EvtId meson = getDaug( 0 );
75 
76  if ( meson == DST0 || meson == DSTP || meson == DSTM || meson == DSTB ) {
77  DecayBDstarpilnuGR( p, getDaug( 0 ), getDaug( 2 ), getDaug( 3 ) );
78  } else {
79  if ( meson == D0 || meson == DP || meson == DM || meson == D0B ) {
80  DecayBDpilnuGR( p, getDaug( 0 ), getDaug( 2 ), getDaug( 3 ) );
81  } else {
82  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
83  << "Wrong daugther in EvtGoityRoberts!\n";
84  }
85  }
86  return;
87 }
88 
89 void EvtGoityRoberts::DecayBDstarpilnuGR( EvtParticle* pb, EvtId ndstar,
90  EvtId nlep, EvtId /*nnu*/ )
91 {
92  pb->initializePhaseSpace( getNDaug(), getDaugs() );
93 
94  //added by Lange Jan4,2000
95  static EvtId EM = EvtPDL::getId( "e-" );
96  static EvtId EP = EvtPDL::getId( "e+" );
97  static EvtId MUM = EvtPDL::getId( "mu-" );
98  static EvtId MUP = EvtPDL::getId( "mu+" );
99 
100  EvtParticle *dstar, *pion, *lepton, *neutrino;
101 
102  // pb->makeDaughters(getNDaug(),getDaugs());
103  dstar = pb->getDaug( 0 );
104  pion = pb->getDaug( 1 );
105  lepton = pb->getDaug( 2 );
106  neutrino = pb->getDaug( 3 );
107 
108  EvtVector4C l1, l2, et0, et1, et2;
109 
110  EvtVector4R v, vp, p4_pi;
111  double w;
112 
113  v.set( 1.0, 0.0, 0.0, 0.0 ); //4-velocity of B meson
114  vp = ( 1.0 / dstar->getP4().mass() ) * dstar->getP4(); //4-velocity of D*
115  p4_pi = pion->getP4();
116 
117  w = v * vp; //four velocity transfere.
118 
119  EvtTensor4C omega;
120 
121  double mb = EvtPDL::getMeanMass( pb->getId() ); //B mass
122  double md = EvtPDL::getMeanMass( ndstar ); //D* mass
123 
124  EvtComplex dmb( 0.0460, -0.5 * 0.00001 ); // B*-B mass splitting ?
125  EvtComplex dmd( 0.1421, -0.5 * 0.00006 );
126  // The last two sets of numbers should
127  // be correctly calculated from the
128  // dstar and pion charges.
129  double g = 0.5; // EvtAmplitude proportional to these coupling constants
130  double alpha3 = 0.690; // See table I in G&R's paper
131  double alpha1 = -1.430;
132  double alpha2 = -0.140;
133  double f0 = 0.093; // The pion decay constants set to 93 MeV
134 
135  EvtComplex dmt3( 0.563, -0.5 * 0.191 ); // Mass splitting = dmt - iGamma/2
136  EvtComplex dmt1( 0.392, -0.5 * 1.040 );
137  EvtComplex dmt2( 0.709, -0.5 * 0.405 );
138 
139  double betas = 0.285; // magic number for meson wave function ground state
140  double betap = 0.280; // magic number for meson wave function state "1"
141  double betad = 0.260; // magic number for meson wave function state "2"
142  double betasp = betas * betas + betap * betap;
143  double betasd = betas * betas + betad * betad;
144 
145  double lambdabar = 0.750; //M(0-,1-) - mQ From Goity&Roberts's code
146 
147  // Isgur&Wise fct
148  double xi = exp( lambdabar * lambdabar * ( 1.0 - w * w ) /
149  ( 4 * betas * betas ) );
150  double xi1 =
151  -1.0 * sqrt( 2.0 / 3.0 ) *
152  ( lambdabar * lambdabar * ( w * w - 1.0 ) / ( 4 * betas * betas ) ) *
153  exp( lambdabar * lambdabar * ( 1.0 - w * w ) / ( 4 * betas * betas ) );
154  double rho1 = sqrt( 1.0 / 2.0 ) * ( lambdabar / betas ) *
155  pow( ( 2 * betas * betap / ( betasp ) ), 2.5 ) *
156  exp( lambdabar * lambdabar * ( 1.0 - w * w ) / ( 2 * betasp ) );
157  double rho2 = sqrt( 1.0 / 8.0 ) *
158  ( lambdabar * lambdabar / ( betas * betas ) ) *
159  pow( ( 2 * betas * betad / ( betasd ) ), 3.5 ) *
160  exp( lambdabar * lambdabar * ( 1.0 - w * w ) / ( 2 * betasd ) );
161 
162  //EvtGenReport(EVTGEN_INFO,"EvtGen") <<"rho's:"<<rho1<<rho2<<endl;
163 
164  EvtComplex h1nr, h2nr, h3nr, f1nr, f2nr;
165  EvtComplex f3nr, f4nr, f5nr, f6nr, knr, g1nr, g2nr, g3nr, g4nr, g5nr;
166  EvtComplex h1r, h2r, h3r, f1r, f2r, f3r, f4r, f5r, f6r, kr, g1r, g2r, g3r,
167  g4r, g5r;
168  EvtComplex h1, h2, h3, f1, f2, f3, f4, f5, f6, k, g1, g2, g3, g4, g5;
169 
170  // Non-resonance part
171  h1nr = -g * xi * ( p4_pi * v ) /
172  ( f0 * mb * md * ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) );
173  h2nr = -g * xi / ( f0 * mb * ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) );
174  h3nr = -( g * xi / ( f0 * md ) ) *
175  ( 1.0 / ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) -
176  EvtComplex( ( 1.0 + w ) / ( p4_pi * vp ), 0.0 ) );
177 
178  f1nr = -( g * xi / ( 2 * f0 * mb ) ) *
179  ( 1.0 / ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) -
180  1.0 / ( EvtComplex( p4_pi * vp, 0.0 ) + dmd ) );
181  f2nr = f1nr * mb / md;
182  f3nr = EvtComplex( 0.0 );
183  f4nr = EvtComplex( 0.0 );
184  f5nr = ( g * xi / ( 2 * f0 * mb * md ) ) *
185  ( EvtComplex( 1.0, 0.0 ) +
186  ( p4_pi * v ) / ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) );
187  f6nr = ( g * xi / ( 2 * f0 * mb ) ) *
188  ( 1.0 / ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) -
189  EvtComplex( 1.0 / ( p4_pi * vp ), 0.0 ) );
190 
191  knr = ( g * xi / ( 2 * f0 ) ) *
192  ( ( p4_pi * ( vp - w * v ) ) / ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) +
193  EvtComplex( ( p4_pi * ( v - w * vp ) ) / ( p4_pi * vp ), 0.0 ) );
194 
195  g1nr = EvtComplex( 0.0 );
196  g2nr = EvtComplex( 0.0 );
197  g3nr = EvtComplex( 0.0 );
198  g4nr = ( g * xi ) / ( f0 * md * EvtComplex( p4_pi * vp ) );
199  g5nr = EvtComplex( 0.0 );
200 
201  // Resonance part (D** removed by hand - alainb)
202  h1r = -alpha1 * rho1 * ( p4_pi * v ) /
203  ( f0 * mb * md * ( EvtComplex( p4_pi * v, 0.0 ) + dmt1 ) ) +
204  alpha2 * rho2 * ( p4_pi * ( v + 2.0 * w * v - vp ) ) /
205  ( 3 * f0 * mb * md * ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) -
206  alpha3 * xi1 * ( p4_pi * v ) /
207  ( f0 * mb * md * EvtComplex( p4_pi * v, 0.0 ) + dmt3 );
208  h2r = -alpha2 * ( 1 + w ) * rho2 /
209  ( 3 * f0 * mb * ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) -
210  alpha3 * xi1 / ( f0 * mb * ( EvtComplex( p4_pi * v, 0.0 ) + dmt3 ) );
211  h3r = alpha2 * rho2 * ( 1 + w ) /
212  ( 3 * f0 * md * ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) -
213  alpha3 * xi1 / ( f0 * md * ( EvtComplex( p4_pi * v, 0.0 ) + dmt3 ) );
214 
215  f1r = -alpha2 * rho2 * ( w - 1.0 ) /
216  ( 6 * f0 * mb * ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) -
217  alpha3 * xi1 /
218  ( 2 * f0 * mb * ( EvtComplex( p4_pi * v, 0.0 ) + dmt3 ) );
219  f2r = f1r * mb / md;
220  f3r = EvtComplex( 0.0 );
221  f4r = EvtComplex( 0.0 );
222  f5r = alpha1 * rho1 * ( p4_pi * v ) /
223  ( 2 * f0 * mb * md * ( EvtComplex( p4_pi * v, 0.0 ) + dmt1 ) ) +
224  alpha2 * rho2 * ( p4_pi * ( vp - v / 3.0 - 2.0 / 3.0 * w * v ) ) /
225  ( 2 * f0 * mb * md * ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) +
226  alpha3 * xi1 * ( p4_pi * v ) /
227  ( 2 * f0 * mb * md * ( EvtComplex( p4_pi * v, 0.0 ) + dmt3 ) );
228  f6r = alpha2 * rho2 * ( w - 1.0 ) /
229  ( 6 * f0 * mb * ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) +
230  alpha3 * xi1 /
231  ( 2 * f0 * mb * ( EvtComplex( p4_pi * v, 0.0 ) + dmt3 ) );
232 
233  kr = -alpha1 * rho1 * ( w - 1.0 ) * ( p4_pi * v ) /
234  ( 2 * f0 * ( EvtComplex( p4_pi * v, 0.0 ) + dmt1 ) ) -
235  alpha2 * rho2 * ( w - 1.0 ) * ( p4_pi * ( vp - w * v ) ) /
236  ( 3 * f0 * ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) +
237  alpha3 * xi1 * ( p4_pi * ( vp - w * v ) ) /
238  ( 2 * f0 * ( EvtComplex( p4_pi * v, 0.0 ) + dmt3 ) );
239 
240  g1r = EvtComplex( 0.0 );
241  g2r = EvtComplex( 0.0 );
242  g3r = -g2r;
243  g4r = 2.0 * alpha2 * rho2 /
244  ( 3 * f0 * md * ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) );
245  g5r = EvtComplex( 0.0 );
246 
247  //Sum
248  h1 = h1nr + h1r;
249  h2 = h2nr + h2r;
250  h3 = h3nr + h3r;
251 
252  f1 = f1nr + f1r;
253  f2 = f2nr + f2r;
254  f3 = f3nr + f3r;
255  f4 = f4nr + f4r;
256  f5 = f5nr + f5r;
257  f6 = f6nr + f6r;
258 
259  k = knr + kr;
260 
261  g1 = g1nr + g1r;
262  g2 = g2nr + g2r;
263  g3 = g3nr + g3r;
264  g4 = g4nr + g4r;
265  g5 = g5nr + g5r;
266 
267  EvtTensor4C g_metric;
268  g_metric.setdiag( 1.0, -1.0, -1.0, -1.0 );
269 
270  if ( nlep == EM || nlep == MUM ) {
271  omega =
272  EvtComplex( 0.0, 0.5 ) *
273  dual( h1 * mb * md * EvtGenFunctions::directProd( v, vp ) +
274  h2 * mb * EvtGenFunctions::directProd( v, p4_pi ) +
275  h3 * md * EvtGenFunctions::directProd( vp, p4_pi ) ) +
276  f1 * mb * EvtGenFunctions::directProd( v, p4_pi ) +
277  f2 * md * EvtGenFunctions::directProd( vp, p4_pi ) +
278  f3 * EvtGenFunctions::directProd( p4_pi, p4_pi ) +
279  f4 * mb * mb * EvtGenFunctions::directProd( v, v ) +
280  f5 * mb * md * EvtGenFunctions::directProd( vp, v ) +
281  f6 * mb * EvtGenFunctions::directProd( p4_pi, v ) + k * g_metric +
282  EvtComplex( 0.0, 0.5 ) *
283  EvtGenFunctions::directProd(
284  dual( EvtGenFunctions::directProd( vp, p4_pi ) ).cont2( v ),
285  ( g1 * p4_pi + g2 * mb * v ) ) +
286  EvtComplex( 0.0, 0.5 ) *
287  EvtGenFunctions::directProd(
288  ( g3 * mb * v + g4 * md * vp + g5 * p4_pi ),
289  dual( EvtGenFunctions::directProd( vp, p4_pi ) ).cont2( v ) );
290 
291  l1 = EvtLeptonVACurrent( lepton->spParent( 0 ),
292  neutrino->spParentNeutrino() );
293  l2 = EvtLeptonVACurrent( lepton->spParent( 1 ),
294  neutrino->spParentNeutrino() );
295  } else {
296  if ( nlep == EP || nlep == MUP ) {
297  omega =
298  EvtComplex( 0.0, -0.5 ) *
299  dual( h1 * mb * md * EvtGenFunctions::directProd( v, vp ) +
300  h2 * mb * EvtGenFunctions::directProd( v, p4_pi ) +
301  h3 * md * EvtGenFunctions::directProd( vp, p4_pi ) ) +
302  f1 * mb * EvtGenFunctions::directProd( v, p4_pi ) +
303  f2 * md * EvtGenFunctions::directProd( vp, p4_pi ) +
304  f3 * EvtGenFunctions::directProd( p4_pi, p4_pi ) +
305  f4 * mb * mb * EvtGenFunctions::directProd( v, v ) +
306  f5 * mb * md * EvtGenFunctions::directProd( vp, v ) +
307  f6 * mb * EvtGenFunctions::directProd( p4_pi, v ) + k * g_metric +
308  EvtComplex( 0.0, -0.5 ) *
309  EvtGenFunctions::directProd(
310  dual( EvtGenFunctions::directProd( vp, p4_pi ) ).cont2( v ),
311  ( g1 * p4_pi + g2 * mb * v ) ) +
312  EvtComplex( 0.0, -0.5 ) *
313  EvtGenFunctions::directProd(
314  ( g3 * mb * v + g4 * md * vp + g5 * p4_pi ),
315  dual( EvtGenFunctions::directProd( vp, p4_pi ) ).cont2( v ) );
316 
317  l1 = EvtLeptonVACurrent( neutrino->spParentNeutrino(),
318  lepton->spParent( 0 ) );
319  l2 = EvtLeptonVACurrent( neutrino->spParentNeutrino(),
320  lepton->spParent( 1 ) );
321  } else {
322  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
323  << "42387dfs878w wrong lepton number\n";
324  }
325  }
326 
327  et0 = omega.cont2( dstar->epsParent( 0 ).conj() );
328  et1 = omega.cont2( dstar->epsParent( 1 ).conj() );
329  et2 = omega.cont2( dstar->epsParent( 2 ).conj() );
330 
331  vertex( 0, 0, l1.cont( et0 ) );
332  vertex( 0, 1, l2.cont( et0 ) );
333 
334  vertex( 1, 0, l1.cont( et1 ) );
335  vertex( 1, 1, l2.cont( et1 ) );
336 
337  vertex( 2, 0, l1.cont( et2 ) );
338  vertex( 2, 1, l2.cont( et2 ) );
339 
340  return;
341 }
342 
343 void EvtGoityRoberts::DecayBDpilnuGR( EvtParticle* pb, EvtId nd, EvtId nlep,
344  EvtId /*nnu*/ )
345 
346 {
347  //added by Lange Jan4,2000
348  static EvtId EM = EvtPDL::getId( "e-" );
349  static EvtId EP = EvtPDL::getId( "e+" );
350  static EvtId MUM = EvtPDL::getId( "mu-" );
351  static EvtId MUP = EvtPDL::getId( "mu+" );
352 
353  EvtParticle *d, *pion, *lepton, *neutrino;
354 
355  pb->initializePhaseSpace( getNDaug(), getDaugs() );
356  d = pb->getDaug( 0 );
357  pion = pb->getDaug( 1 );
358  lepton = pb->getDaug( 2 );
359  neutrino = pb->getDaug( 3 );
360 
361  EvtVector4C l1, l2, et0, et1, et2;
362 
363  EvtVector4R v, vp, p4_pi;
364  double w;
365 
366  v.set( 1.0, 0.0, 0.0, 0.0 ); //4-velocity of B meson
367  vp = ( 1.0 / d->getP4().mass() ) * d->getP4(); //4-velocity of D
368  p4_pi = pion->getP4(); //4-momentum of pion
369  w = v * vp; //four velocity transfer.
370 
371  double mb = EvtPDL::getMeanMass( pb->getId() ); //B mass
372  double md = EvtPDL::getMeanMass( nd ); //D* mass
373  EvtComplex dmb( 0.0460, -0.5 * 0.00001 ); //B mass splitting ?
374  //The last two numbers should be
375  //correctly calculated from the
376  //dstar and pion particle number.
377 
378  double g = 0.5; // Amplitude proportional to these coupling constants
379  double alpha3 = 0.690; // See table I in G&R's paper
380  double alpha1 = -1.430;
381  double alpha2 = -0.140;
382  double f0 = 0.093; // The pion decay constant set to 93 MeV
383 
384  EvtComplex dmt3( 0.563, -0.5 * 0.191 ); // Mass splitting = dmt - iGamma/2
385  EvtComplex dmt1( 0.392, -0.5 * 1.040 );
386  EvtComplex dmt2( 0.709, -0.5 * 0.405 );
387 
388  double betas = 0.285; // magic number for meson wave function ground state
389  double betap = 0.280; // magic number for meson wave function state "1"
390  double betad = 0.260; // magic number for meson wave function state "2"
391  double betasp = betas * betas + betap * betap;
392  double betasd = betas * betas + betad * betad;
393 
394  double lambdabar = 0.750; //M(0-,1-) - mQ From Goity&Roberts's code
395 
396  // Isgur&Wise fct
397  double xi = exp( lambdabar * lambdabar * ( 1.0 - w * w ) /
398  ( 4 * betas * betas ) );
399  double xi1 =
400  -1.0 * sqrt( 2.0 / 3.0 ) *
401  ( lambdabar * lambdabar * ( w * w - 1.0 ) / ( 4 * betas * betas ) ) *
402  exp( lambdabar * lambdabar * ( 1.0 - w * w ) / ( 4 * betas * betas ) );
403  double rho1 = sqrt( 1.0 / 2.0 ) * ( lambdabar / betas ) *
404  pow( ( 2 * betas * betap / ( betasp ) ), 2.5 ) *
405  exp( lambdabar * lambdabar * ( 1.0 - w * w ) / ( 2 * betasp ) );
406  double rho2 = sqrt( 1.0 / 8.0 ) *
407  ( lambdabar * lambdabar / ( betas * betas ) ) *
408  pow( ( 2 * betas * betad / ( betasd ) ), 3.5 ) *
409  exp( lambdabar * lambdabar * ( 1.0 - w * w ) / ( 2 * betasd ) );
410 
411  EvtComplex h, a1, a2, a3;
412  EvtComplex hnr, a1nr, a2nr, a3nr;
413  EvtComplex hr, a1r, a2r, a3r;
414 
415  // Non-resonance part (D* and D** removed by hand - alainb)
416  hnr = g * xi * ( 1.0 / ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) ) /
417  ( 2 * f0 * mb * md );
418  a1nr = -1.0 * g * xi * ( 1 + w ) *
419  ( 1.0 / ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) ) / ( 2 * f0 );
420  a2nr = g * xi *
421  ( ( p4_pi * ( v + vp ) ) / ( EvtComplex( p4_pi * v, 0.0 ) + dmb ) ) /
422  ( 2 * f0 * mb );
423  a3nr = EvtComplex( 0.0, 0.0 );
424 
425  // Resonance part (D** remove by hand - alainb)
426  hr = alpha2 * rho2 * ( w - 1 ) *
427  ( 1.0 / ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) /
428  ( 6 * f0 * mb * md ) +
429  alpha3 * xi1 * ( 1.0 / ( EvtComplex( p4_pi * v, 0.0 ) + dmt3 ) ) /
430  ( 2 * f0 * mb * md );
431  a1r = -1.0 * alpha2 * rho2 * ( w * w - 1 ) *
432  ( 1.0 / ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) / ( 6 * f0 ) -
433  alpha3 * xi1 * ( 1 + w ) *
434  ( 1.0 / ( EvtComplex( p4_pi * v, 0.0 ) + dmt3 ) ) / ( 2 * f0 );
435  a2r = alpha1 * rho1 *
436  ( ( p4_pi * v ) / ( EvtComplex( p4_pi * v, 0.0 ) + dmt1 ) ) /
437  ( 2 * f0 * mb ) +
438  alpha2 * rho2 *
439  ( 0.5 * p4_pi * ( w * vp - v ) + p4_pi * ( vp - w * v ) ) /
440  ( 3 * f0 * mb * ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) +
441  alpha3 * xi1 *
442  ( ( p4_pi * ( v + vp ) ) / ( EvtComplex( p4_pi * v, 0.0 ) + dmt3 ) ) /
443  ( 2 * f0 * mb );
444  a3r = -1.0 * alpha1 * rho1 *
445  ( ( p4_pi * v ) / ( EvtComplex( p4_pi * v, 0.0 ) + dmt1 ) ) /
446  ( 2 * f0 * md ) -
447  alpha2 * rho2 *
448  ( ( p4_pi * ( vp - w * v ) ) /
449  ( EvtComplex( p4_pi * v, 0.0 ) + dmt2 ) ) /
450  ( 2 * f0 * md );
451 
452  // Sum
453  h = hnr + hr;
454  a1 = a1nr + a1r;
455  a2 = a2nr + a2r;
456  a3 = a3nr + a3r;
457 
458  EvtVector4C omega;
459 
460  if ( nlep == EM || nlep == MUM ) {
461  omega = EvtComplex( 0.0, -1.0 ) * h * mb * md *
462  dual( EvtGenFunctions::directProd( vp, p4_pi ) ).cont2( v ) +
463  a1 * p4_pi + a2 * mb * v + a3 * md * vp;
464  l1 = EvtLeptonVACurrent( lepton->spParent( 0 ),
465  neutrino->spParentNeutrino() );
466  l2 = EvtLeptonVACurrent( lepton->spParent( 1 ),
467  neutrino->spParentNeutrino() );
468  } else {
469  if ( nlep == EP || nlep == MUP ) {
470  omega = EvtComplex( 0.0, 1.0 ) * h * mb * md *
471  dual( EvtGenFunctions::directProd( vp, p4_pi ) ).cont2( v ) +
472  a1 * p4_pi + a2 * mb * v + a3 * md * vp;
473  l1 = EvtLeptonVACurrent( neutrino->spParentNeutrino(),
474  lepton->spParent( 0 ) );
475  l2 = EvtLeptonVACurrent( neutrino->spParentNeutrino(),
476  lepton->spParent( 1 ) );
477  } else {
478  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
479  << "42387dfs878w wrong lepton number\n";
480  }
481  }
482 
483  vertex( 0, l1 * omega );
484  vertex( 1, l2 * omega );
485 
486  return;
487 }
EvtGoityRoberts::DecayBDstarpilnuGR
void DecayBDstarpilnuGR(EvtParticle *pb, EvtId ndstar, EvtId nlep, EvtId nnu)
Definition: EvtGoityRoberts.cpp:89
EvtGenFunctions::directProd
EvtTensor3C directProd(const EvtVector3C &c1, const EvtVector3C &c2)
Definition: EvtTensor3C.cpp:154
EvtTensor4C::setdiag
void setdiag(double t00, double t11, double t22, double t33)
Definition: EvtTensor4C.cpp:203
dual
EvtTensor4C dual(const EvtTensor4C &t2)
Definition: EvtTensor4C.cpp:364
a3
const double a3
Definition: EvtLi2Spence.cpp:25
EvtGenReport
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
EvtGoityRoberts::decay
void decay(EvtParticle *p) override
Definition: EvtGoityRoberts.cpp:62
EvtVector4R::mass
double mass() const
Definition: EvtVector4R.cpp:49
a2
const double a2
Definition: EvtLi2Spence.cpp:24
EvtTensor4C::cont2
EvtVector4C cont2(const EvtVector4C &v4) const
Definition: EvtTensor4C.cpp:461
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition: EvtPDL.cpp:314
EvtGoityRoberts::DecayBDpilnuGR
void DecayBDpilnuGR(EvtParticle *pb, EvtId nd, EvtId nlep, EvtId nnu)
Definition: EvtGoityRoberts.cpp:343
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cpp:130
a1
const double a1
Definition: EvtLi2Spence.cpp:23
EvtVector4R::set
void set(int i, double d)
Definition: EvtVector4R.hh:167
EvtParticle::spParent
virtual EvtDiracSpinor spParent(int) const
Definition: EvtParticle.cpp:665
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:66
EvtGoityRoberts::getName
std::string getName() override
Definition: EvtGoityRoberts.cpp:35
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition: EvtDecayBase.cpp:321
EvtId
Definition: EvtId.hh:27
EvtDecayAmp::vertex
void vertex(const EvtComplex &amp)
Definition: EvtDecayAmp.hh:37
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::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition: EvtDecayBase.cpp:542
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
EvtGoityRoberts::initProbMax
void initProbMax() override
Definition: EvtGoityRoberts.cpp:57
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition: EvtDecayBase.cpp:573
EvtParticle::spParentNeutrino
virtual EvtDiracSpinor spParentNeutrino() const
Definition: EvtParticle.cpp:689
EvtDecayBase::getNDaug
int getNDaug() const
Definition: EvtDecayBase.hh:65
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:240
EvtLeptonVACurrent
EvtVector4C EvtLeptonVACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:179
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
EvtVector4C::conj
EvtVector4C conj() const
Definition: EvtVector4C.hh:202
EvtVector4C::cont
EvtComplex cont(const EvtVector4C &v4) const
Definition: EvtVector4C.hh:140
EvtGoityRoberts::init
void init() override
Definition: EvtGoityRoberts.cpp:45
EvtGoityRoberts::clone
EvtDecayBase * clone() override
Definition: EvtGoityRoberts.cpp:40
EvtDecayBase::getDaug
EvtId getDaug(int i) const
Definition: EvtDecayBase.hh:67
Generated by   doxygen 1.8.15