Page MenuHomeHEPForge
Differential D87 Diff 389 src/EvtGenModels/EvtbTosllVectorAmpNP.cpp

Changeset View

Standalone View

View Options

src/EvtGenModels/EvtbTosllVectorAmpNP.cpp

  • This file was added.
/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
* *
* This file is part of EvtGen. *
* *
* EvtGen is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* EvtGen is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/
#include "EvtGenModels/EvtbTosllVectorAmpNP.hh"
#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenModels/EvtbTosllAmp.hh"
#include "EvtGenModels/EvtbTosllFF.hh"
#include <algorithm>
#include <ccomplex>
static std::complex<double> h( double q2, double mq )
{
const double mu = 4.8; // mu = m_b = 4.8 GeV
if ( mq == 0.0 )
return std::complex<double>( 8. / 27 + 4. / 9 * log( mu * mu / q2 ),
4. / 9 * M_PI );
double z = 4 * mq * mq / q2;
std::complex<double> t;
if ( z > 1 ) {
t = atan( 1 / sqrt( z - 1 ) );
} else {
t = std::complex<double>( log( ( 1 + sqrt( 1 - z ) ) / sqrt( z ) ),
-M_PI / 2 );
}
return -4 / 9. * ( 2 * log( mq / mu ) - 2 / 3. - z ) -
( 4 / 9. * ( 2 + z ) * sqrt( fabs( z - 1 ) ) ) * t;
}
static EvtComplex C9( double q2,
std::complex<double> hReso = std::complex<double>( 0, 0 ) )
{
double m_c = 1.3, m_b = 4.8; // quark masses
double C1 = -0.257, C2 = 1.009, C3 = -0.005, C4 = -0.078, C5 = 0, C6 = 0.001;
std::complex<double> Y = ( h( q2, m_c ) + hReso ) *
( 4 / 3. * C1 + C2 + 6 * C3 + 60 * C5 );
Y -= 0.5 * h( q2, m_b ) * ( 7 * C3 + 4 / 3. * C4 + 76 * C5 + 64 / 3. * C6 );
Y -= 0.5 * h( q2, 0.0 ) * ( C3 + 4 / 3. * C4 + 16 * C5 + 64 / 3. * C6 );
Y += 4 / 3. * C3 + 64 / 9. * C5 + 64 / 27. * C6;
return EvtComplex( 4.211 + real( Y ), imag( Y ) );
}
static std::complex<double> interpol(
const std::vector<std::pair<double, std::complex<double>>>& v, double s )
{
if ( !v.size() )
return 0;
int j = std::lower_bound( v.begin(), v.end(), s,
[]( const std::pair<double, std::complex<double>>& a,
float b ) { return a.first < b; } ) -
v.begin();
double dx = v[j].first - v[j - 1].first;
auto dy = v[j].second - v[j - 1].second;
return v[j - 1].second + ( s - v[j - 1].first ) * ( dy / dx );
}
static EvtId IdB0, IdaB0, IdBp, IdBm, IdBs, IdaBs, IdRhop, IdRhom, IdRho0,
IdOmega, IdKst0, IdaKst0, IdKstp, IdKstm, IdK0, IdaK0, IdKp, IdKm;
static bool cafirst = true;
void EvtbTosllVectorAmpNP::CalcAmp( EvtParticle* parent, EvtAmp& amp,
EvtbTosllFF* formFactors )
{
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-" );
}
EvtId dId = parent->getDaug( 0 )->getId();
if ( dId == IdKst0 || dId == IdaKst0 || dId == IdKstp || dId == IdKstm ) {
CalcVAmp( parent, amp, formFactors );
}
if ( dId == IdK0 || dId == IdaK0 || dId == IdKp || dId == IdKm ) {
CalcSAmp( parent, amp, formFactors );
}
}
void EvtbTosllVectorAmpNP::CalcVAmp( EvtParticle* parent, EvtAmp& amp,
EvtbTosllFF* formFactors )
{
// Add the lepton and neutrino 4 momenta to find q2
EvtId pId = parent->getId();
EvtId dId = parent->getDaug( 0 )->getId();
EvtVector4R q = parent->getDaug( 1 )->getP4() + parent->getDaug( 2 )->getP4();
double q2 = q.mass2();
double mesonmass = parent->getDaug( 0 )->mass();
double a1, a2, a0, v, t1, t2, t3; // form factors
formFactors->getVectorFF( pId, dId, q2, 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 = ( 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;
EvtTensor4C tt0( EvtGenFunctions::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;
T1.addScaled( -b, gmn );
T1.addScaled( c, tt1 );
T1.addScaled( d, tt2 );
T2 *= e;
T2.addScaled( -f, gmn );
T2.addScaled( g, tt1 );
T2.addScaled( 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;
T1.addScaled( -b, gmn );
T1.addScaled( c, tt1 );
T1.addScaled( d, tt2 );
T2 *= -e;
T2.addScaled( -f, gmn );
T2.addScaled( g, tt1 );
T2.addScaled( 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 );
}
}
void getScalarFF( double q2, double& fp, double& f0, double& fT )
{
// The form factors are taken from:
// B -> K and D -> K form factors from fully relativistic lattice QCD
// W. G. Parrott, C. Bouchard, C. T. H. Davies
// https://arxiv.org/abs/2207.12468
static const double MK = 0.495644, MB = 5.279495108051249,
MBs0 = 5.729495108051249, MBsstar = 5.415766151925566,
logsB = 1.3036556717286227;
static const int N = 3;
static const double a0B[] = { 0.2546162729845652, 0.211016713841977,
0.02690776720598433, 0.2546162729845652,
-0.7084710010870424, 0.3096901516968882,
0.2549078414069112, -0.6549384905373116,
0.36265904141973127 },
*apB = a0B + N, *aTB = apB + N;
double pole0 = 1 / ( 1 - q2 / ( MBs0 * MBs0 ) );
double polep = 1 / ( 1 - q2 / ( MBsstar * MBsstar ) );
double B = MB + MK, A = sqrt( B * B - q2 ), z = ( A - B ) / ( A + B ),
zN = z * z * z / N, zn = z;
f0 = a0B[0];
fp = apB[0];
fT = aTB[0];
for ( int i = 1; i < N; i++ ) {
f0 += a0B[i] * zn;
fp += apB[i] * zn;
fT += aTB[i] * zn;
double izN = i * zN;
if ( ( i - N ) & 1 ) {
fp += apB[i] * izN;
fT += aTB[i] * izN;
} else {
fp -= apB[i] * izN;
fT -= aTB[i] * izN;
}
zn *= z;
}
f0 *= pole0 * logsB;
fp *= polep * logsB;
fT *= polep * logsB;
}
void EvtbTosllVectorAmpNP::CalcSAmp( EvtParticle* parent, EvtAmp& amp,
EvtbTosllFF* formFactors )
{
// Add the lepton and neutrino 4 momenta to find q2
EvtId pId = parent->getId();
EvtId dId = parent->getDaug( 0 )->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 );
}