Evtbs2llGammaFFMNT.cpp

Go to the documentation of this file.

/***********************************************************************
* 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/Evtbs2llGammaFFMNT.hh"

#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4C.hh"

#include <cmath>
#include <cstdlib>

Evtbs2llGammaFFMNT::Evtbs2llGammaFFMNT()
{
}

/*                                                              *
 *  decay_id = 0 for b \bar q -> l^+ l^- \gamma transitions     *
 *             1 for q \bar b -> l^+ l^- \gamma transitions;    * 
 *  fb    - leptonic decay constant of the B_q - meson;         *
 *  mb    - the mass of the b-quark;                            *
 *  mq    - the mass of the light quark (d or s);               *
 *  c7gam - Wilson coefficient C_{7\gamma};                     *
 *  a1 = c1 + c2/3.0 - linear combination of the Wils. Coeff.;  *
 *  lambda_qu = V^*_{uq}*V_{ub}/V^*_{tq}*V_{tb}, where q={d,s}; *
 *  lambda_qc = V^*_{cq}*V_{cb}/V^*_{tq}*V_{tb}, where q={d,s}. *
 *                                                              */
void Evtbs2llGammaFFMNT::getPhotonFF( int decay_id, double fb, EvtId parent,
                                      double q2, double M1, double mb,
                                      double mq, EvtComplex c7gam, EvtComplex a1,
                                      EvtComplex lambda_qu, EvtComplex lambda_qc,
                                      EvtComplex& Fv, EvtComplex& Fa,
                                      EvtComplex& Ftv, EvtComplex& Fta )
{
    EvtComplex unit1( 1.0, 0.0 );    // real unit
    EvtComplex uniti( 0.0, 1.0 );    // imaginary unit
    EvtComplex unit0( 0.0, 0.0 );    // complex zero unit

    // characteristics of resonances rho, omega, phi
    double M_res[] = {0.7758, 0.78259, 1.019456};    // particle masses, Gev
    double Gamma[] = {0.1503, 0.00849, 0.00426};     // particle widthes, Gev

    double f_lept[] = {5.04, 17.1, -13.2};     // decay constants f_i
    double g_plus[] = {0.27, -0.27, -0.38};    // and form-factors g+(0)
    g_plus[0] = g_plus[0] / sqrt( 2.0 );    // by D.Melikhov, N.Nikitin, K.Toms,
    g_plus[1] = g_plus[1] / sqrt( 2.0 );    // Phys.At.Nucl. 68, p.1842 (2005)

    double hatq2 = q2 / pow( M1, 2 );
    // E - photon energy in the B-meson rest frame
    double E = 0.5 * M1 * ( 1 - hatq2 );

    // parametrs for form-factors Fv, Ftv, Fa, Fta
    //(by D.Melikhov, N.Nikitin, K.Toms, Yad. Fiz. 62, No 11)
    double beta[] = {0.28, 0.30, 0.26, 0.33};     // beta, Gev^(-1)
    double Delta[] = {0.04, 0.04, 0.30, 0.30};    // Delta, Gev

    // form-factors
    EvtComplex Ftvq0, Ftaq0, Ftv00, Fta00;
    Fv = unit1 * beta[0] * fb * M1 / ( Delta[0] + E );       // Fv(q^2)
    Ftvq0 = unit1 * beta[1] * fb * M1 / ( Delta[1] + E );    // Ftv(q^2,0)
    Fa = unit1 * beta[2] * fb * M1 / ( Delta[2] + E );       // Fa(q^2)
    Ftaq0 = unit1 * beta[3] * fb * M1 / ( Delta[3] + E );    // Fta(q^2,0)
    Ftv00 = unit1 * beta[1] * fb * M1 / ( Delta[1] + 0.5 * M1 );    // Ftv(0,0)
    Fta00 = unit1 * beta[3] * fb * M1 / ( Delta[3] + 0.5 * M1 );    // Fta(0,0)
    EvtComplex Ftv_WA( 0.0, 0.0 );    // the weak annihilation contribution

    // Resonant contribution to the form-factors Ftv(0,q^2) and Fta(0,q^2)
    EvtComplex ResSum( 0.0, 0.0 );

    if ( parent == EvtPDL::getId( std::string( "B_s0" ) ) ||
         parent == EvtPDL::getId( std::string( "anti-B_s0" ) ) ) {
        // only \phi-resonant contribution to the Bs-decays
        ResSum = 2.0 * g_plus[2] * q2 /
                 ( f_lept[2] * ( unit1 * ( q2 - pow( M_res[2], 2 ) ) +
                                 uniti * M_res[2] * Gamma[2] ) );
    }
    if ( parent == EvtPDL::getId( std::string( "B_d0" ) ) ||
         parent == EvtPDL::getId( std::string( "anti-B_d0" ) ) ) {
        // \rho- and \omega-resonant contribution to the Bd-decays
        for ( int i = 0; i < 2; i++ ) {
            ResSum = ResSum +
                     2.0 * g_plus[i] * q2 /
                         ( f_lept[i] * ( unit1 * ( q2 - pow( M_res[i], 2 ) ) +
                                         uniti * M_res[i] * Gamma[i] ) );
        }
    }

    EvtComplex Ftv0q = Ftv00 - ResSum;    // form-factor Ftv(0,q^2)
    EvtComplex Fta0q = Fta00 - ResSum;    // form-factor Fta(0,q^2)

    // Ftv(q^2,q^2) = Ftv(q^2,0)+Ftv(0,q^2)
    Ftv = Ftvq0 + Ftv0q;
    // Fta(q^2,q^2) = Fta(q^2,0)+Fta(0,q^2)
    Fta = Ftaq0 + Fta0q;

    // Weak annihilation
    if ( abs( c7gam ) < 0.0000001 ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "\n\n The function Evtbs2llGammaFFMNT::getPhotonFF"
            << "\n Error: the Wilson coefficient C7gamma = 0!"
            << " c7gam = " << c7gam << std::endl;
        ::abort();
    }
    if ( mb < 0.001 ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "\n\n The function Evtbs2llGammaFFMNT::getPhotonFF"
            << " mb = " << mb << " << 5 GeV!" << std::endl;
        ::abort();
    }

    switch ( decay_id ) {
        /* b \bar q -> l^+ l^- \gamma transitions */
        case 0:
            Ftv_WA = ( 16.0 / 3.0 ) * ( lambda_qu + lambda_qc ) *
                     ( a1 / c7gam ) * ( fb / mb );
            Ftv = ( 1.0 + mq / mb ) * Ftv - Ftv_WA;
            Fta = ( 1.0 - mq / mb ) * Fta;
            //Fv  = Fv;
            //Fa  = Fa;
            break;
        /* q \bar b -> l^+ l^- \gamma transitions */
        case 1:
            Ftv_WA = ( 16.0 / 3.0 ) * conj( lambda_qu + lambda_qc ) *
                     ( a1 / c7gam ) * ( fb / mb );
            Ftv = ( 1.0 + mq / mb ) * Ftv + Ftv_WA;
            Fta = ( 1.0 - mq / mb ) * Fta;    // The change of the sign
            //Fv  = Fv;                         // is included in the
            //Fa  = Fa;                         // amplitudes definition!

            break;
    };
}

// Getting the quark mass (in GeV) using to the dispersion quark model
// of D.Melikhov, B.Stech, PRD62, 014006 (2000).
//
// i=1 => return m_u;
// i=2 => return m_d;
// i=3 => return m_s;
// i=4 => return m_c;
// i=5 => return m_b;
double Evtbs2llGammaFFMNT::getQuarkMass( int i )
{
    double qm = 0.0;

    switch ( i ) {
        case 1:
            qm = 0.23;    // m_u
            break;
        case 2:
            qm = 0.23;    // m_d = m_u
            break;
        case 3:
            qm = 0.35;    // m_s
            break;
        case 4:
            qm = 1.45;    // m_c
            break;
        case 5:
            qm = 4.85;    // m_b
            break;
        default:
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "In the function EvtbTosllMSFF::getQuarkMass   \n"
                << "the parametr i not equal 1, 2,  3, 4 or 5! \n"
                << "i =" << i << std::endl;
            ::abort();
    }

    return qm;
}