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

EvtLinSample *	m_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 790 of file EvtbTosllNPR.cc.

{
  for (BW_t* t : m_rs) delete t;
  if (m_ls) delete m_ls;
}

Member Function Documentation

◆ CalcAmp()

void CalcAmp	(	EvtParticle *	parent,
		EvtAmp &	amp )

private

The method to evaluate the decay amplitude.

Definition at line 836 of file EvtbTosllNPR.cc.

{
  EvtId dId = parent->getDaug(0)->getId();
  if (dId == IdKst0 || dId == IdaKst0 || dId == IdKstp || dId == IdKstm) {
    CalcVAmp(parent, amp);
  }
  if (dId == IdK0 || dId == IdaK0 || dId == IdKp || dId == IdKm) {
    CalcSAmp(parent, amp);
  }
}

◆ CalcSAmp()

void CalcSAmp	(	EvtParticle *	parent,
		EvtAmp &	amp )

private

The method to evaluate the scalar kaon decay amplitude.

Definition at line 1159 of file EvtbTosllNPR.cc.

{
  // Add the lepton and neutrino 4 momenta to find q2
  EvtId pId = parent->getId();
 
  EvtVector4R q = parent->getDaug(1)->getP4() + parent->getDaug(2)->getP4();
  double q2 = q.mass2();
  double dmass = parent->getDaug(0)->mass();
 
  double fp, f0, ft;    // form factors
  getScalarFF(q2, fp, f0, ft);
  const EvtVector4R& k = parent->getDaug(0)->getP4();
  double pmass = parent->mass();
  const EvtVector4R p(pmass, 0.0, 0.0, 0.0);
  EvtVector4R pk = p + k;
 
  EvtComplex c7 = -0.304;
  EvtComplex c9 = C9(q2, interpol(m_reso, q2));
  EvtComplex c10 = -4.103;
  c7 += m_dc7;
  c9 += m_dc9;
  c10 += m_dc10;
 
  double mb = 4.8 /*GeV/c^2*/, ms = 0.093 /*GeV/c^2*/;
 
  int charge1 = EvtPDL::chg3(parent->getDaug(1)->getId());
 
  EvtParticle* lepPos = (charge1 > 0) ? parent->getDaug(1)
                        : parent->getDaug(2);
  EvtParticle* lepNeg = (charge1 < 0) ? parent->getDaug(1)
                        : parent->getDaug(2);
 
  EvtDiracSpinor lp0(lepPos->spParent(0)), lp1(lepPos->spParent(1));
  EvtDiracSpinor lm0(lepNeg->spParent(0)), lm1(lepNeg->spParent(1));
 
  EvtVector4C l11, l12, l21, l22, a11, a12, a21, a22;
  EvtComplex s11, s12, s21, s22, p11, p12, p21, p22;
 
  if (pId == IdBm || pId == IdaB0 || pId == IdaBs) {
    EvtLeptonVandACurrents(l11, a11, lp0, lm0);
    EvtLeptonVandACurrents(l21, a21, lp1, lm0);
    EvtLeptonVandACurrents(l12, a12, lp0, lm1);
    EvtLeptonVandACurrents(l22, a22, lp1, lm1);
 
    s11 = EvtLeptonSCurrent(lp0, lm0);
    p11 = EvtLeptonPCurrent(lp0, lm0);
    s21 = EvtLeptonSCurrent(lp1, lm0);
    p21 = EvtLeptonPCurrent(lp1, lm0);
    s12 = EvtLeptonSCurrent(lp0, lm1);
    p12 = EvtLeptonPCurrent(lp0, lm1);
    s22 = EvtLeptonSCurrent(lp1, lm1);
    p22 = EvtLeptonPCurrent(lp1, lm1);
  } else if (pId == IdBp || pId == IdB0 || pId == IdBs) {
    EvtLeptonVandACurrents(l11, a11, lp1, lm1);
    EvtLeptonVandACurrents(l21, a21, lp0, lm1);
    EvtLeptonVandACurrents(l12, a12, lp1, lm0);
    EvtLeptonVandACurrents(l22, a22, lp0, lm0);
 
    s11 = EvtLeptonSCurrent(lp1, lm1);
    p11 = EvtLeptonPCurrent(lp1, lm1);
    s21 = EvtLeptonSCurrent(lp0, lm1);
    p21 = EvtLeptonPCurrent(lp0, lm1);
    s12 = EvtLeptonSCurrent(lp1, lm0);
    p12 = EvtLeptonPCurrent(lp1, lm0);
    s22 = EvtLeptonSCurrent(lp0, lm0);
    p22 = EvtLeptonPCurrent(lp0, lm0);
  } else {
    EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Wrong lepton number\n";
  }
  double dm2 = pmass * pmass - dmass * dmass, t0 = dm2 / q2,
         t1 = 2 * mb * ft / (pmass + dmass);
  c7 += m_c7p;
  c9 += m_c9p;
  c10 += m_c10p;
  EvtVector4C E1 = (c9 * fp + c7 * t1) * pk +
                   (t0 * (c9 * (f0 - fp) - c7 * t1)) * q;
  EvtVector4C E2 = (c10 * fp) * pk + (t0 * (f0 - fp)) * q;
  double s = dm2 / (mb - ms) * f0;
  amp.vertex(0, 0, l11 * E1 + a11 * E2 + (m_cS * s11 + m_cP * p11) * s);
  amp.vertex(0, 1, l12 * E1 + a12 * E2 + (m_cS * s12 + m_cP * p12) * s);
  amp.vertex(1, 0, l21 * E1 + a21 * E2 + (m_cS * s21 + m_cP * p21) * s);
  amp.vertex(1, 1, l22 * E1 + a22 * E2 + (m_cS * s22 + m_cP * p22) * s);
}

