Belle II Software development
Public Member Functions | Private Member Functions | Private Attributes | List of all members
EvtBToDstarlnuNP Class Reference

Implementation of $B \to D^* \ell \bar{\nu}_\ell$ with new physics contributions and the BGL hadronic form factor parameterization. More...

#include <EvtBToDstarlnuNP.h>

Inheritance diagram for EvtBToDstarlnuNP:

Public Member Functions

std::string getName () override
 The function which returns the name of the model.
 
EvtDecayBase * clone () override
 The function which makes a copy of the model.
 
void decay (EvtParticle *p) override
 The method to calculate a decay amplitude.
 
void init () override
 Initialization method.
 
void initProbMax () override
 The method to evaluate the maximum amplitude.
 

Private Member Functions

double CalcMaxProb ()
 The method to evaluate the maximum decay amplitude.
 
void CalcAmp (EvtParticle *, EvtAmp &)
 The method to evaluate the decay amplitude.
 

Private Attributes

EvtComplex _Cvl
 Strength of the left handed part of the vector hadronic current.
 
EvtComplex _Cvr
 Strength of the Right handed part of the vector hadronic current.
 
EvtComplex _Csl
 (g_S - g_P) strength of pseudoscalar current
 
EvtComplex _Csr
 (g_S + g_P) strength of pseudoscalar current with constraint _Csl = -_Csr
 
EvtComplex _cT
 C_T – tensor current.
 

Detailed Description

Implementation of $B \to D^* \ell \bar{\nu}_\ell$ with new physics contributions and the BGL hadronic form factor parameterization.

Based on the publication: Phys.Rev.D 107 (2023) 1, 015011; e-Print:2206.11283 [hep-ph]

Definition at line 22 of file EvtBToDstarlnuNP.h.

Member Function Documentation

◆ CalcAmp()

void CalcAmp ( EvtParticle * parent,
EvtAmp & amp )
private

The method to evaluate the decay amplitude.

Definition at line 273 of file EvtBToDstarlnuNP.cc.

