diff --git a/EvtGenModels/EvtBLLNuL.hh b/EvtGenModels/EvtBLLNuL.hh index ee377d3..b3ff38d 100644 --- a/EvtGenModels/EvtBLLNuL.hh +++ b/EvtGenModels/EvtBLLNuL.hh @@ -1,49 +1,49 @@ //-------------------------------------------------------------------------- // // Copyright Information: See EvtGen/COPYRIGHT // // Module: EvtBLLNuL.hh // // Description: The header file for the model "BLLNUL" which simulates // the rare four-leptonic B-decays // B^-(p) -> ell^+(k_1) ell^-(k_2) neu (k_3) ell^-(k_4) // // Modification history: // // Anna Danilina (anna.danilina@cern.ch) and // Nikolai Nikitin (Nikolai.Nikitine@cern.ch) Nov 2018 Module created // //------------------------------------------------------------------------ #ifndef EVTBLLNUL_HH #define EVTBLLNUL_HH #include "EvtGenBase/EvtDecayAmp.hh" +#include "EvtGenModels/EvtBLLNuLAmp.hh" #include class EvtParticle; class EvtbTosllMSFF; // Form factor class -class EvtBLLNuLAmp; // Amplitude class class EvtBLLNuL: public EvtDecayAmp { public: EvtBLLNuL(); virtual ~EvtBLLNuL(); virtual std::string getName(); virtual EvtDecayBase* clone(); virtual void init(); virtual void initProbMax(); virtual void decay(EvtParticle *p); private: - EvtBLLNuLAmp* calcAmp_; + EvtBLLNuLAmp calcAmp_; }; #endif diff --git a/EvtGenModels/EvtBLLNuLAmp.hh b/EvtGenModels/EvtBLLNuLAmp.hh index 4ee80e1..2b5d70e 100644 --- a/EvtGenModels/EvtBLLNuLAmp.hh +++ b/EvtGenModels/EvtBLLNuLAmp.hh @@ -1,132 +1,119 @@ //-------------------------------------------------------------------------- // // // Module: EvtB2MuMuENuAmp.hh // // Description: Header file for the amplitude calculation for the "BLLNUL" // model which generates rare four-leptonic B-decays // B^-(p) -> ell^+(k_1) ell^-(k_2) neu (k_3) ell^-(k_4) // // Modification history: // // Anna Danilina (anna.danilina@cern.ch) and // Nikolai Nikitin (Nikolai.Nikitine@cern.ch) Nov 2018 Module created // John B Code optimisations // //------------------------------------------------------------------------ #ifndef EVTBLLNUL_AMP_HH #define EVTBLLNUL_AMP_HH #include "EvtGenBase/EvtAmp.hh" #include "EvtGenBase/EvtComplex.hh" #include "EvtGenBase/EvtId.hh" #include "EvtGenBase/EvtTensor4C.hh" #include "EvtGenBase/EvtVector4R.hh" #include class EvtParticle; class EvtBLLNuLAmp { public: + EvtBLLNuLAmp(double Vub = 4.09e-3); EvtBLLNuLAmp(double qSqMin, double kSqMin, bool symmetry, double Vub = 4.09e-3); virtual ~EvtBLLNuLAmp(); void CalcAmp(EvtParticle *parent, EvtAmp& amp) const; + void setParameters(double qSqMin, double kSqMin, bool symmetry); // Resonance poles class ResPole { public: ResPole(double mass, double width, double coupling); virtual ~ResPole() {;} EvtComplex propagator(double qSq, int numForm = 0) const; double getMass() const {return m0_;} double getMassSq() const {return m0Sq_;} double getWidth() const {return w0_;} double getCoupling() const {return c_;} private: double m0_; // pole mass double m0Sq_; double w0_; // width double c_; // coupling constant EvtComplex I_; EvtComplex Imw_; }; - // Kinematics - class KinInfo { - public: - KinInfo(const EvtVector4R& q, const EvtVector4R& k, - double qSq, double kSq, double MB, - int sign); - virtual ~KinInfo() {;} - - EvtVector4R getQ() const {return q_;} - EvtVector4R getK() const {return k_;} - double getQSq() const {return qSq_;} - double getKSq() const {return kSq_;} - double getMB() const {return MB_;} - int getSign() const {return sign_;} - - private: - EvtVector4R q_; - EvtVector4R k_; - double qSq_; - double kSq_; - double MB_; - int sign_; - }; - protected: - EvtTensor4C getHadronTensor(const EvtBLLNuLAmp::KinInfo& info) const; + EvtTensor4C getHadronTensor(const EvtVector4R& q, const EvtVector4R& k, + const double qSq, const double kSq, const double MB, + const int sign) const; std::vector getVMDTerms(double qSq, double kSq, double MB) const; EvtComplex getBStarTerm(double qSq, double kSq, double MB) const; double FF_B2Bstar(double qSq) const; double FF_V(double kSq) const; double FF_A1(double kSq) const; double FF_A2(double kSq) const; private: // Kinematic cut-offs double qSqMin_; double kSqMin_; // If we have identical charged lepton flavours bool symmetry_; // B+, B- Ids EvtId BpId_, BnId_; // Form factor constants double coupling_, sqrt2_; double fBu_; // Resonance poles EvtBLLNuLAmp::ResPole* Bstar_; EvtBLLNuLAmp::ResPole* Upsilon_; std::vector resPoles_; int nPoles_; // Complex number constants EvtComplex zero_, unitI_; }; +inline void EvtBLLNuLAmp::setParameters(double qSqMin, double kSqMin, bool symmetry) +{ + qSqMin_ = qSqMin; + kSqMin_ = kSqMin; + symmetry_ = symmetry; +} + #endif diff --git a/src/EvtGenModels/EvtBLLNuL.cpp b/src/EvtGenModels/EvtBLLNuL.cpp index 9f49ef5..0608ab9 100644 --- a/src/EvtGenModels/EvtBLLNuL.cpp +++ b/src/EvtGenModels/EvtBLLNuL.cpp @@ -1,158 +1,156 @@ //-------------------------------------------------------------------------- // // Copyright Information: See EvtGen/COPYRIGHT // // Module: EvtBLLNuL.cpp // // Description: The main file for the model "BLLNUL" which simulates // the rare four-leptonic B-decays // B^-(p) -> ell^+(k_1) ell^-(k_2) neu (k_3) ell^-(k_4) // // Modification history: // // Anna Danilina (anna.danilina@cern.ch) and // Nikolai Nikitin (Nikolai.Nikitine@cern.ch) Nov 2018 Module created // John B Code optimisations // //------------------------------------------------------------------------ // #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtSpinType.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtIdSet.hh" #include "EvtGenModels/EvtBLLNuL.hh" -#include "EvtGenModels/EvtBLLNuLAmp.hh" EvtBLLNuL::EvtBLLNuL() : - calcAmp_(0) + calcAmp_() { } EvtBLLNuL::~EvtBLLNuL() { - delete calcAmp_; } std::string EvtBLLNuL::getName() { // The model name return "BLLNUL"; } EvtDecayBase* EvtBLLNuL::clone() { return new EvtBLLNuL(); } void EvtBLLNuL::init() { // check that there are 4 daughters checkNDaug(4); // We expect that the parent to be a scalar (B meson) and // the daughters to be ell+ (k1), ell- (k2), neutrino (k3) // and the last lepton ell- (k4) // Check spin types checkSpinParent(EvtSpinType::SCALAR); checkSpinDaughter(0, EvtSpinType::DIRAC); // ell+(k_1) checkSpinDaughter(1, EvtSpinType::DIRAC); // ell-(k_2) checkSpinDaughter(2, EvtSpinType::NEUTRINO); // neu (k_3) checkSpinDaughter(3, EvtSpinType::DIRAC); // ell-(k_4) // Check that we have a charged B parent static EvtIdSet BMesons("B-", "B+"); if (!BMesons.contains(getParentId())) { EvtGenReport(EVTGEN_ERROR, "EvtBLLNuL") << "Expecting the parent to be a charged B. Found PDG = " << EvtPDL::getStdHep(getParentId()) << std::endl; ::abort(); } // Make sure the first two leptons are charged conjugates of each other int id1 = EvtPDL::getStdHep(getDaug(0)); int id2 = EvtPDL::getStdHep(getDaug(1)); if (id1 != -id2) { EvtGenReport(EVTGEN_ERROR, "EvtBLLNuL") << "Expecting the first 2 leptons, with PDG codes " << id1 << " and " << id2 << ", to be charged conjugates of each other" << std::endl; ::abort(); } // Check that the last lepton has the same charge as the B parent int q3 = EvtPDL::chg3(getDaug(3))/3; int qB = EvtPDL::chg3(getParentId())/3; if (q3 != qB) { EvtGenReport(EVTGEN_ERROR, "EvtBLLNuL") << "The 3rd lepton charge " << q3 << " does not match the B charge " << qB << std::endl; ::abort(); } // Also check that the 2nd lepton has the same charge as the 3rd one int q2 = EvtPDL::chg3(getDaug(1))/3; if (q2 != q3) { EvtGenReport(EVTGEN_ERROR, "EvtBLLNuL") << "The 2nd lepton charge " << q2 << " does not match the 3rd lepton charge " << q3 << std::endl; ::abort(); } // Identify if the decay has 3 charged leptons with the same flavour. // If so, then we need to include amplitude terms where the 2nd and 3rd // same-sign leptons are interchanged: k2 <-> k4 bool symmetry(false); int id3 = EvtPDL::getStdHep(getDaug(3)); if (abs(id1) == abs(id2) && abs(id1) == abs(id3)) { symmetry = true; } // Specify the qSq minimum cut-off as 4*(muon mass)^2 = 0.044655 and the // kSq minimum cut off as 4*(electron mass)^2 = 1.044e-6 double muMass = EvtPDL::getMeanMass(EvtPDL::getId("mu+")); double eMass = EvtPDL::getMeanMass(EvtPDL::getId("e+")); double qSqMin = 4.0*muMass*muMass; double kSqMin = 4.0*eMass*eMass; // Optionally set these cut-offs using two decay file parameters. We may // have a 3rd parameter (max prob), so check for at least 2 parameters if (getNArg() >= 2) { qSqMin = getArg(0); kSqMin = getArg(1); } - // Define the amplitude calculation pointer with the qSq and kSq cut-offs, - // also specifying if the decay mode has flavour symmetry - calcAmp_ = new EvtBLLNuLAmp(qSqMin, kSqMin, symmetry); + // Define the amplitude qSq and kSq cut-offs, also + // specifying if the decay mode has flavour symmetry + calcAmp_.setParameters(qSqMin, kSqMin, symmetry); } void EvtBLLNuL::initProbMax() { // Set the maximum probability of the decay double maxProb(3.2); // Optional 3rd decay file parameter, e.g. if qSq and/or kSq min have changed. // Note that both qSq and kSq parameters must still be specified in the decay // file to ensure that the maximum probability value is the 3rd parameter! if (getNArg() == 3) { maxProb = getArg(2); } setProbMax(maxProb); } void EvtBLLNuL::decay(EvtParticle *p) { p->initializePhaseSpace(getNDaug(), getDaugs()); - calcAmp_->CalcAmp(p, _amp2); + calcAmp_.CalcAmp(p, _amp2); } diff --git a/src/EvtGenModels/EvtBLLNuLAmp.cpp b/src/EvtGenModels/EvtBLLNuLAmp.cpp index cf03c88..6ba313b 100644 --- a/src/EvtGenModels/EvtBLLNuLAmp.cpp +++ b/src/EvtGenModels/EvtBLLNuLAmp.cpp @@ -1,426 +1,440 @@ //-------------------------------------------------------------------------- // // Module: EvtBLLNuLAmp.cpp // // Description: Amplitude calculation class for the "BLLNUL" model of // rare four-leptonic B-decays: // B^-(p) -> ell^+(k_1) ell^-(k_2) neu (k_3) ell^-(k_4) // // Modification history: // // Anna Danilina (anna.danilina@cern.ch) and // Nikolai Nikitin (Nikolai.Nikitine@cern.ch) Nov 2018 Module created // John B Code optimisations // //----------------------------------------------------------------------------------------- // #include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtConst.hh" #include "EvtGenBase/EvtDiracSpinor.hh" #include "EvtGenBase/EvtIdSet.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtVector4C.hh" #include "EvtGenModels/EvtBLLNuLAmp.hh" #include +EvtBLLNuLAmp::EvtBLLNuLAmp(double Vub) : + qSqMin_(0.0), + kSqMin_(0.0), + symmetry_(false), + BpId_(EvtPDL::getId("B+")), + BnId_(EvtPDL::getId("B-")), + coupling_(0.0), + sqrt2_(sqrt(2.0)), + fBu_(0.191), // leptonic constant (GeV) + Bstar_(new EvtBLLNuLAmp::ResPole(5.32, 0.00658, 0.183/3.0)), + Upsilon_(new EvtBLLNuLAmp::ResPole(9.64, 0.0, 0.0)), + resPoles_(), + nPoles_(0), + zero_(EvtComplex(0.0, 0.0)), + unitI_(EvtComplex(0.0, 1.0)) +{ + double GF = 1.166371e-5; // GeV^{-2} + double alphaEM = 1.0/137.0; + + // Normalisation constant, multiplied by 1e4 to increase probability scale + coupling_ = 400.0*GF*EvtConst::pi*alphaEM*Vub*1e4/sqrt2_; + + // Define VMD resonance poles using PDG 2016 values with constants from + // D.Melikhov, N.Nikitin and K.Toms, Phys. Atom. Nucl. 68, 1842 (2005) + EvtBLLNuLAmp::ResPole* rho = new EvtBLLNuLAmp::ResPole(0.77526, 0.1491, 1.0/5.04); + resPoles_.push_back(rho); + + EvtBLLNuLAmp::ResPole* omega = new EvtBLLNuLAmp::ResPole(0.78265, 0.00849, 1.0/17.1); + resPoles_.push_back(omega); + + nPoles_ = resPoles_.size(); +} + EvtBLLNuLAmp::EvtBLLNuLAmp(double qSqMin, double kSqMin, bool symmetry, double Vub) : qSqMin_(qSqMin), kSqMin_(kSqMin), symmetry_(symmetry), BpId_(EvtPDL::getId("B+")), BnId_(EvtPDL::getId("B-")), coupling_(0.0), sqrt2_(sqrt(2.0)), fBu_(0.191), // leptonic constant (GeV) Bstar_(new EvtBLLNuLAmp::ResPole(5.32, 0.00658, 0.183/3.0)), Upsilon_(new EvtBLLNuLAmp::ResPole(9.64, 0.0, 0.0)), resPoles_(), nPoles_(0), zero_(EvtComplex(0.0, 0.0)), unitI_(EvtComplex(0.0, 1.0)) { double GF = 1.166371e-5; // GeV^{-2} double alphaEM = 1.0/137.0; - double MyPi = EvtConst::pi; + // Normalisation constant, multiplied by 1e4 to increase probability scale - coupling_ = 400.0*GF*MyPi*alphaEM*Vub*1e4/sqrt2_; + coupling_ = 400.0*GF*EvtConst::pi*alphaEM*Vub*1e4/sqrt2_; // Define VMD resonance poles using PDG 2016 values with constants from // D.Melikhov, N.Nikitin and K.Toms, Phys. Atom. Nucl. 68, 1842 (2005) EvtBLLNuLAmp::ResPole* rho = new EvtBLLNuLAmp::ResPole(0.77526, 0.1491, 1.0/5.04); resPoles_.push_back(rho); EvtBLLNuLAmp::ResPole* omega = new EvtBLLNuLAmp::ResPole(0.78265, 0.00849, 1.0/17.1); resPoles_.push_back(omega); nPoles_ = resPoles_.size(); } EvtBLLNuLAmp::~EvtBLLNuLAmp() { delete Bstar_; delete Upsilon_; // Delete VMD poles std::vector::iterator iter; for (iter = resPoles_.begin(); iter != resPoles_.end(); ++iter) { delete *iter; } resPoles_.clear(); } // Storing resonance pole information EvtBLLNuLAmp::ResPole::ResPole(double mass, double width, double coupling) : m0_(mass), m0Sq_(mass*mass), w0_(width), c_(coupling), I_(EvtComplex(0.0, 1.0)), Imw_(I_*mass*width) { } EvtComplex EvtBLLNuLAmp::ResPole::propagator(double qSq, int numForm) const { // Numerator term: mass-squared (default) or mass double num(m0Sq_); if (numForm == 1) {num = m0_;} EvtComplex result = num*c_/((qSq - m0Sq_) + Imw_); return result; } -// Store kinematic information that can be passed around functions -EvtBLLNuLAmp::KinInfo::KinInfo(const EvtVector4R& q, const EvtVector4R& k, - double qSq, double kSq, double MB, int sign) : - q_(q), - k_(k), - qSq_(qSq), - kSq_(kSq), - MB_(MB), - sign_(sign) -{ -} - // Amplitude calculation void EvtBLLNuLAmp::CalcAmp(EvtParticle *parent, EvtAmp& amp) const { // Check for 4 daughters and an existing parent if (!parent || parent->getNDaug() != 4) {return;} // The first two charged leptons. The 2nd one will have // the same charge as the 3rd charged lepton EvtParticle* lepA = parent->getDaug(0); EvtParticle* lepB = parent->getDaug(1); // The neutrino EvtParticle* neu = parent->getDaug(2); // The third charged lepton EvtParticle* lepC = parent->getDaug(3); // Kinematics double MB = parent->mass(); // B-meson mass, GeV // 4-momenta of the leptons in the B rest frame. The daughters will already // be in the correct order since this check is done in EvtBLLNuL::init() // when initialising the model using the decay file EvtVector4R p1 = lepA->getP4(); EvtVector4R p2 = lepB->getP4(); EvtVector4R p3 = neu->getP4(); EvtVector4R p4 = lepC->getP4(); // 4-momenta sums EvtVector4R q12 = p1 + p2; EvtVector4R k34 = p3 + p4; // Mandelstam variables: q^2 and k^2 double q12Sq = q12.mass2(); double k34Sq = k34.mass2(); // Check if we are above mass thresholds bool threshold(true); if (q12Sq < qSqMin_ || k34Sq < kSqMin_) { threshold = false; } // For the symmetric terms when we exchange the // 2nd and 3rd charged leptons: p2 <-> p4 EvtVector4R q14, k23; double q14Sq(0.0), k23Sq(0.0); if (symmetry_) { q14 = p1 + p4; k23 = p2 + p3; q14Sq = q14.mass2(); k23Sq = k23.mass2(); if (q14Sq < qSqMin_ || k23Sq < kSqMin_) { threshold = false; } } // B meson id EvtId parId = parent->getId(); // B+ or B- decays int sign(1); if (parId == BnId_) {sign = -1;} // Hadronic tensors - EvtBLLNuLAmp::KinInfo infoA(q12, k34, q12Sq, k34Sq, MB, sign); - EvtTensor4C THadronA = getHadronTensor(infoA); + EvtTensor4C THadronA = getHadronTensor(q12, k34, q12Sq, k34Sq, MB, sign); // When we need to include the symmetric terms EvtTensor4C THadronB; if (symmetry_) { - EvtBLLNuLAmp::KinInfo infoB(q14, k23, q14Sq, k23Sq, MB, sign); - THadronB = getHadronTensor(infoB); + THadronB = getHadronTensor(q14, k23, q14Sq, k23Sq, MB, sign); } // Leptonic currents: A for normal terms, B for symmetric terms EvtVector4C L1A, L2A, L1B, L2B; int leptonSpins[4]; // array for saving the leptonic spin configuration // Loop over lepton spin states for (int i2 = 0; i2 < 2; i2++) { leptonSpins[0] = i2; for (int i1 = 0; i1 < 2; i1++) { leptonSpins[1] = i1; if (sign == -1) { // B- currents // L2^{\nu} = \bar mu(k_2) \gamma^{\nu} mu(- k_1) L2A = EvtLeptonVCurrent(lepB->spParent(i2), lepA->spParent(i1)); if (symmetry_) { // Swapping the 2nd and 3rd charged leptons L1B = EvtLeptonVACurrent(lepB->spParent(i2), neu->spParentNeutrino()); } } else { // B+ currents // L2^{\nu} = \bar mu(k_1) \gamma^{\nu} mu(- k_2) L2A = EvtLeptonVCurrent(lepA->spParent(i1), lepB->spParent(i2)); if (symmetry_) { // Swapping the 2nd and 3rd charged leptons L1B = EvtLeptonVACurrent(neu->spParentNeutrino(), lepB->spParent(i2)); } } // Production: Tfi^{\mu} = THadron^{\mu \nu} L_{2 \nu} EvtVector4C THL2A = THadronA.cont2(L2A); for (int i4 = 0; i4 < 2; i4++) { leptonSpins[2] = i4; leptonSpins[3] = 0; // neutrino handedness if (sign == -1) { // B- currents // L1^{\mu} = \bar e(k_4) \gamma^{\mu} (1 - \gamma^5) nu_e(- k_3) L1A = EvtLeptonVACurrent(lepC->spParent(i4), neu->spParentNeutrino()); if (symmetry_) { // Swapping the 2nd and 3rd charged leptons L2B = EvtLeptonVCurrent(lepC->spParent(i4), lepA->spParent(i1)); } } else { // B+ currents // L1^{\mu} = \bar nu_e(k_3) \gamma^{\mu} (1 - \gamma^5) e(- k_4) L1A = EvtLeptonVACurrent(neu->spParentNeutrino(), lepC->spParent(i4)); if (symmetry_) { // Swapping the 2nd and 3rd charged leptons L2B = EvtLeptonVCurrent(lepA->spParent(i1), lepC->spParent(i4)); } } if (threshold == false) { // Below kinematic thresholds amp.vertex(leptonSpins, zero_); } else { // Decay amplitude calculation: L_1^{\mu} Tfi_{\mu} EvtComplex decAmp = L1A*THL2A; // If we also need to swap the 2nd and 3rd charged leptons if (symmetry_) { // Hadronic current production term. L2B depends on i4 so we need // it here instead of inside the i2 loop as was the case for THL2A EvtVector4C THL2B = THadronB.cont2(L2B); // The symmetric amplitude EvtComplex ampB = L1B*THL2B; // Subtract this from the total amplitude decAmp -= ampB; } amp.vertex(leptonSpins, decAmp); } } // i4 loop } // i1 loop } // i2 loop } -EvtTensor4C EvtBLLNuLAmp::getHadronTensor(const EvtBLLNuLAmp::KinInfo& info) const +EvtTensor4C EvtBLLNuLAmp::getHadronTensor(const EvtVector4R& q, const EvtVector4R& k, + const double qSq, const double kSq, const double MB, + const int sign) const { - // Hadronic tensor calculation. - // First retrieve kinematic variables - EvtVector4R q = info.getQ(); - EvtVector4R k = info.getK(); - double qSq = info.getQSq(); - double kSq = info.getKSq(); - double MB = info.getMB(); - int sign = info.getSign(); + // Hadronic tensor calculation EvtTensor4C epskq = dual(EvtGenFunctions::directProd(k, q)); EvtTensor4C qk = EvtGenFunctions::directProd(q, k); EvtComplex BstarAmp = getBStarTerm(qSq, kSq, MB); std::vector VMDAmps = getVMDTerms(qSq, kSq, MB); EvtComplex FF_ekq = BstarAmp + VMDAmps[0]; EvtComplex FF_g = VMDAmps[1] - fBu_; EvtComplex FF_qk = VMDAmps[2]; // Full hadronic tensor EvtTensor4C THadron = sign*2.0*FF_ekq*epskq + unitI_*(2.0*FF_qk*qk - FF_g*EvtTensor4C::g()); // Kinematic cuts double coeffcut(0.0); if (qSq > qSqMin_ && kSq > kSqMin_) { coeffcut = 1.0/qSq; } // Normalisation constant THadron *= coeffcut*coupling_; return THadron; } std::vector EvtBLLNuLAmp::getVMDTerms(double qSq, double kSq, double MB) const { // Find the 3 VMD form factors: epsilon*k*q, g(uv) and q*k terms EvtComplex VMD1(0.0, 0.0), VMD2(0.0, 0.0), VMD3(0.0, 0.0); // Loop over the VMD poles for (int iPole = 0; iPole < nPoles_; iPole++) { EvtBLLNuLAmp::ResPole* pole = resPoles_[iPole]; // Propagator term, common for all factors EvtComplex prop = pole->propagator(qSq); double mSum = MB + pole->getMass(); VMD1 += prop/mSum; VMD2 += mSum*prop; } // Third pole summation term is the same as the first one VMD3 = VMD1; // Multiply by couplings for the given kSq VMD1 *= FF_V(kSq); VMD2 *= FF_A1(kSq); VMD3 *= FF_A2(kSq); // Return the factors as a vector std::vector factors; factors.push_back(VMD1); factors.push_back(VMD2); factors.push_back(VMD3); return factors; } EvtComplex EvtBLLNuLAmp::getBStarTerm(double qSq, double kSq, double MB) const { EvtComplex amplitude = Bstar_->propagator(kSq, 1)*FF_B2Bstar(qSq)/(MB + Bstar_->getMass()); return amplitude; } double EvtBLLNuLAmp::FF_B2Bstar(double qSq) const { // Electromagnetic FF for B -> B* transition, when gamma is emitted from the b quark // D.Melikhov, private communication double y = qSq/Upsilon_->getMassSq(); double denom = (1.0 - y)*(1.0 - 0.81*y); double V(0.0); if (fabs(denom) > 1e-10) { V = 1.044/denom; } return V; } double EvtBLLNuLAmp::FF_V(double kSq) const { // D. Melikhov and B. Stech, PRD 62, 014006 (2000) Table XV double y = kSq/Bstar_->getMassSq(); double denom = sqrt2_*(1.0 - y)*(1.0 - 0.59*y); double V(0.0); if (fabs(denom) > 1e-10) { V = 0.31/denom; } return V; } double EvtBLLNuLAmp::FF_A1(double kSq) const { // D. Melikhov and B. Stech, PRD 62, 014006 (2000) Table XV double y = kSq/Bstar_->getMassSq(); double denom = ((0.1*y - 0.73)*y + 1.0)*sqrt2_; double A1(0.0); if (fabs(denom) > 1e-10) { A1 = 0.26/denom; } return A1; } double EvtBLLNuLAmp::FF_A2(double kSq) const { // D. Melikhov and B. Stech, PRD 62, 014006 (2000) Table XV double y = kSq/Bstar_->getMassSq(); double denom = ((0.5*y - 1.4)*y + 1.0)*sqrt2_; double A2(0.0); if (fabs(denom) > 1e-10) { A2 = 0.24/denom; } return A2; }