◆ CalcVAmp()

void CalcVAmp	(	EvtParticle *	parent,
		EvtAmp &	amp )

private

The method to evaluate the vector kaon decay amplitude.

Definition at line 972 of file EvtbTosllNPR.cc.

{
  // Add the lepton and neutrino 4 momenta to find q2
  EvtId pId = parent->getId();
 
  EvtVector4R q = parent->getDaug(1)->getP4() + parent->getDaug(2)->getP4();
  double q2 = q.mass2();
  double mesonmass = parent->getDaug(0)->mass();
  double parentmass = parent->mass();
 
  double a1, a2, a0, v, t1, t2, t3;    // form factors
  getVectorFF(q2, parentmass, mesonmass, a1, a2, a0, v, t1, t2, t3);
 
  const EvtVector4R& p4meson = parent->getDaug(0)->getP4();
  double pmass = parent->mass(), ipmass = 1 / pmass;
  const EvtVector4R p4b(pmass, 0.0, 0.0, 0.0);
  EvtVector4R pbhat = p4b * ipmass;
  EvtVector4R qhat = q * ipmass;
  EvtVector4R pkstarhat = p4meson * ipmass;
  EvtVector4R phat = pbhat + pkstarhat;
 
  EvtComplex c7 = -0.304;
  EvtComplex c9 = C9(q2, interpol(m_reso, q2));
  EvtComplex c10 = -4.103;
  c7 += m_dc7;
  c9 += m_dc9;
  c10 += m_dc10;
 
  EvtComplex I(0.0, 1.0);
 
  double mb = 4.8 /*GeV/c^2*/ * ipmass, ms = 0.093 /*GeV/c^2*/ * ipmass;
  double mH = mesonmass * ipmass, oamH = 1 + mH, osmH = 1 - mH,
         osmH2 = oamH * osmH, iosmH2 = 1 / osmH2;    // mhatkstar
  double is = pmass * pmass / q2;                    // 1/shat
  a1 *= oamH;
  a2 *= osmH;
  a0 *= 2 * mH;
  double cs0 = (a1 - a2 - a0) * is;
  a2 *= iosmH2;
  v *= 2 / oamH;
 
  EvtComplex a = (c9 + m_c9p) * v + (c7 + m_c7p) * (4 * mb * is * t1);
  EvtComplex b = (c9 - m_c9p) * a1 +
                 (c7 - m_c7p) * (2 * mb * is * osmH2 * t2);
  EvtComplex c = (c9 - m_c9p) * a2 +
                 (c7 - m_c7p) * (2 * mb * (t3 * iosmH2 + t2 * is));
  EvtComplex d = (c9 - m_c9p) * cs0 - (c7 - m_c7p) * (2 * mb * is * t3);
  EvtComplex e = (c10 + m_c10p) * v;
  EvtComplex f = (c10 - m_c10p) * a1;
  EvtComplex g = (c10 - m_c10p) * a2;
  EvtComplex h = (c10 - m_c10p) * cs0;
  double sscale = a0 / (mb + ms);
  EvtComplex hs = m_cS * sscale, hp = m_cP * sscale;
 
  int charge1 = EvtPDL::chg3(parent->getDaug(1)->getId());
 
  EvtParticle* lepPos = parent->getDaug(2 - (charge1 > 0));
  EvtParticle* lepNeg = parent->getDaug(1 + (charge1 > 0));
 
  EvtDiracSpinor lp0(lepPos->spParent(0)), lp1(lepPos->spParent(1));
  EvtDiracSpinor lm0(lepNeg->spParent(0)), lm1(lepNeg->spParent(1));
 
  EvtVector4C l11, l12, l21, l22, a11, a12, a21, a22;
  EvtComplex s11, s12, s21, s22, p11, p12, p21, p22;
  EvtTensor4C tt0 = asymProd(pbhat, pkstarhat);
 
  EvtTensor4C T1(tt0), T2(tt0);
  const EvtTensor4C& gmn = EvtTensor4C::g();
  EvtTensor4C tt1(EvtGenFunctions::directProd(pbhat, phat));
  EvtTensor4C tt2(EvtGenFunctions::directProd(pbhat, qhat));
 
  b *= I;
  c *= I;
  d *= I;
  f *= I;
  g *= I;
  h *= I;
  if (pId == IdBm || pId == IdaB0 || pId == IdaBs) {
    T1 *= a;
    addScaled(T1, -b, gmn);
    addScaled(T1, c, tt1);
    addScaled(T1, d, tt2);
    T2 *= e;
    addScaled(T2, -f, gmn);
    addScaled(T2, g, tt1);
    addScaled(T2, h, tt2);
 
    EvtLeptonVandACurrents(l11, a11, lp0, lm0);
    EvtLeptonVandACurrents(l21, a21, lp1, lm0);
    EvtLeptonVandACurrents(l12, a12, lp0, lm1);
    EvtLeptonVandACurrents(l22, a22, lp1, lm1);
 
    s11 = EvtLeptonSCurrent(lp0, lm0);
    p11 = EvtLeptonPCurrent(lp0, lm0);
    s21 = EvtLeptonSCurrent(lp1, lm0);
    p21 = EvtLeptonPCurrent(lp1, lm0);
    s12 = EvtLeptonSCurrent(lp0, lm1);
    p12 = EvtLeptonPCurrent(lp0, lm1);
    s22 = EvtLeptonSCurrent(lp1, lm1);
    p22 = EvtLeptonPCurrent(lp1, lm1);
  } else if (pId == IdBp || pId == IdB0 || pId == IdBs) {
    T1 *= -a;
    addScaled(T1, -b, gmn);
    addScaled(T1, c, tt1);
    addScaled(T1, d, tt2);
    T2 *= -e;
    addScaled(T2, -f, gmn);
    addScaled(T2, g, tt1);
    addScaled(T2, h, tt2);
 
    EvtLeptonVandACurrents(l11, a11, lp1, lm1);
    EvtLeptonVandACurrents(l21, a21, lp0, lm1);
    EvtLeptonVandACurrents(l12, a12, lp1, lm0);
    EvtLeptonVandACurrents(l22, a22, lp0, lm0);
 
    s11 = EvtLeptonSCurrent(lp1, lm1);
    p11 = EvtLeptonPCurrent(lp1, lm1);
    s21 = EvtLeptonSCurrent(lp0, lm1);
    p21 = EvtLeptonPCurrent(lp0, lm1);
    s12 = EvtLeptonSCurrent(lp1, lm0);
    p12 = EvtLeptonPCurrent(lp1, lm0);
    s22 = EvtLeptonSCurrent(lp0, lm0);
    p22 = EvtLeptonPCurrent(lp0, lm0);
  } else {
    EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Wrong lepton number\n";
    T1.zero();
    T2.zero();    // Set all tensor terms to zero.
  }
  for (int i = 0; i < 3; i++) {
    EvtVector4C eps = parent->getDaug(0)->epsParent(i).conj();
 
    EvtVector4C E1 = T1.cont1(eps);
    EvtVector4C E2 = T2.cont1(eps);
 
    EvtComplex epsq = I * (eps * q), epsqs = epsq * hs, epsqp = epsq * hp;
 
    amp.vertex(i, 0, 0, l11 * E1 + a11 * E2 - s11 * epsqs - p11 * epsqp);
    amp.vertex(i, 0, 1, l12 * E1 + a12 * E2 - s12 * epsqs - p12 * epsqp);
    amp.vertex(i, 1, 0, l21 * E1 + a21 * E2 - s21 * epsqs - p21 * epsqp);
    amp.vertex(i, 1, 1, l22 * E1 + a22 * E2 - s22 * epsqs - p22 * epsqp);
  }
}