274{
275 if (cafirst) {
276 cafirst = false;
277 EM = EvtPDL::getId("e-"); EP = EvtPDL::getId("e+");
278 MUM = EvtPDL::getId("mu-"); MUP = EvtPDL::getId("mu+");
279 TAUM = EvtPDL::getId("tau-"); TAUP = EvtPDL::getId("tau+");
280
281 D0 = EvtPDL::getId("D0"); D0B = EvtPDL::getId("anti-D0");
282 DP = EvtPDL::getId("D+"); DM = EvtPDL::getId("D-");
283 DSP = EvtPDL::getId("D_s+"); DSM = EvtPDL::getId("D_s-");
284 }
285
286 EvtParticle
287 *pvm = parent->getDaug(0), // vector meson
288 *plp = parent->getDaug(1), // charged lepton
289 *pnu = parent->getDaug(2); // neutrino
290 EvtId pId = parent->getId(), lId = plp->getId();
291
292 //Add the lepton and neutrino 4 momenta to find q2
293 EvtVector4R q = plp->getP4() + pnu->getP4();
294 double q2 = q.mass2(), m_V = pvm->mass(), m_B = parent->mass();
295 double m_BaV = m_B + m_V, m_BsV = m_B - m_V;
296
297 double a1f, a2f, vf, a0f, a3f;
298 /* Tensor form factors */
299 double T1, T2, T3;
300 getFF(m_B, m_V, q2, vf, a0f, a1f, a2f, T1, T2, T3);
301 a3f = (m_BaV * a1f - m_BsV * a2f) / (2.0 * m_V);
302
303 double c1 = (m_BaV * m_BsV * (T1 - T2)) / q2;
304 double c2 = (T1 - T2 - T3 * q2 / (m_BaV * m_BsV)) * 2 / q2;
305
306 if (pId == D0 || pId == D0B || pId == DP || pId == DM || pId == DSP || pId == DSM)
307 vf = -vf;
308
309 EvtVector4R p(m_B, 0., 0., 0.);
310 const EvtVector4R& k = pvm->getP4();
311 EvtVector4R pak = p + k;
312
313 EvtDiracSpinor
314 nu = pnu->spParentNeutrino(), lp0 = plp->spParent(0), lp1 = plp->spParent(1);
315
316 /* Various lepton currents */
317 EvtComplex p1, p2;
318 EvtVector4C l1, l2, a1, a2;
319 EvtTensor4C z1, z2, w1, w2;
320 if (lId == EM || lId == MUM || lId == TAUM) {
321 l1 = EvtLeptonVCurrent(nu, lp0);
322 l2 = EvtLeptonVCurrent(nu, lp1);
323 a1 = EvtLeptonACurrent(nu, lp0);
324 a2 = EvtLeptonACurrent(nu, lp1);
325 p1 = EvtLeptonPCurrent(nu, lp0);
326 p2 = EvtLeptonPCurrent(nu, lp1);
327 z1 = EvtLeptonTCurrent(nu, lp0);
328 z2 = EvtLeptonTCurrent(nu, lp1);
329 } else if (lId == EP || lId == MUP || lId == TAUP) {
330 vf = -vf; /*This takes care of sign flip in vector hadronic current when you conjugate it*/
331 l1 = EvtLeptonVCurrent(lp0, nu);
332 l2 = EvtLeptonVCurrent(lp1, nu);
333 a1 = EvtLeptonACurrent(lp0, nu);
334 a2 = EvtLeptonACurrent(lp1, nu);
335 p1 = EvtLeptonPCurrent(lp0, nu);
336 p2 = EvtLeptonPCurrent(lp1, nu);
337 z1 = EvtLeptonTCurrent(lp0, nu);
338 z2 = EvtLeptonTCurrent(lp1, nu);
339 }
340
341 /* Quark masses from PDG */
342 double m_c = 1.27, m_b = 4.18;
343
344 /* New Physics */
345 EvtComplex C_V_L = 1, //Standard Model
346 C_V_R = _Cvr,
347 C_S_L = _Csl,
348 C_S_R = _Csr,
349 C_T = _cT;
350
351 C_V_L += _Cvl;
352
353 /*Eq 12 Axial vector hadronic tensor*/
354 EvtTensor4C tA = (a1f * m_BaV) * EvtTensor4C::g();
355 tA.addDirProd(-(a2f / m_BaV) * p, pak);
356 tA.addDirProd(-2 * m_V / q2 * (a3f - a0f) * p, q);
357
358 /*Eq 11 Vector hadronic tensor*/
359 EvtTensor4C tV = EvtComplex(0.0, (-2 * vf / m_BaV)) * asymProd(p, k);
360
361 EvtTensor4C tVA = 0.5 * (C_V_L + C_V_R) * tV + 0.5 * (C_V_R - C_V_L) * tA;
362
363 EvtComplex tS = 0.5 * (C_S_R - C_S_L) * 2 * m_V / (m_b + m_c) * a0f;
364
365 pak *= T1;
366 EvtVector4C qc1(q); qc1 *= c1;
367 EvtTensor4C tt, tt5;
368 EvtComplex sum1, sum2, TT1, TT2;
369 for (int i = 0; i < 3; i++) {
370 EvtVector4C eps = pvm->epsParent(i).conj(), epsTVA = tVA.cont1(eps);
371 EvtComplex epsq = eps * q, epsqc2 = epsq * c2;
372
373 /*Eq 14 Tensor component */
374 for (int i2 = 0; i2 < 4; i2++) {
375 for (int i3 = 0; i3 < 4; i3++) {
376 tt.set(i2, i3, eps.get(i2) * (-pak.get(i3) + qc1.get(i3)) + p.get(i2)*k.get(i3)*epsqc2);
377 }
378 }
379 if (lId == EM || lId == MUM || lId == TAUM) {
380 TT1 = EvtComplex(0.0, 1.0) * C_T * (asymProd(z1, tt));
381 TT2 = EvtComplex(0.0, 1.0) * C_T * (asymProd(z2, tt));
382 sum1 = l1.cont(epsTVA) - a1.cont(epsTVA) - epsq * tS * p1 + TT1;
383 sum2 = l2.cont(epsTVA) - a2.cont(epsTVA) - epsq * tS * p2 + TT2;
384 } else if (lId == EP || lId == MUP || lId == TAUP) {
385 TT1 = EvtComplex(0.0, -1.0) * conj(C_T) * (asymProd(z1, tt));
386 TT2 = EvtComplex(0.0, -1.0) * conj(C_T) * (asymProd(z2, tt));
387 sum1 = l1.cont(epsTVA) - a1.cont(epsTVA) - epsq * tS * p1 + TT1;
388 sum2 = l2.cont(epsTVA) - a2.cont(epsTVA) - epsq * tS * p2 + TT2;
389 }
390
391 amp.vertex(i, 0, sum1);
392 amp.vertex(i, 1, sum2);
393 }
394}
EvtBToDstarlnuNP::_cT
EvtComplex _cT
C_T – tensor current.
Definition EvtBToDstarlnuNP.h:63
EvtBToDstarlnuNP::_Cvl
EvtComplex _Cvl
Strength of the left handed part of the vector hadronic current.
Definition EvtBToDstarlnuNP.h:51
EvtBToDstarlnuNP::_Cvr
EvtComplex _Cvr
Strength of the Right handed part of the vector hadronic current.
Definition EvtBToDstarlnuNP.h:54
EvtBToDstarlnuNP::_Csr
EvtComplex _Csr
(g_S + g_P) strength of pseudoscalar current with constraint _Csl = -_Csr
Definition EvtBToDstarlnuNP.h:60
EvtBToDstarlnuNP::_Csl
EvtComplex _Csl
(g_S - g_P) strength of pseudoscalar current
Definition EvtBToDstarlnuNP.h:57

