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