◆ clone()

EvtDecayBase * clone ( )

override

The function which makes a copy of the model.

Definition at line 46 of file EvtbTosllNPR.cc.

{
  return new EvtbTosllNPR;
}

◆ decay()

void decay ( EvtParticle * p )

override

The method to calculate a decay amplitude.

Definition at line 166 of file EvtbTosllNPR.cc.

{
  static EvtVector4R p4[3];
  static double mass[3];
  double m_b = p->mass();
  for (int i = 0; i < 3; i++)
    mass[i] = p->getDaug(i)->mass();
  EvtId* daughters = getDaugs();
  double weight = PhaseSpacePole2(m_b, mass[2], mass[1], mass[0], p4, m_ls);
  p->getDaug(0)->init(daughters[0], p4[2]);
  p->getDaug(1)->init(daughters[1], p4[1]);
  p->getDaug(2)->init(daughters[2], p4[0]);
 
  setWeight(weight);
  CalcAmp(p, _amp2);
}

◆ getName()

std::string getName ( )

override

The function which returns the name of the model.

Definition at line 41 of file EvtbTosllNPR.cc.

{
  return "BTOSLLNPR";
}

◆ init()

void init ( )

override

Initialization method.

Definition at line 626 of file EvtbTosllNPR.cc.

{
  if (cafirst) {
    cafirst = false;
    IdB0 = EvtPDL::getId("B0");
    IdaB0 = EvtPDL::getId("anti-B0");
    IdBp = EvtPDL::getId("B+");
    IdBm = EvtPDL::getId("B-");
    IdBs = EvtPDL::getId("B_s0");
    IdaBs = EvtPDL::getId("anti-B_s0");
    IdRhop = EvtPDL::getId("rho+");
    IdRhom = EvtPDL::getId("rho-");
    IdRho0 = EvtPDL::getId("rho0");
    IdOmega = EvtPDL::getId("omega");
    IdKst0 = EvtPDL::getId("K*0");
    IdaKst0 = EvtPDL::getId("anti-K*0");
    IdKstp = EvtPDL::getId("K*+");
    IdKstm = EvtPDL::getId("K*-");
    IdK0 = EvtPDL::getId("K0");
    IdaK0 = EvtPDL::getId("anti-K0");
    IdKp = EvtPDL::getId("K+");
    IdKm = EvtPDL::getId("K-");
  }
 
  // First choose format of NP coefficients from the .DEC file
  // Cartesian(x,y)(0) or Polar(R,phase)(1)
  int n = getNArg();
  checkNDaug(3);
 
  //We expect the parent to be a scalar
  //and the daughters to be K(*) lepton+ lepton-
 
  checkSpinParent(EvtSpinType::SCALAR);
  checkSpinDaughter(1, EvtSpinType::DIRAC);
  checkSpinDaughter(2, EvtSpinType::DIRAC);
 
  EvtId mId = getDaug(0);
  if (mId != IdKst0 && mId != IdaKst0 && mId != IdKstp && mId != IdKstm &&
      mId != IdK0 && mId != IdaK0 && mId != IdKp && mId != IdKm) {
    EvtGenReport(EVTGEN_ERROR, "EvtGen")
        << "EvtbTosllNPR generator expected a K(*) 1st daughter, found: "
        << EvtPDL::name(getDaug(0)) << endl;
    EvtGenReport(EVTGEN_ERROR, "EvtGen")
        << "Will terminate execution!" << endl;
    ::abort();
  }
 
  auto getInteger = [this](int i) -> int {
    double a = getArg(i);
    if (a - static_cast<int>(a) != 0)
    {
      EvtGenReport(EVTGEN_ERROR, "EvtGen")
          << "Parameters is not integer in the BTOSLLNPR decay model: " << i
          << " " << a << endl;
      EvtGenReport(EVTGEN_ERROR, "EvtGen")
          << "Will terminate execution!" << endl;
      ::abort();
    }
    return static_cast<int>(a);
  };
 
  double phi = 0.0, n1 = 0.0, d1 = 1.0, eta = 0.0, theta_Jpsi = 0.0, theta_psi2S = 0.0, Nr = 1.0;
  if (n > 0) {      // parse arguments
    int i = 0, coordsyst = getInteger(i++);
    auto getcomplex = [this, &i, coordsyst]() -> EvtComplex {
      double a0 = getArg(i++);
      double a1 = getArg(i++);
      return (coordsyst) ? EvtComplex(a0 * cos(a1), a0 * sin(a1))
      : EvtComplex(a0, a1);
    };
    auto getreal = [this, &i]() -> double { return getArg(i++); };
    while (i < n) {
      int id = getInteger(i++);      // New Physics cooeficient Id
      if (id == 0)
        m_dc7 = getcomplex();    // delta C_7eff
      if (id == 1)
        m_dc9 = getcomplex();    // delta C_9eff
      if (id == 2)
        m_dc10 = getcomplex();    // delta C_10eff
      if (id == 3)
        m_c7p = getcomplex();    // C'_7eff -- right hand polarizations
      if (id == 4)
        m_c9p = getcomplex();    // C'_9eff -- right hand polarizations
      if (id == 5)
        m_c10p = getcomplex();    // c'_10eff -- right hand polarizations
      if (id == 6)
        m_cS = getcomplex();    // (C_S - C'_S) -- scalar right and left polarizations
      if (id == 7)
        m_cP = getcomplex();    // (C_P - C'_P) -- pseudo-scalar right and left polarizations
      if (id == 10)
        m_flag = getInteger(i++); // include resonances or not
      if (id == 11)
        phi = getreal(); // phase of the non-factorizable contribution
      if (id == 12)
        n1 = getreal(); // amplitude of the non-factorizable contribution
      if (id == 13)
        d1 = getreal(); // width of the non-factorizable contribution
      if (id == 14)
        eta = getreal(); // Eq. 47
      if (id == 15)
        theta_Jpsi = getreal(); // Eq. 47
      if (id == 16)
        theta_psi2S = getreal(); // Eq. 47
      if (id == 17)
        Nr = getreal(); // Eq. 47
    }
  }
  m_ls = new EvtLinSample;
 
  if (!m_flag) return;
  // BW_t jpsi(3.0969, 0.0926 / 1000, 0.0926 / 1000 * 0.877, 0.05971);
  // BW_t psi2s(3.6861, 0.294 / 1000, 0.294 / 1000 * 0.9785, 0.00793);
  // BW_t psi3770(3773.7 / 1000, 27.2 / 1000, 27.2 / 1000, 1e-5);
  // BW_t psi4040(4039 / 1000., 80 / 1000., 80 / 1000., 1.07e-5);
  // BW_t psi4160(4191 / 1000., 70 / 1000., 70 / 1000., 6.9e-6);
  // BW_t psi4230(4220 / 1000., 60 / 1000., 60 / 1000., 1e-5);
 
  TParticlePDG* p_jpsi = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("J/psi");
  // B2WARNING("extracting leptonic width"<<std::flush);
  // TDecayChannel *c_jpsi = p_jpsi->DecayChannel(1);
  // B2WARNING("Br = "<<c_jpsi->BranchingRatio());
  BW_t* jpsi = new BW_t(p_jpsi->Mass(), p_jpsi->Width(),  p_jpsi->Width() * 0.877, 0.05971);
  m_rs.push_back(jpsi);
 
  TParticlePDG* p_psi2s = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("psi(2S)");
  BW_t* psi2s = new BW_t(p_psi2s->Mass(), p_psi2s->Width(),  p_psi2s->Width() * 0.9785, 0.00793);
  m_rs.push_back(psi2s);
 
  TParticlePDG* p_psi3770 = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("psi(3770)");
  BW_t* psi3770 = new BW_t(p_psi3770->Mass(), p_psi3770->Width(), p_psi3770->Width(), 1e-5);
  m_rs.push_back(psi3770);
 
  TParticlePDG* p_psi4040 = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("psi(4040)");
  BW_t* psi4040 = new BW_t(p_psi4040->Mass(), p_psi4040->Width(), p_psi4040->Width(), 1.07e-5);
  m_rs.push_back(psi4040);
 
  TParticlePDG* p_psi4160 = Belle2::EvtGenDatabasePDG::Instance()->GetParticle("psi(4160)");
  BW_t* psi4160 = new BW_t(p_psi4160->Mass(), p_psi4160->Width(), p_psi4160->Width(), 6.9e-6);
  m_rs.push_back(psi4160);
 
  BW_t* psi4230 = new BW_t(4220 / 1000., 60 / 1000., 60 / 1000., 1e-5);
  m_rs.push_back(psi4230);
 
  std::vector<double> v = getgrid(m_rs);
  std::complex<double> eiphi = exp(1i * phi), pJpsi = Nr * exp(1i * theta_Jpsi), ppsi2S = Nr * exp(1i * theta_psi2S);
 
  const double C1 = -0.257, C2 = 1.009, C3 = -0.005, C5 = 0;
  double sc = 1.0 / (4.0 / 3.0 * C1 + C2 + 6.0 * C3 + 60.0 * C5);
  const double m_b = 4.8, m_c = 1.3;
  for (auto t : v) {
    // 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
    double re = -8.0 / 9.0 * log(m_c / m_b) - 4.0 / 9.0 +
                t / 3.0 * DRInt(t, real(pJpsi), real(ppsi2S), Nr, m_rs) - M_PI / 3 * uR(t, imag(pJpsi), imag(ppsi2S), 0, m_rs);
    double im =
      t / 3.0 * DRInt(t, imag(pJpsi), imag(ppsi2S), 0, m_rs) + M_PI / 3 * uR(t, real(pJpsi), real(ppsi2S), Nr, m_rs);
    std::complex<double> r(re, im);
    if (eta == 0.0) {
      double z = t - jpsi->Mv2;
      r += eiphi * (n1 * sc / (z * z + d1));    // bump under J/psi as a non-factorizable contribution
    }
    m_reso.push_back(std::make_pair(t, r));
  }
}