◆ CalcMaxProb()

double CalcMaxProb ( )
private

The method to evaluate the maximum decay amplitude.

Definition at line 396 of file EvtBToDstarlnuNP.cc.

397{
398 //This routine takes the arguements parent, meson, and lepton
399 //number, and a form factor model, and returns a maximum
400 //probability for this semileptonic form factor model. A
401 //brute force method is used. The 2D cos theta lepton and
402 //q2 phase space is probed.
403
404 //Start by declaring a particle at rest.
405
406 //It only makes sense to have a scalar parent. For now.
407 //This should be generalized later.
408
409 EvtId pId = getParentId(), mId = getDaug(0), l1Id = getDaug(1), l2Id = getDaug(2);
410
411 EvtScalarParticle* scalar_part;
412 EvtParticle* root_part;
413
414 scalar_part = new EvtScalarParticle;
415
416 //cludge to avoid generating random numbers!
417 scalar_part->noLifeTime();
418
419 EvtVector4R p_init;
420
421 p_init.set(EvtPDL::getMass(pId), 0.0, 0.0, 0.0);
422 scalar_part->init(pId, p_init);
423 root_part = (EvtParticle*)scalar_part;
424 root_part->setDiagonalSpinDensity();
425
426 EvtId listdaug[3];
427 listdaug[0] = mId;
428 listdaug[1] = l1Id;
429 listdaug[2] = l2Id;
430
431 EvtAmp amp;
432 amp.init(pId, 3, listdaug);
433
434 root_part->makeDaughters(3, listdaug);
435 EvtParticle* daughter, *lep, *trino;
436 daughter = root_part->getDaug(0);
437 lep = root_part->getDaug(1);
438 trino = root_part->getDaug(2);
439
440 //cludge to avoid generating random numbers!
441 daughter->noLifeTime();
442 lep->noLifeTime();
443 trino->noLifeTime();
444
445 //Initial particle is unpolarized, well it is a scalar so it is
446 //trivial
447 EvtSpinDensity rho;
448 rho.setDiag(root_part->getSpinStates());
449
450 EvtVector4R p4meson, p4lepton, p4nu, p4w;
451 double mass[3], m = root_part->mass(), maxfoundprob = 0.0;
452
453 for (int massiter = 0; massiter < 3; massiter++) {
454 mass[0] = EvtPDL::getMeanMass(mId);
455 mass[1] = EvtPDL::getMeanMass(l1Id);
456 mass[2] = EvtPDL::getMeanMass(l2Id);
457 if (massiter == 1) {
458 mass[0] = EvtPDL::getMinMass(mId);
459 }
460 if (massiter == 2) {
461 mass[0] = EvtPDL::getMaxMass(mId);
462 if ((mass[0] + mass[1] + mass[2]) > m)
463 mass[0] = m - mass[1] - mass[2] - 0.00001;
464 }
465
466 double q2min = mass[1] * mass[1], q2max = (m - mass[0]) * (m - mass[0]);
467
468 for (int imax = 40, i = 0; i < imax; i++) {
469 double q2 = q2min + ((i + 0.5) * (q2max - q2min)) / imax;
470 double erho = (m * m + mass[0] * mass[0] - q2) / (2.0 * m);
471 double prho = sqrt(erho * erho - mass[0] * mass[0]);
472 p4meson.set(erho, 0.0, 0.0, -1.0 * prho);
473 p4w.set(m - erho, 0.0, 0.0, prho);
474
475 //This is in the W rest frame
476 double elepton = (q2 + mass[1] * mass[1]) / (2.0 * sqrt(q2));
477 double plepton = sqrt(elepton * elepton - mass[1] * mass[1]);
478
479 double probctl[3];
480 for (int j = 0; j < 3; j++) {
481 double costl = 0.99 * (j - 1.0);
482
483 //These are in the W rest frame. Need to boost out into
484 //the B frame.
485 p4lepton.set(elepton, 0.0, plepton * sqrt(1.0 - costl * costl),
486 plepton * costl);
487 p4nu.set(plepton, 0.0,
488 -1.0 * plepton * sqrt(1.0 - costl * costl),
489 -1.0 * plepton * costl);
490
491 EvtVector4R boost((m - erho), 0.0, 0.0, 1.0 * prho);
492 p4lepton = boostTo(p4lepton, boost);
493 p4nu = boostTo(p4nu, boost);
494
495 //Now initialize the daughters...
496
497 daughter->init(mId, p4meson);
498 lep->init(l1Id, p4lepton);
499 trino->init(l2Id, p4nu);
500
501 CalcAmp(root_part, amp);
502
503 //Now find the probability at this q2 and cos theta lepton point
504 //and compare to maxfoundprob.
505
506 //Do a little magic to get the probability!!
507 probctl[j] = rho.normalizedProb(amp.getSpinDensity());
508 }
509
510 //probclt contains prob at ctl=-1,0,1.
511 //prob=a+b*ctl+c*ctl^2
512
513 double a = probctl[1];
514 double b = 0.5 * (probctl[2] - probctl[0]);
515 double c = 0.5 * (probctl[2] + probctl[0]) - probctl[1];
516
517 double prob = probctl[0];
518 if (probctl[1] > prob)
519 prob = probctl[1];
520 if (probctl[2] > prob)
521 prob = probctl[2];
522
523 if (fabs(c) > 1e-20) {
524 double ctlx = -0.5 * b / c;
525 if (fabs(ctlx) < 1.0) {
526 double probtmp = a + b * ctlx + c * ctlx * ctlx;
527 if (probtmp > prob)
528 prob = probtmp;
529 }
530 }
531
532 if (prob > maxfoundprob) {
533 maxfoundprob = prob;
534 }
535 }
536 if (EvtPDL::getWidth(mId) <= 0.0) {
537 // if the particle is narrow don't bother with changing the mass.
538 massiter = 4;
539 }
540 }
541 root_part->deleteTree();
542
543 maxfoundprob *= 1.1;
544 return maxfoundprob;
545}
EvtBToDstarlnuNP::CalcAmp
void CalcAmp(EvtParticle *, EvtAmp &)
The method to evaluate the decay amplitude.
Definition EvtBToDstarlnuNP.cc:273
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28

