Belle II Software development
EvtbTosllNPR Class Reference

Description: Implementation of the B->K(*)l^+l^- decays with New Physics contributions and c\bar{c} resonances included at the amplitude level according to the note by Rahul Sinha and Rusa Mandal. More...

#include <EvtbTosllNPR.h>

Inheritance diagram for EvtbTosllNPR:

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.
 
virtual ~EvtbTosllNPR ()
 The destructor to clean up objects created during the initialization.
 

Private Member Functions

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

Private Attributes

EvtComplex m_dc7
 delta C_7eff – addition to NNLO SM value
 
EvtComplex m_dc9
 delta C_9eff – addition to NNLO SM value
 
EvtComplex m_dc10
 delta C_10eff – addition to NNLO SM value
 
EvtComplex m_c7p
 C'_7eff – right hand polarizations.
 
EvtComplex m_c9p
 C'_9eff – right hand polarizations.
 
EvtComplex m_c10p
 c'_10eff – right hand polarizations
 
EvtComplex m_cS
 (C_S - C'_S) – scalar right and left polarizations
 
EvtComplex m_cP
 (C_P - C'_P) – pseudo-scalar right and left polarizations
 
int m_flag {0}
 flag is set nonzero to include resonances
 
std::vector< BW_t * > m_rs
 vector with c\bar{c} resonance lineshapes
 
EvtLinSamplem_ls {0}
 piece-wise interpolation of maximum of the matrix element depend on q^2
 
hvec_t m_reso
 tabulated resonance contribution
 

Detailed Description

Description: Implementation of the B->K(*)l^+l^- decays with New Physics contributions and c\bar{c} resonances included at the amplitude level according to the note by Rahul Sinha and Rusa Mandal.

Definition at line 29 of file EvtbTosllNPR.h.

Constructor & Destructor Documentation

◆ ~EvtbTosllNPR()

~EvtbTosllNPR ( )
virtual

The destructor to clean up objects created during the initialization.

Definition at line 778 of file EvtbTosllNPR.cc.

779{
780 for (BW_t* t : m_rs) delete t;
781 if (m_ls) delete m_ls;
782}
std::vector< BW_t * > m_rs
vector with c\bar{c} resonance lineshapes
EvtLinSample * m_ls
piece-wise interpolation of maximum of the matrix element depend on q^2

Member Function Documentation

◆ CalcAmp()

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

The method to evaluate the decay amplitude.

Definition at line 824 of file EvtbTosllNPR.cc.

825{
826 EvtId dId = parent->getDaug(0)->getId();
827 if (dId == IdKst0 || dId == IdaKst0 || dId == IdKstp || dId == IdKstm) {
828 CalcVAmp(parent, amp);
829 }
830 if (dId == IdK0 || dId == IdaK0 || dId == IdKp || dId == IdKm) {
831 CalcSAmp(parent, amp);
832 }
833}
void CalcSAmp(EvtParticle *, EvtAmp &)
The method to evaluate the scalar kaon decay amplitude.
void CalcVAmp(EvtParticle *, EvtAmp &)
The method to evaluate the vector kaon decay amplitude.

◆ CalcSAmp()

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

The method to evaluate the scalar kaon decay amplitude.

Definition at line 1147 of file EvtbTosllNPR.cc.

1148{
1149 // Add the lepton and neutrino 4 momenta to find q2
1150 EvtId pId = parent->getId();
1151
1152 EvtVector4R q = parent->getDaug(1)->getP4() + parent->getDaug(2)->getP4();
1153 double q2 = q.mass2();
1154 double dmass = parent->getDaug(0)->mass();
1155
1156 double fp, f0, ft; // form factors
1157 getScalarFF(q2, fp, f0, ft);
1158 const EvtVector4R& k = parent->getDaug(0)->getP4();
1159 double pmass = parent->mass();
1160 const EvtVector4R p(pmass, 0.0, 0.0, 0.0);
1161 EvtVector4R pk = p + k;
1162
1163 EvtComplex c7 = -0.304;
1164 EvtComplex c9 = C9(q2, interpol(m_reso, q2));
1165 EvtComplex c10 = -4.103;
1166 c7 += m_dc7;
1167 c9 += m_dc9;
1168 c10 += m_dc10;
1169
1170 double mb = 4.8 /*GeV/c^2*/, ms = 0.093 /*GeV/c^2*/;
1171
1172 int charge1 = EvtPDL::chg3(parent->getDaug(1)->getId());
1173
1174 EvtParticle* lepPos = (charge1 > 0) ? parent->getDaug(1)
1175 : parent->getDaug(2);
1176 EvtParticle* lepNeg = (charge1 < 0) ? parent->getDaug(1)
1177 : parent->getDaug(2);
1178
1179 EvtDiracSpinor lp0(lepPos->spParent(0)), lp1(lepPos->spParent(1));
1180 EvtDiracSpinor lm0(lepNeg->spParent(0)), lm1(lepNeg->spParent(1));
1181
1182 EvtVector4C l11, l12, l21, l22, a11, a12, a21, a22;
1183 EvtComplex s11, s12, s21, s22, p11, p12, p21, p22;
1184
1185 if (pId == IdBm || pId == IdaB0 || pId == IdaBs) {
1186 EvtLeptonVandACurrents(l11, a11, lp0, lm0);
1187 EvtLeptonVandACurrents(l21, a21, lp1, lm0);
1188 EvtLeptonVandACurrents(l12, a12, lp0, lm1);
1189 EvtLeptonVandACurrents(l22, a22, lp1, lm1);
1190
1191 s11 = EvtLeptonSCurrent(lp0, lm0);
1192 p11 = EvtLeptonPCurrent(lp0, lm0);
1193 s21 = EvtLeptonSCurrent(lp1, lm0);
1194 p21 = EvtLeptonPCurrent(lp1, lm0);
1195 s12 = EvtLeptonSCurrent(lp0, lm1);
1196 p12 = EvtLeptonPCurrent(lp0, lm1);
1197 s22 = EvtLeptonSCurrent(lp1, lm1);
1198 p22 = EvtLeptonPCurrent(lp1, lm1);
1199 } else if (pId == IdBp || pId == IdB0 || pId == IdBs) {
1200 EvtLeptonVandACurrents(l11, a11, lp1, lm1);
1201 EvtLeptonVandACurrents(l21, a21, lp0, lm1);
1202 EvtLeptonVandACurrents(l12, a12, lp1, lm0);
1203 EvtLeptonVandACurrents(l22, a22, lp0, lm0);
1204
1205 s11 = EvtLeptonSCurrent(lp1, lm1);
1206 p11 = EvtLeptonPCurrent(lp1, lm1);
1207 s21 = EvtLeptonSCurrent(lp0, lm1);
1208 p21 = EvtLeptonPCurrent(lp0, lm1);
1209 s12 = EvtLeptonSCurrent(lp1, lm0);
1210 p12 = EvtLeptonPCurrent(lp1, lm0);
1211 s22 = EvtLeptonSCurrent(lp0, lm0);
1212 p22 = EvtLeptonPCurrent(lp0, lm0);
1213 } else {
1214 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Wrong lepton number\n";
1215 }
1216 double dm2 = pmass * pmass - dmass * dmass, t0 = dm2 / q2,
1217 t1 = 2 * mb * ft / (pmass + dmass);
1218 c7 += m_c7p;
1219 c9 += m_c9p;
1220 c10 += m_c10p;
1221 EvtVector4C E1 = (c9 * fp + c7 * t1) * pk +
1222 (t0 * (c9 * (f0 - fp) - c7 * t1)) * q;
1223 EvtVector4C E2 = (c10 * fp) * pk + (t0 * (f0 - fp)) * q;
1224 double s = dm2 / (mb - ms) * f0;
1225 amp.vertex(0, 0, l11 * E1 + a11 * E2 + (m_cS * s11 + m_cP * p11) * s);
1226 amp.vertex(0, 1, l12 * E1 + a12 * E2 + (m_cS * s12 + m_cP * p12) * s);
1227 amp.vertex(1, 0, l21 * E1 + a21 * E2 + (m_cS * s21 + m_cP * p21) * s);
1228 amp.vertex(1, 1, l22 * E1 + a22 * E2 + (m_cS * s22 + m_cP * p22) * s);
1229}
EvtComplex m_dc10
delta C_10eff – addition to NNLO SM value
hvec_t m_reso
tabulated resonance contribution
EvtComplex m_cS
(C_S - C'_S) – scalar right and left polarizations
EvtComplex m_cP
(C_P - C'_P) – pseudo-scalar right and left polarizations
EvtComplex m_c9p
C'_9eff – right hand polarizations.
EvtComplex m_dc7
delta C_7eff – addition to NNLO SM value
EvtComplex m_dc9
delta C_9eff – addition to NNLO SM value
EvtComplex m_c7p
C'_7eff – right hand polarizations.
EvtComplex m_c10p
c'_10eff – right hand polarizations

◆ CalcVAmp()

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

The method to evaluate the vector kaon decay amplitude.

Definition at line 960 of file EvtbTosllNPR.cc.

961{
962 // Add the lepton and neutrino 4 momenta to find q2
963 EvtId pId = parent->getId();
964
965 EvtVector4R q = parent->getDaug(1)->getP4() + parent->getDaug(2)->getP4();
966 double q2 = q.mass2();
967 double mesonmass = parent->getDaug(0)->mass();
968 double parentmass = parent->mass();
969
970 double a1, a2, a0, v, t1, t2, t3; // form factors
971 getVectorFF(q2, parentmass, mesonmass, a1, a2, a0, v, t1, t2, t3);
972
973 const EvtVector4R& p4meson = parent->getDaug(0)->getP4();
974 double pmass = parent->mass(), ipmass = 1 / pmass;
975 const EvtVector4R p4b(pmass, 0.0, 0.0, 0.0);
976 EvtVector4R pbhat = p4b * ipmass;
977 EvtVector4R qhat = q * ipmass;
978 EvtVector4R pkstarhat = p4meson * ipmass;
979 EvtVector4R phat = pbhat + pkstarhat;
980
981 EvtComplex c7 = -0.304;
982 EvtComplex c9 = C9(q2, interpol(m_reso, q2));
983 EvtComplex c10 = -4.103;
984 c7 += m_dc7;
985 c9 += m_dc9;
986 c10 += m_dc10;
987
988 EvtComplex I(0.0, 1.0);
989
990 double mb = 4.8 /*GeV/c^2*/ * ipmass, ms = 0.093 /*GeV/c^2*/ * ipmass;
991 double mH = mesonmass * ipmass, oamH = 1 + mH, osmH = 1 - mH,
992 osmH2 = oamH * osmH, iosmH2 = 1 / osmH2; // mhatkstar
993 double is = pmass * pmass / q2; // 1/shat
994 a1 *= oamH;
995 a2 *= osmH;
996 a0 *= 2 * mH;
997 double cs0 = (a1 - a2 - a0) * is;
998 a2 *= iosmH2;
999 v *= 2 / oamH;
1000
1001 EvtComplex a = (c9 + m_c9p) * v + (c7 + m_c7p) * (4 * mb * is * t1);
1002 EvtComplex b = (c9 - m_c9p) * a1 +
1003 (c7 - m_c7p) * (2 * mb * is * osmH2 * t2);
1004 EvtComplex c = (c9 - m_c9p) * a2 +
1005 (c7 - m_c7p) * (2 * mb * (t3 * iosmH2 + t2 * is));
1006 EvtComplex d = (c9 - m_c9p) * cs0 - (c7 - m_c7p) * (2 * mb * is * t3);
1007 EvtComplex e = (c10 + m_c10p) * v;
1008 EvtComplex f = (c10 - m_c10p) * a1;
1009 EvtComplex g = (c10 - m_c10p) * a2;
1010 EvtComplex h = (c10 - m_c10p) * cs0;
1011 double sscale = a0 / (mb + ms);
1012 EvtComplex hs = m_cS * sscale, hp = m_cP * sscale;
1013
1014 int charge1 = EvtPDL::chg3(parent->getDaug(1)->getId());
1015
1016 EvtParticle* lepPos = parent->getDaug(2 - (charge1 > 0));
1017 EvtParticle* lepNeg = parent->getDaug(1 + (charge1 > 0));
1018
1019 EvtDiracSpinor lp0(lepPos->spParent(0)), lp1(lepPos->spParent(1));
1020 EvtDiracSpinor lm0(lepNeg->spParent(0)), lm1(lepNeg->spParent(1));
1021
1022 EvtVector4C l11, l12, l21, l22, a11, a12, a21, a22;
1023 EvtComplex s11, s12, s21, s22, p11, p12, p21, p22;
1024 EvtTensor4C tt0 = asymProd(pbhat, pkstarhat);
1025
1026 EvtTensor4C T1(tt0), T2(tt0);
1027 const EvtTensor4C& gmn = EvtTensor4C::g();
1028 EvtTensor4C tt1(EvtGenFunctions::directProd(pbhat, phat));
1029 EvtTensor4C tt2(EvtGenFunctions::directProd(pbhat, qhat));
1030
1031 b *= I;
1032 c *= I;
1033 d *= I;
1034 f *= I;
1035 g *= I;
1036 h *= I;
1037 if (pId == IdBm || pId == IdaB0 || pId == IdaBs) {
1038 T1 *= a;
1039 addScaled(T1, -b, gmn);
1040 addScaled(T1, c, tt1);
1041 addScaled(T1, d, tt2);
1042 T2 *= e;
1043 addScaled(T2, -f, gmn);
1044 addScaled(T2, g, tt1);
1045 addScaled(T2, h, tt2);
1046
1047 EvtLeptonVandACurrents(l11, a11, lp0, lm0);
1048 EvtLeptonVandACurrents(l21, a21, lp1, lm0);
1049 EvtLeptonVandACurrents(l12, a12, lp0, lm1);
1050 EvtLeptonVandACurrents(l22, a22, lp1, lm1);
1051
1052 s11 = EvtLeptonSCurrent(lp0, lm0);
1053 p11 = EvtLeptonPCurrent(lp0, lm0);
1054 s21 = EvtLeptonSCurrent(lp1, lm0);
1055 p21 = EvtLeptonPCurrent(lp1, lm0);
1056 s12 = EvtLeptonSCurrent(lp0, lm1);
1057 p12 = EvtLeptonPCurrent(lp0, lm1);
1058 s22 = EvtLeptonSCurrent(lp1, lm1);
1059 p22 = EvtLeptonPCurrent(lp1, lm1);
1060 } else if (pId == IdBp || pId == IdB0 || pId == IdBs) {
1061 T1 *= -a;
1062 addScaled(T1, -b, gmn);
1063 addScaled(T1, c, tt1);
1064 addScaled(T1, d, tt2);
1065 T2 *= -e;
1066 addScaled(T2, -f, gmn);
1067 addScaled(T2, g, tt1);
1068 addScaled(T2, h, tt2);
1069
1070 EvtLeptonVandACurrents(l11, a11, lp1, lm1);
1071 EvtLeptonVandACurrents(l21, a21, lp0, lm1);
1072 EvtLeptonVandACurrents(l12, a12, lp1, lm0);
1073 EvtLeptonVandACurrents(l22, a22, lp0, lm0);
1074
1075 s11 = EvtLeptonSCurrent(lp1, lm1);
1076 p11 = EvtLeptonPCurrent(lp1, lm1);
1077 s21 = EvtLeptonSCurrent(lp0, lm1);
1078 p21 = EvtLeptonPCurrent(lp0, lm1);
1079 s12 = EvtLeptonSCurrent(lp1, lm0);
1080 p12 = EvtLeptonPCurrent(lp1, lm0);
1081 s22 = EvtLeptonSCurrent(lp0, lm0);
1082 p22 = EvtLeptonPCurrent(lp0, lm0);
1083 } else {
1084 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Wrong lepton number\n";
1085 T1.zero();
1086 T2.zero(); // Set all tensor terms to zero.
1087 }
1088 for (int i = 0; i < 3; i++) {
1089 EvtVector4C eps = parent->getDaug(0)->epsParent(i).conj();
1090
1091 EvtVector4C E1 = T1.cont1(eps);
1092 EvtVector4C E2 = T2.cont1(eps);
1093
1094 EvtComplex epsq = I * (eps * q), epsqs = epsq * hs, epsqp = epsq * hp;
1095
1096 amp.vertex(i, 0, 0, l11 * E1 + a11 * E2 - s11 * epsqs - p11 * epsqp);
1097 amp.vertex(i, 0, 1, l12 * E1 + a12 * E2 - s12 * epsqs - p12 * epsqp);
1098 amp.vertex(i, 1, 0, l21 * E1 + a21 * E2 - s21 * epsqs - p21 * epsqp);
1099 amp.vertex(i, 1, 1, l22 * E1 + a22 * E2 - s22 * epsqs - p22 * epsqp);
1100 }
1101}

◆ clone()

EvtDecayBase * clone ( )
override

The function which makes a copy of the model.

Definition at line 46 of file EvtbTosllNPR.cc.

47{
48 return new EvtbTosllNPR;
49}

◆ decay()

void decay ( EvtParticle * p)
override

The method to calculate a decay amplitude.

Definition at line 166 of file EvtbTosllNPR.cc.

167{
168 static EvtVector4R p4[3];
169 static double mass[3];
170 double m_b = p->mass();
171 for (int i = 0; i < 3; i++)
172 mass[i] = p->getDaug(i)->mass();
173 EvtId* daughters = getDaugs();
174 double weight = PhaseSpacePole2(m_b, mass[2], mass[1], mass[0], p4, m_ls);
175 p->getDaug(0)->init(daughters[0], p4[2]);
176 p->getDaug(1)->init(daughters[1], p4[1]);
177 p->getDaug(2)->init(daughters[2], p4[0]);
178
179 setWeight(weight);
180 CalcAmp(p, _amp2);
181}
void CalcAmp(EvtParticle *, EvtAmp &)
The method to evaluate the decay amplitude.

◆ getName()

std::string getName ( )
override

The function which returns the name of the model.

Definition at line 41 of file EvtbTosllNPR.cc.

42{
43 return "BTOSLLNPR";
44}

◆ init()

void init ( )
override

Initialization method.

Definition at line 614 of file EvtbTosllNPR.cc.

615{
616 if (cafirst) {
617 cafirst = false;
618 IdB0 = EvtPDL::getId("B0");
619 IdaB0 = EvtPDL::getId("anti-B0");
620 IdBp = EvtPDL::getId("B+");
621 IdBm = EvtPDL::getId("B-");
622 IdBs = EvtPDL::getId("B_s0");
623 IdaBs = EvtPDL::getId("anti-B_s0");
624 IdRhop = EvtPDL::getId("rho+");
625 IdRhom = EvtPDL::getId("rho-");
626 IdRho0 = EvtPDL::getId("rho0");
627 IdOmega = EvtPDL::getId("omega");
628 IdKst0 = EvtPDL::getId("K*0");
629 IdaKst0 = EvtPDL::getId("anti-K*0");
630 IdKstp = EvtPDL::getId("K*+");
631 IdKstm = EvtPDL::getId("K*-");
632 IdK0 = EvtPDL::getId("K0");
633 IdaK0 = EvtPDL::getId("anti-K0");
634 IdKp = EvtPDL::getId("K+");
635 IdKm = EvtPDL::getId("K-");
636 }
637
638 // First choose format of NP coefficients from the .DEC file
639 // Cartesian(x,y)(0) or Polar(R,phase)(1)
640 int n = getNArg();
641 checkNDaug(3);
642
643 //We expect the parent to be a scalar
644 //and the daughters to be K(*) lepton+ lepton-
645
646 checkSpinParent(EvtSpinType::SCALAR);
647 checkSpinDaughter(1, EvtSpinType::DIRAC);
648 checkSpinDaughter(2, EvtSpinType::DIRAC);
649
650 EvtId mId = getDaug(0);
651 if (mId != IdKst0 && mId != IdaKst0 && mId != IdKstp && mId != IdKstm &&
652 mId != IdK0 && mId != IdaK0 && mId != IdKp && mId != IdKm) {
653 EvtGenReport(EVTGEN_ERROR, "EvtGen")
654 << "EvtbTosllNPR generator expected a K(*) 1st daughter, found: "
655 << EvtPDL::name(getDaug(0)) << endl;
656 EvtGenReport(EVTGEN_ERROR, "EvtGen")
657 << "Will terminate execution!" << endl;
658 ::abort();
659 }
660
661 auto getInteger = [this](int i) -> int {
662 double a = getArg(i);
663 if (a - static_cast<int>(a) != 0)
664 {
665 EvtGenReport(EVTGEN_ERROR, "EvtGen")
666 << "Parameters is not integer in the BTOSLLNPR decay model: " << i
667 << " " << a << endl;
668 EvtGenReport(EVTGEN_ERROR, "EvtGen")
669 << "Will terminate execution!" << endl;
670 ::abort();
671 }
672 return static_cast<int>(a);
673 };
674
675 double phi = 0.0, n1 = 0.0, d1 = 1.0, eta = 0.0, theta_Jpsi = 0.0, theta_psi2S = 0.0, Nr = 1.0;
676 if (n > 0) { // parse arguments
677 int i = 0, coordsyst = getInteger(i++);
678 auto getcomplex = [this, &i, coordsyst]() -> EvtComplex {
679 double a0 = getArg(i++);
680 double a1 = getArg(i++);
681 return (coordsyst) ? EvtComplex(a0 * cos(a1), a0 * sin(a1))
682 : EvtComplex(a0, a1);
683 };
684 auto getreal = [this, &i]() -> double { return getArg(i++); };
685 while (i < n) {
686 int id = getInteger(i++); // New Physics cooeficient Id
687 if (id == 0)
688 m_dc7 = getcomplex(); // delta C_7eff
689 if (id == 1)
690 m_dc9 = getcomplex(); // delta C_9eff
691 if (id == 2)
692 m_dc10 = getcomplex(); // delta C_10eff
693 if (id == 3)
694 m_c7p = getcomplex(); // C'_7eff -- right hand polarizations
695 if (id == 4)
696 m_c9p = getcomplex(); // C'_9eff -- right hand polarizations
697 if (id == 5)
698 m_c10p = getcomplex(); // c'_10eff -- right hand polarizations
699 if (id == 6)
700 m_cS = getcomplex(); // (C_S - C'_S) -- scalar right and left polarizations
701 if (id == 7)
702 m_cP = getcomplex(); // (C_P - C'_P) -- pseudo-scalar right and left polarizations
703 if (id == 10)
704 m_flag = getInteger(i++); // include resonances or not
705 if (id == 11)
706 phi = getreal(); // phase of the non-factorizable contribution
707 if (id == 12)
708 n1 = getreal(); // amplitude of the non-factorizable contribution
709 if (id == 13)
710 d1 = getreal(); // width of the non-factorizable contribution
711 if (id == 14)
712 eta = getreal(); // Eq. 47
713 if (id == 15)
714 theta_Jpsi = getreal(); // Eq. 47
715 if (id == 16)
716 theta_psi2S = getreal(); // Eq. 47
717 if (id == 17)
718 Nr = getreal(); // Eq. 47
719 }
720 }
721 m_ls = new EvtLinSample;
722
723 if (!m_flag) return;
724 // BW_t jpsi(3.0969, 0.0926 / 1000, 0.0926 / 1000 * 0.877, 0.05971);
725 // BW_t psi2s(3.6861, 0.294 / 1000, 0.294 / 1000 * 0.9785, 0.00793);
726 // BW_t psi3770(3773.7 / 1000, 27.2 / 1000, 27.2 / 1000, 1e-5);
727 // BW_t psi4040(4039 / 1000., 80 / 1000., 80 / 1000., 1.07e-5);
728 // BW_t psi4160(4191 / 1000., 70 / 1000., 70 / 1000., 6.9e-6);
729 // BW_t psi4230(4220 / 1000., 60 / 1000., 60 / 1000., 1e-5);
730
731 TParticlePDG* p_jpsi = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("J/psi");
732 // B2WARNING("extracting leptonic width"<<std::flush);
733 // TDecayChannel *c_jpsi = p_jpsi->DecayChannel(1);
734 // B2WARNING("Br = "<<c_jpsi->BranchingRatio());
735 BW_t* jpsi = new BW_t(p_jpsi->Mass(), p_jpsi->Width(), p_jpsi->Width() * 0.877, 0.05971);
736 m_rs.push_back(jpsi);
737
738 TParticlePDG* p_psi2s = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("psi(2S)");
739 BW_t* psi2s = new BW_t(p_psi2s->Mass(), p_psi2s->Width(), p_psi2s->Width() * 0.9785, 0.00793);
740 m_rs.push_back(psi2s);
741
742 TParticlePDG* p_psi3770 = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("psi(3770)");
743 BW_t* psi3770 = new BW_t(p_psi3770->Mass(), p_psi3770->Width(), p_psi3770->Width(), 1e-5);
744 m_rs.push_back(psi3770);
745
746 TParticlePDG* p_psi4040 = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("psi(4040)");
747 BW_t* psi4040 = new BW_t(p_psi4040->Mass(), p_psi4040->Width(), p_psi4040->Width(), 1.07e-5);
748 m_rs.push_back(psi4040);
749
750 TParticlePDG* p_psi4160 = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("psi(4160)");
751 BW_t* psi4160 = new BW_t(p_psi4160->Mass(), p_psi4160->Width(), p_psi4160->Width(), 6.9e-6);
752 m_rs.push_back(psi4160);
753
754 BW_t* psi4230 = new BW_t(4220 / 1000., 60 / 1000., 60 / 1000., 1e-5);
755 m_rs.push_back(psi4230);
756
757 std::vector<double> v = getgrid(m_rs);
758 std::complex<double> eiphi = exp(1i * phi), pJpsi = Nr * exp(1i * theta_Jpsi), ppsi2S = Nr * exp(1i * theta_psi2S);
759
760 const double C1 = -0.257, C2 = 1.009, C3 = -0.005, C5 = 0;
761 double sc = 1.0 / (4.0 / 3.0 * C1 + C2 + 6.0 * C3 + 60.0 * C5);
762 const double m_b = 4.8, m_c = 1.3;
763 for (auto t : v) {
764 // Eq. 13, 14, and 15 in Emi's note. For the default theta_Jpsi=0 and theta_psi2S=0 it is identical to the Rahul's dispersion relation
765 double re = -8.0 / 9.0 * log(m_c / m_b) - 4.0 / 9.0 +
766 t / 3.0 * DRInt(t, real(pJpsi), real(ppsi2S), Nr, m_rs) - M_PI / 3 * uR(t, imag(pJpsi), imag(ppsi2S), 0, m_rs);
767 double im =
768 t / 3.0 * DRInt(t, imag(pJpsi), imag(ppsi2S), 0, m_rs) + M_PI / 3 * uR(t, real(pJpsi), real(ppsi2S), Nr, m_rs);
769 std::complex<double> r(re, im);
770 if (eta == 0.0) {
771 double z = t - jpsi->Mv2;
772 r += eiphi * (n1 * sc / (z * z + d1)); // bump under J/psi as a non-factorizable contribution
773 }
774 m_reso.push_back(std::make_pair(t, r));
775 }
776}
static EvtGenDatabasePDG * Instance()
Instance method that loads the EvtGen table.
int m_flag
flag is set nonzero to include resonances
double Mv2
resonance mass squared

◆ initProbMax()

void initProbMax ( )
override

The method to evaluate the maximum amplitude.

Definition at line 484 of file EvtbTosllNPR.cc.

485{
486 EvtId pId = getParentId(), mId = getDaug(0), l1Id = getDaug(1), l2Id = getDaug(2);
487
488 std::vector<double> v;
489 std::vector<EvtPointf> pp;
490
491 EvtScalarParticle* scalar_part;
492 EvtParticle* root_part;
493
494 scalar_part = new EvtScalarParticle;
495
496 //includge to avoid generating random numbers!
497 scalar_part->noLifeTime();
498
499 EvtVector4R p_init;
500
501 p_init.set(EvtPDL::getMass(pId), 0.0, 0.0, 0.0);
502 scalar_part->init(pId, p_init);
503 root_part = (EvtParticle*)scalar_part;
504 root_part->setDiagonalSpinDensity();
505
506 EvtParticle* daughter, *lep1, *lep2;
507
508 EvtAmp amp;
509
510 EvtId listdaug[3];
511 listdaug[0] = mId;
512 listdaug[1] = l1Id;
513 listdaug[2] = l2Id;
514
515 amp.init(pId, 3, listdaug);
516
517 root_part->makeDaughters(3, listdaug);
518 daughter = root_part->getDaug(0);
519 lep1 = root_part->getDaug(1);
520 lep2 = root_part->getDaug(2);
521
522 //cludge to avoid generating random numbers!
523 daughter->noLifeTime();
524 lep1->noLifeTime();
525 lep2->noLifeTime();
526
527 //Initial particle is unpolarized, well it is a scalar so it is trivial
528 EvtSpinDensity rho;
529 rho.setDiag(root_part->getSpinStates());
530
531 double mass[3];
532
533 double m = root_part->mass();
534
535 EvtVector4R p4meson, p4lepton1, p4lepton2, p4w;
536
537 double maxprobfound = 0.0;
538
539 mass[1] = EvtPDL::getMeanMass(l1Id);
540 mass[2] = EvtPDL::getMeanMass(l2Id);
541 std::vector<double> mH;
542 mH.push_back(EvtPDL::getMeanMass(mId));
543 //if the particle is narrow dont bother with changing the mass.
544 double g = EvtPDL::getWidth(mId);
545 if (g > 0) {
546 mH.push_back(EvtPDL::getMinMass(mId));
547 mH.push_back(
548 std::min(EvtPDL::getMaxMass(mId), m - mass[1] - mass[2]));
549 mH.push_back(mH.front() - g);
550 mH.push_back(mH.front() + g);
551 }
552
553 double q2min = (mass[1] + mass[2]) * (mass[1] + mass[2]);
554
555 double m0 = EvtPDL::getMinMass(mId);
556 double q2max = (m - m0) * (m - m0);
557 v = getgrid(m_rs);
558 m0 = mH[0];
559 //loop over q2
560 for (double q2 : v) {
561 // want to avoid picking up the tail of the photon propagator
562 if (!(q2min <= q2 && q2 < q2max))
563 continue;
564 double Erho = (m * m + m0 * m0 - q2) / (2.0 * m);
565 if (Erho < m0) {
566 m0 = EvtPDL::getMinMass(mId);
567 Erho = (m * m + m0 * m0 - q2) / (2.0 * m);
568 }
569 double Prho = sqrt((Erho - m0) * (Erho + m0));
570
571 p4meson.set(Erho, 0.0, 0.0, -Prho);
572 p4w.set(m - Erho, 0.0, 0.0, Prho);
573
574 //This is in the W rest frame
575 double Elep = sqrt(q2) * 0.5,
576 Plep = sqrt((Elep - mass[1]) * (Elep + mass[1]));
577
578 const int nj = 3 + 2 + 4 + 8 + 32;
579 double cmin = -1, cmax = 1, dc = (cmax - cmin) / (nj - 1);
580 double maxprob = 0;
581 for (int j = 0; j < nj; j++) {
582 double c = cmin + dc * j;
583
584 //These are in the W rest frame. Need to boost out into the B frame.
585 double Py = Plep * sqrt(1.0 - c * c), Pz = Plep * c;
586 p4lepton1.set(Elep, 0.0, Py, Pz);
587 p4lepton2.set(Elep, 0.0, -Py, -Pz);
588
589 p4lepton1 = boostTo(p4lepton1, p4w);
590 p4lepton2 = boostTo(p4lepton2, p4w);
591
592 //Now initialize the daughters...
593 daughter->init(mId, p4meson);
594 lep1->init(l1Id, p4lepton1);
595 lep2->init(l2Id, p4lepton2);
596 CalcAmp(root_part, amp);
597 double prob = rho.normalizedProb(amp.getSpinDensity());
598 maxprob = std::max(maxprob, prob);
599 }
600 pp.push_back({ (float)q2, (float)maxprob });
601 maxprobfound = std::max(maxprobfound, maxprob);
602 }
603 root_part->deleteTree();
604
605 m_ls->init(pp);
606 maxprobfound = 4;
607 setProbMax(maxprobfound);
608}
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28

Member Data Documentation

◆ m_c10p

EvtComplex m_c10p
private

c'_10eff – right hand polarizations

Definition at line 79 of file EvtbTosllNPR.h.

◆ m_c7p

EvtComplex m_c7p
private

C'_7eff – right hand polarizations.

Definition at line 73 of file EvtbTosllNPR.h.

◆ m_c9p

EvtComplex m_c9p
private

C'_9eff – right hand polarizations.

Definition at line 76 of file EvtbTosllNPR.h.

◆ m_cP

EvtComplex m_cP
private

(C_P - C'_P) – pseudo-scalar right and left polarizations

Definition at line 85 of file EvtbTosllNPR.h.

◆ m_cS

EvtComplex m_cS
private

(C_S - C'_S) – scalar right and left polarizations

Definition at line 82 of file EvtbTosllNPR.h.

◆ m_dc10

EvtComplex m_dc10
private

delta C_10eff – addition to NNLO SM value

Definition at line 70 of file EvtbTosllNPR.h.

◆ m_dc7

EvtComplex m_dc7
private

delta C_7eff – addition to NNLO SM value

Definition at line 64 of file EvtbTosllNPR.h.

◆ m_dc9

EvtComplex m_dc9
private

delta C_9eff – addition to NNLO SM value

Definition at line 67 of file EvtbTosllNPR.h.

◆ m_flag

int m_flag {0}
private

flag is set nonzero to include resonances

Definition at line 88 of file EvtbTosllNPR.h.

88{0};

◆ m_ls

EvtLinSample* m_ls {0}
private

piece-wise interpolation of maximum of the matrix element depend on q^2

Definition at line 94 of file EvtbTosllNPR.h.

94{0};

◆ m_reso

hvec_t m_reso
private

tabulated resonance contribution

Definition at line 97 of file EvtbTosllNPR.h.

◆ m_rs

std::vector<BW_t*> m_rs
private

vector with c\bar{c} resonance lineshapes

Definition at line 91 of file EvtbTosllNPR.h.


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