◆ initProbMax()

void initProbMax ( )

override

The method to evaluate the maximum amplitude.

Definition at line 496 of file EvtbTosllNPR.cc.

{
  EvtId pId = getParentId(), mId = getDaug(0), l1Id = getDaug(1), l2Id = getDaug(2);
 
  std::vector<double> v;
  std::vector<EvtPointf> pp;
 
  EvtScalarParticle* scalar_part;
  EvtParticle* root_part;
 
  scalar_part = new EvtScalarParticle;
 
  //includge to avoid generating random numbers!
  scalar_part->noLifeTime();
 
  EvtVector4R p_init;
 
  p_init.set(EvtPDL::getMass(pId), 0.0, 0.0, 0.0);
  scalar_part->init(pId, p_init);
  root_part = (EvtParticle*)scalar_part;
  root_part->setDiagonalSpinDensity();
 
  EvtParticle* daughter, *lep1, *lep2;
 
  EvtAmp amp;
 
  EvtId listdaug[3];
  listdaug[0] = mId;
  listdaug[1] = l1Id;
  listdaug[2] = l2Id;
 
  amp.init(pId, 3, listdaug);
 
  root_part->makeDaughters(3, listdaug);
  daughter = root_part->getDaug(0);
  lep1 = root_part->getDaug(1);
  lep2 = root_part->getDaug(2);
 
  //cludge to avoid generating random numbers!
  daughter->noLifeTime();
  lep1->noLifeTime();
  lep2->noLifeTime();
 
  //Initial particle is unpolarized, well it is a scalar so it is trivial
  EvtSpinDensity rho;
  rho.setDiag(root_part->getSpinStates());
 
  double mass[3];
 
  double m = root_part->mass();
 
  EvtVector4R p4meson, p4lepton1, p4lepton2, p4w;
 
  double maxprobfound = 0.0;
 
  mass[1] = EvtPDL::getMeanMass(l1Id);
  mass[2] = EvtPDL::getMeanMass(l2Id);
  std::vector<double> mH;
  mH.push_back(EvtPDL::getMeanMass(mId));
  //if the particle is narrow dont bother with changing the mass.
  double g = EvtPDL::getWidth(mId);
  if (g > 0) {
    mH.push_back(EvtPDL::getMinMass(mId));
    mH.push_back(
      std::min(EvtPDL::getMaxMass(mId), m - mass[1] - mass[2]));
    mH.push_back(mH.front() - g);
    mH.push_back(mH.front() + g);
  }
 
  double q2min = (mass[1] + mass[2]) * (mass[1] + mass[2]);
 
  double m0 = EvtPDL::getMinMass(mId);
  double q2max = (m - m0) * (m - m0);
  v = getgrid(m_rs);
  m0 = mH[0];
  //loop over q2
  for (double q2 : v) {
    // want to avoid picking up the tail of the photon propagator
    if (!(q2min <= q2 && q2 < q2max))
      continue;
    double Erho = (m * m + m0 * m0 - q2) / (2.0 * m);
    if (Erho < m0) {
      m0 = EvtPDL::getMinMass(mId);
      Erho = (m * m + m0 * m0 - q2) / (2.0 * m);
    }
    double Prho = sqrt((Erho - m0) * (Erho + m0));
 
    p4meson.set(Erho, 0.0, 0.0, -Prho);
    p4w.set(m - Erho, 0.0, 0.0, Prho);
 
    //This is in the W rest frame
    double Elep = sqrt(q2) * 0.5,
           Plep = sqrt((Elep - mass[1]) * (Elep + mass[1]));
 
    const int nj = 3 + 2 + 4 + 8 + 32;
    double cmin = -1, cmax = 1, dc = (cmax - cmin) / (nj - 1);
    double maxprob = 0;
    for (int j = 0; j < nj; j++) {
      double c = cmin + dc * j;
 
      //These are in the W rest frame. Need to boost out into the B frame.
      double Py = Plep * sqrt(1.0 - c * c), Pz = Plep * c;
      p4lepton1.set(Elep, 0.0, Py, Pz);
      p4lepton2.set(Elep, 0.0, -Py, -Pz);
 
      p4lepton1 = boostTo(p4lepton1, p4w);
      p4lepton2 = boostTo(p4lepton2, p4w);
 
      //Now initialize the daughters...
      daughter->init(mId, p4meson);
      lep1->init(l1Id, p4lepton1);
      lep2->init(l2Id, p4lepton2);
      CalcAmp(root_part, amp);
      double prob = rho.normalizedProb(amp.getSpinDensity());
      maxprob = std::max(maxprob, prob);
    }
    pp.push_back({ (float)q2, (float)maxprob });
    maxprobfound = std::max(maxprobfound, maxprob);
  }
  root_part->deleteTree();
 
  m_ls->init(pp);
  maxprobfound = 4;
  setProbMax(maxprobfound);
}

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:

generators/evtgen/models/include/EvtbTosllNPR.h
generators/evtgen/models/src/EvtbTosllNPR.cc

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ ~EvtbTosllNPR()

Member Function Documentation

◆ CalcAmp()

◆ CalcSAmp()

◆ CalcVAmp()

◆ clone()

◆ decay()

◆ getName()

◆ init()

◆ initProbMax()

Member Data Documentation

◆ m_c10p

◆ m_c7p

◆ m_c9p

◆ m_cP

◆ m_cS

◆ m_dc10

◆ m_dc7

◆ m_dc9

◆ m_flag

◆ m_ls

◆ m_reso

◆ m_rs