◆ clone()

EvtDecayBase * clone ( )
override

The function which makes a copy of the model.

Definition at line 30 of file EvtBToDstarlnuNP.cc.

31{
32 return new EvtBToDstarlnuNP;
33}

◆ decay()

void decay ( EvtParticle * p)
override

The method to calculate a decay amplitude.

Definition at line 35 of file EvtBToDstarlnuNP.cc.

36{
37 p->initializePhaseSpace(getNDaug(), getDaugs());
38 CalcAmp(p, _amp2);
39}

◆ getName()

std::string getName ( )
override

The function which returns the name of the model.

Definition at line 25 of file EvtBToDstarlnuNP.cc.

26{
27 return "BTODSTARLNUNP";
28}

◆ init()

void init ( )
override

Initialization method.

Definition at line 46 of file EvtBToDstarlnuNP.cc.

47{
48 // First choose format of NP coefficients from the .DEC file
49 // Cartesian(x,y)(0) or Polar(R,phase)(1)
50 int n = getNArg();
51 if (!(n == 0 || (n - 1) % 3 == 0)) {
52 EvtGenReport(EVTGEN_ERROR, "EvtGen")
53 << "Error in parameters in the BTODSTARLNUNP decay model." << endl;
54 EvtGenReport(EVTGEN_ERROR, "EvtGen")
55 << "Will terminate execution!" << endl;
56 ::abort();
57 }
58
59 checkNDaug(3);
60
61 // We expect the parent to be a scalar and
62 // the daughters to be D* l+ nu_l
63
64 checkSpinParent(EvtSpinType::SCALAR);
65 checkSpinDaughter(1, EvtSpinType::DIRAC);
66 checkSpinDaughter(2, EvtSpinType::NEUTRINO);
67
68 EvtId mId = getDaug(0);
69 EvtId IdDst0 = EvtPDL::getId("D*0"), IdaDst0 = EvtPDL::getId("anti-D*0"),
70 IdDstp = EvtPDL::getId("D*+"), IdDstm = EvtPDL::getId("D*-");
71 if (mId != IdDst0 && mId != IdaDst0 && mId != IdDstp && mId != IdDstm) {
72 EvtGenReport(EVTGEN_ERROR, "EvtGen")
73 << "EvtBToDstarlnuNP generator expected a D* 1st daughter, found: "
74 << EvtPDL::name(getDaug(0)) << endl;
75 EvtGenReport(EVTGEN_ERROR, "EvtGen")
76 << "Will terminate execution!" << endl;
77 ::abort();
78 }
79
80 auto getInteger = [this](int i) -> int {
81 double a = getArg(i);
82 if (a - static_cast<int>(a) != 0)
83 {
84 EvtGenReport(EVTGEN_ERROR, "EvtGen")
85 << "Parameters is not integer in the BTODSTARLNUNP decay model: " << i << " " << a << endl;
86 EvtGenReport(EVTGEN_ERROR, "EvtGen")
87 << "Will terminate execution!" << endl;
88 ::abort();
89 }
90 return static_cast<int>(a);
91 };
92
93 if (n > 0) { // parse arguments
94 int i = 0, coordsyst = getInteger(i++);
95 while (i < n) {
96 int id = getInteger(i++); // New Physics cooeficient Id
97 double a0 = getArg(i++);
98 double a1 = getArg(i++);
99 EvtComplex c = (coordsyst) ? EvtComplex(a0 * cos(a1), a0 * sin(a1)) : EvtComplex(a0, a1);
100 if (id == 0) _Cvl = c;
101 if (id == 1) _Cvr = c;
102 if (id == 2) _Csl = c;
103 if (id == 3) _Csr = c;
104 if (id == 4) _cT = c;
105 }
106 }
107}

◆ initProbMax()

void initProbMax ( )
override

The method to evaluate the maximum amplitude.

Definition at line 41 of file EvtBToDstarlnuNP.cc.

42{
43 setProbMax(CalcMaxProb());
44}
EvtBToDstarlnuNP::CalcMaxProb
double CalcMaxProb()
The method to evaluate the maximum decay amplitude.
Definition EvtBToDstarlnuNP.cc:396

Member Data Documentation

◆ _Csl

EvtComplex _Csl
private

(g_S - g_P) strength of pseudoscalar current

Definition at line 57 of file EvtBToDstarlnuNP.h.

◆ _Csr

EvtComplex _Csr
private

(g_S + g_P) strength of pseudoscalar current with constraint _Csl = -_Csr

Definition at line 60 of file EvtBToDstarlnuNP.h.

◆ _cT

EvtComplex _cT
private

C_T – tensor current.

Definition at line 63 of file EvtBToDstarlnuNP.h.

◆ _Cvl

EvtComplex _Cvl
private

Strength of the left handed part of the vector hadronic current.

Definition at line 51 of file EvtBToDstarlnuNP.h.

◆ _Cvr

EvtComplex _Cvr
private

Strength of the Right handed part of the vector hadronic current.

Definition at line 54 of file EvtBToDstarlnuNP.h.

The documentation for this class was generated from the following files: