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

#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtIdSet.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtScalarParticle.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtVectorParticle.hh"

#include "EvtGenModels/EvtbTosllWilsCoeffNLO.hh"
#include "EvtGenModels/Evtbs2llGammaFFMNT.hh"

#include <cstdlib>

// input:   *parent      - the pointer to the parent particle (B-meson, the
//                                          object of the EvtParticle class);
//          *formFactors - the pointer to instance of EvtbTosllGammaFF class object;
//          *WilsCoeff   - the pointer to the Standart Model Wilson Coefficients class;
//           mu          - the scale parameter, GeV;
//           Nf          - number of "effective" flavors (for b-quark Nf=5);
//           res_swch    - resonant switching parameter:
//                         = 0 the resonant contribution switched OFF,
//                         = 1 the resonant contribution switched ON;
//           ias         - switching parameter for \alpha_s(M_Z) value:
//                         = 0 PDG 1sigma minimal alpha_s(M_Z),
//                         = 1 PDG average value  alpha_s(M_Z),
//                         = 2 PDG 1sigma maximal alpha_s(M_Z).
//           Egamma_min  - photon energy cut, GeV;
//           Wolfenstein parameterization for CKM matrix
//                         CKM_A, CKM_lambda, CKM_barrho, CKM_bareta

void Evtbs2llGammaAmp::CalcAmp( EvtParticle* parent, EvtAmp& amp,
                                Evtbs2llGammaFF* formFactors,
                                EvtbTosllWilsCoeffNLO* WilsCoeff, double mu,
                                int Nf, int res_swch, int ias, double Egamma_min,
                                double CKM_A, double CKM_lambda,
                                double CKM_barrho, double CKM_bareta )
{
    //  FILE *mytest;

    int iG = 0;    // photon is the first daughter particle
    int il1 = 1,
        il2 = 2;    // leptons are the second and thirds daughter particles

    EvtComplex unit1( 1.0, 0.0 );    // real unit
    EvtComplex uniti( 0.0, 1.0 );    // imaginary unit

    double M1 = parent->mass();                    // B - meson mass, GeV
    double ml = parent->getDaug( il1 )->mass();    // leptonic mass, GeV
    double mq = 0.0;    // light quark mass from the dispersion QM, GeV
    double mc = formFactors->getQuarkMass(
        4 );    // m_c mass from the dispersion QM, GeV
    double mb = formFactors->getQuarkMass(
        5 );               // m_b mass from the dispersion QM, GeV
    double Mw = 80.403;    // GeV W-boson mass, GeV
    double mt = 174.2;     // GeV t-quark mass, GeV
    double fb = 0.0;       // leptonic decay constant of B-meson, Gev

    EvtComplex Vtb, Vtq, Vub, Vuq, Vcb,
        Vcq;                  // V_{tb}, V_{tq}, V_{ub}, V_{uq}, V_{cb}, V_{cq}
    EvtComplex CKM_factor;    // V^*_{tq}*V_{tb}, where q={d,s}
    EvtComplex lambda_qu;     // V^*_{uq}*V_{ub}/V^*_{tq}*V_{tb}, where q={d,s}
    EvtComplex lambda_qc;     // V^*_{cq}*V_{cb}/V^*_{tq}*V_{tb}, where q={d,s}
    double Relambda_qu, Imlambda_qu;

    // to find charges of ell^+ and ell^- in the B-meson  daughters
    int charge1 = EvtPDL::chg3( parent->getDaug( il1 )->getId() );
    int charge2 = EvtPDL::chg3( parent->getDaug( il2 )->getId() );
    if ( charge1 == 0 || charge2 == 0 ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "\n\n The function EvtbsTollGammaAmp::CalcAmp(...)"
            << "\n Error in the leptonic charge getting!"
            << "\n charge1             =" << charge1
            << "\n charge2             =" << charge2 << "\n charge gamma        ="
            << EvtPDL::chg3( parent->getDaug( iG )->getId() )
            << "\n number of daughters =" << parent->getNDaug() << std::endl;
        ::abort();
    }

    EvtParticle* lepPlus = 0;
    EvtParticle* lepMinus = 0;

    lepPlus = ( charge1 > charge2 )
                  ? parent->getDaug( il1 )
                  : parent->getDaug( il2 );    // positive charged
    lepMinus = ( charge1 < charge2 )
                   ? parent->getDaug( il1 )
                   : parent->getDaug( il2 );    // negative charged

    EvtVector4R p = parent->getP4Restframe();    // B-meson momentum in the B-rest frame
    EvtVector4R k =
        parent->getDaug( iG )->getP4();    // 4-momentum of photon in the B-rest frame
    EvtVector4R q = p - k;    // transition 4-momentum q=p-k in the B-rest frame
    EvtVector4R p_1;          // 4-momentum of ell^+ in the B-rest frame
    EvtVector4R p_2;          // 4-momentum of ell^- in the B-rest frame

    // the preparation of the leptonic 4-momentums in the B-rest frame
    if ( charge1 > charge2 ) {
        p_1 = parent->getDaug( il1 )->getP4();
        p_2 = parent->getDaug( il2 )->getP4();
    } else {
        p_1 = parent->getDaug( il2 )->getP4();
        p_2 = parent->getDaug( il1 )->getP4();
    }

    EvtVector4R p_minus_p_1 =
        p - p_1;    // transition momentum of the B-meson and antilepton p-p_1
    EvtVector4R p_minus_p_2 =
        p - p_2;    // transition momentum of the B-meson and lepton p-p_2

    double q2 = q.mass2();             // Mandelstam variable s=q^2
    double p2 = p.mass2();             // p^2=M1^2
    double t = p_minus_p_1.mass2();    // Mandelstam variable t=(p-p_1)^2
    double u = p_minus_p_2.mass2();    // Mandelstam variable u=(p-p_2)^2

    // scalar products
    double pk = 0.5 * ( p2 - q2 );                // (p*k)
    double p1k = 0.5 * ( pow( ml, 2.0 ) - u );    // (p1*k)
    double p2k = 0.5 * ( pow( ml, 2.0 ) - t );    // (p2*k)

    double hatq2 = q2 / ( M1 * M1 );    // \hat s = q^2/M_1^2
    double Egam = 0.5 * M1 *
                  ( 1 - hatq2 );    // photon energy in the B-meson rest frame

    EvtVector4R hatp = p / M1;
    EvtVector4R hatk = k / M1;

    //  EvtGenReport(EVTGEN_NOTICE,"EvtGen")
    //         << "\n\n The function EvtbsTollGammaAmp::CalcAmp(...)"
    //         << "\n q =  p-k  =" << p-k     << "   q^2 = " << (p-k).mass2()
    //         << "\n q = p1+p2 =" << p_1+p_2 << "   q^2 = " << (p_1+p_2).mass2()
    //         << "\n m_ell     =" << parent->getDaug(il1)->mass()
    //         << "\n m_ell     =" << parent->getDaug(il2)->mass()
    //         << "\n m_gamma   =" << parent->getDaug(iG)->mass()
    //         << std::endl;

    EvtId idparent = parent->getId();    // B-meson Id

    if ( idparent == EvtPDL::getId( std::string( "B_s0" ) ) ||
         idparent == EvtPDL::getId( std::string( "anti-B_s0" ) ) ) {
        mq = formFactors->getQuarkMass( 3 );    // m_s mass from the dispersion QM
        fb = 0.24;                              // leptonic decay constant
            // V_{ts}
        Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) +
              pow( CKM_lambda, 2.0 ) *
                  ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
                  sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
        Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq;
        // V_{us}
        Vuq = CKM_lambda * unit1;
        // V_{cs}
        Vcq = unit1 - 0.5 * pow( CKM_lambda, 2.0 ) -
              0.125 * pow( CKM_lambda, 4.0 ) * ( 1.0 + 4.0 * pow( CKM_A, 2.0 ) );
    }

    if ( idparent == EvtPDL::getId( std::string( "B0" ) ) ||
         idparent == EvtPDL::getId( std::string( "anti-B0" ) ) ) {
        mq = formFactors->getQuarkMass( 2 );    // m_d mass from the dispersion QM
        fb = 0.20;                              // leptonic decay constant
            // V_{td}
        Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) *
                          ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
                          sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
        Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq;
        // V_{ud}
        Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) -
                        0.125 * pow( CKM_lambda, 4.0 ) );
        // V_{cd}
        Vcq = unit1 *
              ( -CKM_lambda +
                0.5 * pow( CKM_A, 2.0 ) * pow( CKM_lambda, 5.0 ) *
                    ( 1.0 - 2.0 * ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
                                sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ) ) );
    }

    if ( mq < 0.001 ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "\n\n The function EvtbsTollGammaAmp::CalcAmp(..// 4-momentum of ell^+.)"
            << "\n Error in the model set!"
            << " mq = " << mq << std::endl;
        ::abort();
    }

    Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda,
                                     2.0 ) );    // V_{tb}
    Vub = CKM_A * pow( CKM_lambda, 3.0 ) *
          ( CKM_barrho * unit1 - CKM_bareta * uniti ) /
          sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );      // V_{ub}
    Vcb = unit1 * CKM_A * pow( CKM_lambda, 2.0 );    // V_{cb}

    CKM_factor = conj( Vtq ) * Vtb;    // V^*_{tq}*V_{tb}

    lambda_qu = conj( Vuq ) * Vub /
                CKM_factor;    // V^*_{uq}*V_{ub}/V^*_{tq}*V_{tb}
    Relambda_qu = real( lambda_qu );
    Imlambda_qu = imag( lambda_qu );

    lambda_qc = conj( Vcq ) * Vcb /
                CKM_factor;    // V^*_{cq}*V_{cb}/V^*_{tq}*V_{tb}

    // The Wilson Coefficients preparation according to the paper
    // A.J.Buras, M.Munz, Phys.Rev.D52, p.189 (1995)
    double c1, c2;
    EvtComplex a1, c7gam, c9eff_b2q, c9eff_barb2barq, c10a;

    // foton energy cut and removal of the J/psi amd psi' resonant area
    if ( Egam < Egamma_min || ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ) {
        c1 = 0.0;
        c2 = 0.0;
        a1 = unit1 * 0.0;
        c7gam = unit1 * 0.0;
        c9eff_b2q = unit1 * 0.0;
        c9eff_barb2barq = unit1 * 0.0;
        c10a = unit1 * 0.0;
    } else {
        c1 = WilsCoeff->C1( mu, Mw, Nf, ias );
        c2 = WilsCoeff->C2( mu, Mw, Nf, ias );
        a1 = unit1 * ( c1 + c2 / 3.0 );
        c7gam = WilsCoeff->GetC7Eff( mu, Mw, mt, Nf, ias );
        c9eff_b2q = WilsCoeff->GetC9Eff( 0, res_swch, ias, Nf, q2, mb, mq, mc, mu,
                                         mt, Mw, ml, Relambda_qu, Imlambda_qu );
        c9eff_barb2barq = WilsCoeff->GetC9Eff( 1, res_swch, ias, Nf, q2, mb, mq,
                                               mc, mu, mt, Mw, ml, Relambda_qu,
                                               Imlambda_qu );
        c10a = WilsCoeff->GetC10Eff( mt, Mw );
    }

    EvtComplex Fv,
        Fa;    // The change of the sign is included in the amplitudes definition!
    EvtComplex Ftv_b2q, Ftv_barb2barq;
    EvtComplex Fta_b2q, Fta_barb2barq;

    // foton energy cut and removal of the J/psi amd psi' resonant area
    if ( Egam < Egamma_min || ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ) {
        fb = 0.0;
        Fa = unit1 * 0.0;
        Fv = unit1 * 0.0;
        Fta_b2q = unit1 * 0.0;
        Fta_barb2barq = unit1 * 0.0;
        Ftv_b2q = unit1 * 0.0;
        Ftv_barb2barq = unit1 * 0.0;
    } else {
        if ( fb < 0.01 ) {
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "\n\n The function EvtbsTollGammaAmp::CalcAmp(...)"
                << "\n Leptonic decay constant fb is not uninitialized in this function!"
                << " fb = " << fb << std::endl;
            ::abort();
        }

        // For \bar B^0_q -> l^+ l^- gamma
        formFactors->getPhotonFF( 0, fb, parent->getId(), q2, M1, mb, mq, c7gam,
                                  a1, lambda_qu, lambda_qc, Fv, Fa, Ftv_b2q,
                                  Fta_b2q );

        // For B^0_q -> l^+ l^- gamma
        formFactors->getPhotonFF( 1, fb, parent->getId(), q2, M1, mb, mq, c7gam,
                                  a1, lambda_qu, lambda_qc, Fv, Fa,
                                  Ftv_barb2barq, Fta_barb2barq );
    }

    //  EvtGenReport(EVTGEN_NOTICE,"EvtGen")
    //      << "\n ============================================================================"
    //      << "\n ============================================================================"
    //      << "\n\n The function Evtbs2llGammaAmp::CalcAmp(...) passed."
    //      << "\n Particle masses:"
    //      << "\n B - meson mass M1 = " << M1
    //      << "\n photon minimum E  = " << Egamma_min
    //      << "\n q2                = " << q2
    //      << "\n leptonic mass  ml = " << ml
    //      << "\n light quark mass  = " << mq
    //      << "\n c - quark mass mc = " << mc
    //      << "\n b - quark mass mb = " << mb
    //      << "\n t - quark mass mt = " << mt
    //      << "\n W - boson mass Mw = " << Mw
    //      << "\n ============================================================================"
    //      << "\n Input parameters:"
    //      << "\n scale parameter         mu = " << mu
    //      << "\n number of flavors       Nf = " << Nf
    //      << "\n resonant switching         = " << res_swch
    //      << "\n parameter for alpha_s(M_Z) = " << ias
    //      << "\n photon energy cut (GeV)    = " << Egamma_min
    //      << "\n ============================================================================"
    //      << "\n Form-factors"
    //      << "\n Egam          = " << Egam
    //      << "\n Egamma_min    = " << Egamma_min
    //      << "\n Fv            = " << Fv
    //      << "\n Fa            = " << Fa
    //      << "\n Ftv_b2q       = " << Ftv_b2q
    //      << "\n Fta_b2q       = " << Fta_b2q
    //      << "\n Ftv_barb2barq = " << Ftv_barb2barq
    //      << "\n Fta_barb2barq = " << Fta_barb2barq
    //      << "\n ============================================================================"
    //      << "\n Wilson Coefficients:"
    //      << "\n Egam                = " << Egam
    //      << "\n Egamma_min          = " << Egamma_min
    //      << "\n Re(c7gam)           = " << real(c7gam)
    //        << " Im(c7gam)           = " << imag(c7gam)
    //      << "\n Re(c9eff_b2q)       = " << real(c9eff_b2q)
    //        << " Im(c9eff_b2q)       = " << imag(c9eff_b2q)
    //      << "\n Re(c9eff_barb2barq) = " << real(c9eff_barb2barq)
    //        << " Im(c9eff_barb2barq) = " << imag(c9eff_barb2barq)
    //      << "\n Re(c10a)            = " << real(c10a)
    //        << " Im(c10a)            = " << imag(c10a)
    //      << std::endl;

    // Hadronic matrix element coefficients
    EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, e_b2q, e_barb2barq,
        f_b2q, f_barb2barq;
    EvtComplex brammS, brammT;

    a_b2q = c9eff_b2q * Fv + 2.0 * c7gam * Ftv_b2q * mb * M1 / q2;
    a_barb2barq = c9eff_barb2barq * Fv +
                  2.0 * c7gam * Ftv_barb2barq * mb * M1 / q2;

    b_b2q = ( c9eff_b2q * Fa + 2.0 * c7gam * Fta_b2q * mb * M1 / q2 ) * pk /
            ( M1 * M1 );
    b_barb2barq = ( c9eff_barb2barq * Fa +
                    2.0 * c7gam * Fta_barb2barq * mb * M1 / q2 ) *
                  pk / ( M1 * M1 );

    e_b2q = c10a * Fv;
    e_barb2barq = e_b2q;

    f_b2q = c10a * Fa * pk / ( M1 * M1 );
    f_barb2barq = f_b2q;

    brammS = 0.0;    // in the Bq-meson rest frame!
    brammT = 0.5 * c10a * ml * fb *
             ( 1.0 / p2k + 1.0 / p1k );    // for Bramsstrahlung

    EvtTensor4C T1, T2;    // hadronic matrix element tensor structures
    EvtVector4C E1, E2;
    EvtComplex E3;

    EvtVector4C epsG;    // photon polarisation vector

    int i;    // photon polarisations counter

    EvtVector4C lvc11, lvc12;    // spin structures for
    EvtVector4C lvc21, lvc22;    // the leptonic vector current

    EvtVector4C lac11, lac12;    // spin structures for
    EvtVector4C lac21, lac22;    // the leptonic axial current

    EvtComplex lsc11, lsc12;    // spin structures for
    EvtComplex lsc21, lsc22;    // the leptonic scalar current

    EvtTensor4C ltc11, ltc12;    // spin structures for
    EvtTensor4C ltc21, ltc22;    // the leptonic tensor current

    // B - and barB - mesons descriptors
    static EvtIdSet bmesons( "anti-B0", "anti-B_s0" );
    static EvtIdSet bbarmesons( "B0", "B_s0" );

    EvtId parentID = parent->getId();

    if ( bmesons.contains( parentID ) ) {
        // The amplitude for the decay barB -> gamma ell^+ ell^-  or
        // b \bar q -> gamma ell^+ ell^-

        T1 = -a_b2q * unit1 * dual( EvtGenFunctions::directProd( hatp, hatk ) ) -
             b_b2q * uniti * EvtTensor4C::g();

        T2 = -e_b2q * unit1 * dual( EvtGenFunctions::directProd( hatp, hatk ) ) -
             f_b2q * uniti * EvtTensor4C::g();

        // spin combinations for vector lepton current
        lvc11 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
                                   lepMinus->spParent( 0 ) );
        lvc21 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
                                   lepMinus->spParent( 0 ) );
        lvc12 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
                                   lepMinus->spParent( 1 ) );
        lvc22 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
                                   lepMinus->spParent( 1 ) );

        lac11 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
                                   lepMinus->spParent( 0 ) );
        lac21 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
                                   lepMinus->spParent( 0 ) );
        lac12 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
                                   lepMinus->spParent( 1 ) );
        lac22 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
                                   lepMinus->spParent( 1 ) );

        lsc11 = EvtLeptonSCurrent( lepPlus->spParent( 0 ),
                                   lepMinus->spParent( 0 ) );
        lsc21 = EvtLeptonSCurrent( lepPlus->spParent( 1 ),
                                   lepMinus->spParent( 0 ) );
        lsc12 = EvtLeptonSCurrent( lepPlus->spParent( 0 ),
                                   lepMinus->spParent( 1 ) );
        lsc22 = EvtLeptonSCurrent( lepPlus->spParent( 1 ),
                                   lepMinus->spParent( 1 ) );

        // \epsilon^{\alpha\beta\mu\nu}*TCurrent_{\mu\nu}
        ltc11 = dual( EvtLeptonTCurrent( lepPlus->spParent( 0 ),
                                         lepMinus->spParent( 0 ) ) );
        ltc21 = dual( EvtLeptonTCurrent( lepPlus->spParent( 1 ),
                                         lepMinus->spParent( 0 ) ) );
        ltc12 = dual( EvtLeptonTCurrent( lepPlus->spParent( 0 ),
                                         lepMinus->spParent( 1 ) ) );
        ltc22 = dual( EvtLeptonTCurrent( lepPlus->spParent( 1 ),
                                         lepMinus->spParent( 1 ) ) );

        // summing up photon polarisations
        for ( i = 0; i < 2; i++ ) {
            // conjaction of epsG (photon polarization vector)
            EvtVector4C epsG = parent->getDaug( 0 )->epsParentPhoton( i ).conj();

            // de-escalation T with epsG
            E1 = T1.cont2( epsG );
            E2 = T2.cont2( epsG );
            E3 = ( epsG * hatp ) * brammS;

            // foton energy cut and removal of the J/psi amd psi' resonant area
            if ( Egam < Egamma_min ||
                 ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ) {
                CKM_factor = 0.0 * unit1;
            }

            // 1
            amp.vertex( i, 0, 0,
                        CKM_factor *
                            ( lvc11 * E1 + lac11 * E2 + uniti * lsc11 * E3 +
                              uniti * ( ( ltc11.cont2( hatp ) ) * epsG ) *
                                  brammT ) );
            //      EvtGenReport(EVTGEN_NOTICE,"EvtGen")
            //        << "\n 1" << CKM_factor*(lvc11*E1+lac11*E2+uniti*lsc11*E3+uniti*((ltc11.cont2(hatp))*epsG)*brammT)
            //        << std::endl;

            // 2
            amp.vertex( i, 0, 1,
                        CKM_factor *
                            ( lvc12 * E1 + lac12 * E2 + uniti * lsc12 * E3 +
                              uniti * ( ( ltc12.cont2( hatp ) ) * epsG ) *
                                  brammT ) );
            //      EvtGenReport(EVTGEN_NOTICE,"EvtGen")
            //        << "\n 2" << CKM_factor*(lvc12*E1+lac12*E2+uniti*lsc12*E3+uniti*((ltc12.cont2(hatp))*epsG)*brammT)
            //        << std::endl;

            // 3
            amp.vertex( i, 1, 0,
                        CKM_factor *
                            ( lvc21 * E1 + lac21 * E2 + uniti * lsc21 * E3 +
                              uniti * ( ( ltc21.cont2( hatp ) ) * epsG ) *
                                  brammT ) );
            //      EvtGenReport(EVTGEN_NOTICE,"EvtGen")
            //        << "\n 3" << CKM_factor*(lvc21*E1+lac21*E2+uniti*lsc21*E3+uniti*((ltc21.cont2(hatp))*epsG)*brammT)
            //        << std::endl;

            // 4
            amp.vertex( i, 1, 1,
                        CKM_factor *
                            ( lvc22 * E1 + lac22 * E2 + uniti * lsc22 * E3 +
                              uniti * ( ( ltc22.cont2( hatp ) ) * epsG ) *
                                  brammT ) );
            //      EvtGenReport(EVTGEN_NOTICE,"EvtGen")
            //        << "\n 4" << CKM_factor*(lvc22*E1+lac22*E2+uniti*lsc22*E3+uniti*((ltc22.cont2(hatp))*epsG)*brammT)
            //        << std::endl;
        }

    } else {
        if ( bbarmesons.contains( parentID ) ) {
            // The amplitude for the decay B -> gamma ell^+ ell^- or
            // q bar b -> gamma ell^+ ell^-

            T1 = -a_barb2barq * unit1 *
                     dual( EvtGenFunctions::directProd( hatp, hatk ) ) +
                 b_barb2barq * uniti * EvtTensor4C::g();

            T2 = -e_barb2barq * unit1 *
                     dual( EvtGenFunctions::directProd( hatp, hatk ) ) +
                 f_barb2barq * uniti * EvtTensor4C::g();

            lvc11 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
                                       lepMinus->spParent( 1 ) );
            lvc21 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
                                       lepMinus->spParent( 1 ) );
            lvc12 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
                                       lepMinus->spParent( 0 ) );
            lvc22 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
                                       lepMinus->spParent( 0 ) );

            lac11 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
                                       lepMinus->spParent( 1 ) );
            lac21 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
                                       lepMinus->spParent( 1 ) );
            lac12 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
                                       lepMinus->spParent( 0 ) );
            lac22 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
                                       lepMinus->spParent( 0 ) );

            lsc11 = EvtLeptonSCurrent( lepPlus->spParent( 1 ),
                                       lepMinus->spParent( 1 ) );
            lsc21 = EvtLeptonSCurrent( lepPlus->spParent( 0 ),
                                       lepMinus->spParent( 1 ) );
            lsc12 = EvtLeptonSCurrent( lepPlus->spParent( 1 ),
                                       lepMinus->spParent( 0 ) );
            lsc22 = EvtLeptonSCurrent( lepPlus->spParent( 0 ),
                                       lepMinus->spParent( 0 ) );

            // \epsilon^{\alpha\beta\mu\nu}*TCurrent_{\mu\nu}
            ltc11 = dual( EvtLeptonTCurrent( lepPlus->spParent( 1 ),
                                             lepMinus->spParent( 1 ) ) );
            ltc21 = dual( EvtLeptonTCurrent( lepPlus->spParent( 0 ),
                                             lepMinus->spParent( 1 ) ) );
            ltc12 = dual( EvtLeptonTCurrent( lepPlus->spParent( 1 ),
                                             lepMinus->spParent( 0 ) ) );
            ltc22 = dual( EvtLeptonTCurrent( lepPlus->spParent( 0 ),
                                             lepMinus->spParent( 0 ) ) );

            // summing up photon polarisations
            for ( i = 0; i < 2; i++ ) {
                EvtVector4C barepsG = parent->getDaug( 0 )->epsParentPhoton( i );

                E1 = T1.cont2( barepsG );
                E2 = T2.cont2( barepsG );
                E3 = ( barepsG * hatp ) * brammS;

                // foton energy cut and removal of the J/psi amd psi' resonant area
                if ( Egam < Egamma_min ||
                     ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ) {
                    CKM_factor = 0.0 * unit1;
                }

                amp.vertex(
                    i, 1, 1,
                    conj( CKM_factor ) *
                        ( lvc11 * E1 + lac11 * E2 + uniti * lsc11 * E3 +    // -?
                          uniti * ( ( ltc11.cont2( hatp ) ) * epsG ) * brammT ) );
                amp.vertex(
                    i, 1, 0,
                    conj( CKM_factor ) *
                        ( lvc12 * E1 + lac12 * E2 + uniti * lsc12 * E3 +    // -?
                          uniti * ( ( ltc12.cont2( hatp ) ) * epsG ) * brammT ) );
                amp.vertex(
                    i, 0, 1,
                    conj( CKM_factor ) *
                        ( lvc21 * E1 + lac21 * E2 + uniti * lsc21 * E3 +    // -?
                          uniti * ( ( ltc21.cont2( hatp ) ) * epsG ) * brammT ) );
                amp.vertex(
                    i, 0, 0,
                    conj( CKM_factor ) *
                        ( lvc22 * E1 + lac22 * E2 + uniti * lsc22 * E3 +    // -?
                          uniti * ( ( ltc22.cont2( hatp ) ) * epsG ) * brammT ) );
            }

        } else {
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "\n\n The function Evtbs2llGammaAmp::CalcAmp(...)"
                << "\n Wrong B-meson number" << std::endl;
            ::abort();
        }
    }
}

//
// The decays B -> Gamma ell^+ ell^- maximum probability calculation for the
// d^2\Gamma/dq^2 d\cos\theta distribution.
//
// \theta - the angle between the photon and ell^- directions in the
//          B-meson rest frame.
//
// If ias=0 (nonresonant case), the maximum is achieved at q2
//           B0s:                q2 = 4*ml^2, Mphi^2, q_max^2;
//           B0d:                q2 = 4*ml^2, Mrho^2, Momega^2, q_max^2;
// If ias=1 (resonant case), the maximum in the same points, because the
//          resonat area is remove
//
double Evtbs2llGammaAmp::CalcMaxProb( EvtId parnum, EvtId photnum, EvtId l1num,
                                      EvtId l2num, Evtbs2llGammaFF* formFactors,
                                      EvtbTosllWilsCoeffNLO* WilsCoeff,
                                      double mu, int Nf, int res_swch, int ias,
                                      double Egamma_min, double CKM_A,
                                      double CKM_lambda, double CKM_barrho,
                                      double CKM_bareta )
{
    double maxfoundprob = -100.0;    // maximum of the probability
                                     //  int ijk_atmax = 1000000;

    double M1 = EvtPDL::getMeanMass( parnum );    // B - meson mass
    double ml = EvtPDL::getMeanMass( l1num );     // leptonic mass

    double Mrho = EvtPDL::getMeanMass(
        EvtPDL::getId( std::string( "rho0" ) ) );    // mass of the rho-meson, MeV
    double Momega = EvtPDL::getMeanMass( EvtPDL::getId(
        std::string( "omega" ) ) );    // mass of the omega-meson, MeV
    double Mphi = EvtPDL::getMeanMass(
        EvtPDL::getId( std::string( "phi" ) ) );    // mass of the phi-meson, MeV

    // EvtGenReport(EVTGEN_NOTICE,"EvtGen")
    //   << "\n M1         = " << M1
    //   << "\n ml         = " << ml
    //   << "\n Mrho       = " << Mrho
    //   << "\n Momega     = " << Momega
    //   << "\n Mphi       = " << Mphi
    //   << "\n Egamma_min = " << Egamma_min
    //   << std::endl;

    double list_of_max_q2_points[5];
    list_of_max_q2_points[0] = pow( 2.0 * ml, 2.0 );
    list_of_max_q2_points[1] = pow( Mrho, 2.0 );
    list_of_max_q2_points[2] = pow( Momega, 2.0 );
    list_of_max_q2_points[3] = pow( Mphi, 2.0 );
    list_of_max_q2_points[4] =
        pow( M1, 2.0 ) - 2.0 * M1 * Egamma_min;    // q^2_max at photon energy cut

    //  if(list_of_max_points[4]<0){
    //     EvtGenReport(EVTGEN_ERROR,"EvtGen")
    //       << "\n\n In the function EvtbsTollGammaAmp::CalcScalarMaxProb(...)"
    //       << "\n Bad photon energy cut: Egamma_min > M1 in the rest frame of B-meson!"
    //       << "\n q2_max     = " << list_of_max_points[4]
    //       << "\n M1         = " << M1
    //       << "\n Egamma_min = " << Egamma_min
    //       << std::endl;
    //  ::abort();
    //  }

    if ( Egamma_min > Mrho ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "\n\n In the function EvtbsTollGammaAmp::CalcScalarMaxProb(...)"
            << "\n Bad photon energy cut: Egamma_min > M_rho0 in the rest frame of B-meson."
            << "\n Mrho       = " << Mrho << "\n Egamma_min = " << Egamma_min
            << std::endl;
        ::abort();
    }

    if ( Egamma_min <= 0 ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "\n\n In the function EvtbsTollGammaAmp::CalcScalarMaxProb(...)"
            << "\n Bad photon energy cut: Egamma_min <= 0 in the rest frame of B-meson."
            << "\n Egamma_min = " << Egamma_min << std::endl;
        ::abort();
    }

    if ( res_swch == 0 || res_swch == 1 ) {
        int i_list;
        for ( i_list = 0; i_list <= 4; i_list++ ) {
            double s;          // mandelstam variable "s";
            double t_minus;    // minimum and maximum of the mandelstam variable "t"
            double t_plus;    // as function of the mandelstam variable "s=q2";
            double t_for_s;
            int ijk;        // counter for variable "t";
            int max_ijk;    // maximal value of this counter;

            s = list_of_max_q2_points[i_list];

            t_plus = pow( M1, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
            t_plus = t_plus + sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
                                  ( pow( M1, 2.0 ) - s );
            t_plus *= 0.5;

            t_minus = pow( M1, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
            t_minus = t_minus - sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
                                    ( pow( M1, 2.0 ) - s );
            t_minus *= 0.5;

            if ( fabs( t_plus - t_minus ) < 0.000001 )
                t_minus = t_plus;

            max_ijk = 1000;
            double dt = ( t_plus - t_minus ) / ( (double)max_ijk );
            if ( fabs( dt ) < 0.00001 )
                dt = 0.0;

            if ( dt < 0.0 ) {
                EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                    << "\n\n In the function EvtbsTollGammaAmp::CalcScalarMaxProb(...)"
                    << "\n dt      = " << dt << " < 0."
                    << "\n s       = " << s << "\n t_plus  = " << t_plus
                    << "\n t_minus = " << t_minus << "\n M1      = " << M1
                    << "\n ml      = " << ml << std::endl;
                ::abort();
            }

            for ( ijk = 0; ijk <= max_ijk; ijk++ ) {
                t_for_s = t_minus + dt * ( (double)ijk );

                // B-meson rest frame particles and they kinematics inicialization
                double Eg, El2;
                Eg = ( pow( M1, 2.0 ) - s ) / ( 2.0 * M1 );    // photon energy
                El2 = ( s + t_for_s - pow( ml, 2.0 ) ) /
                      ( 2.0 * M1 );    // ell^- energy

                double modl2;
                modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );

                double cosBellminus;    // angle between the B-meson and ell^- directions
                cosBellminus = ( pow( ml, 2.0 ) + 2.0 * Eg * El2 - t_for_s ) /
                               ( 2.0 * Eg * modl2 );

                if ( ( fabs( cosBellminus ) > 1.0 ) &&
                     ( fabs( cosBellminus ) <= 1.0001 ) ) {
                    EvtGenReport( EVTGEN_NOTICE, "EvtGen" )
                        << "\n Debug in the function EvtbsTollGammaAmp::CalcMaxProb(...):"
                        << "\n cos(theta) = " << cosBellminus << std::endl;
                    cosBellminus = cosBellminus / fabs( cosBellminus );
                }
                if ( fabs( cosBellminus ) > 1.0001 ) {
                    EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                        << "\n\n In the function EvtbsTollGammaAmp::CalcMaxProb(...)"
                        << "\n |cos(theta)| = " << fabs( cosBellminus ) << " > 1"
                        << "\n s       = " << s << "\n t_for_s = " << t_for_s
                        << "\n t_plus  = " << t_plus << "\n t_minus = " << t_minus
                        << "\n dt      = " << dt << "\n Eg      = " << Eg
                        << "\n El2     = " << El2 << "\n modl2   = " << modl2
                        << "\n ml      = " << ml << std::endl;
                    ::abort();
                }

                EvtVector4R p, k, p1, p2;
                p.set( M1, 0.0, 0.0, 0.0 );
                k.set( Eg, Eg, 0.0, 0.0 );
                p2.set( El2, modl2 * cosBellminus,
                        -modl2 * sqrt( 1.0 - pow( cosBellminus, 2.0 ) ), 0.0 );
                p1 = p - k - p2;

                // B-meson state preparation at the rest frame of B-meson
                EvtScalarParticle* scalar_part;
                EvtParticle* root_part;
                scalar_part = new EvtScalarParticle;

                scalar_part->noLifeTime();
                scalar_part->init( parnum, p );
                root_part = (EvtParticle*)scalar_part;
                root_part->setDiagonalSpinDensity();

                // Amplitude initialization
                EvtId listdaug[3];
                listdaug[0] = photnum;
                listdaug[1] = l1num;
                listdaug[2] = l2num;

                EvtAmp amp;
                amp.init( parnum, 3, listdaug );

                // Daughters states preparation at the rest frame of B-meson
                root_part->makeDaughters( 3, listdaug );

                EvtParticle *gamm, *lep1, *lep2;
                gamm = root_part->getDaug( 0 );
                lep1 = root_part->getDaug( 1 );
                lep2 = root_part->getDaug( 2 );

                gamm->noLifeTime();
                lep1->noLifeTime();
                lep2->noLifeTime();

                gamm->init( photnum, k );
                lep1->init( l1num, p1 );
                lep2->init( l2num, p2 );

                EvtSpinDensity rho;
                rho.setDiag( root_part->getSpinStates() );

                // The amplitude calculation at the
                // "maximum amplitude" kinematical configuration
                CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf,
                         res_swch, ias, Egamma_min, CKM_A, CKM_lambda,
                         CKM_barrho, CKM_bareta );

                // Now find the probability at this q2 and cos theta lepton point
                double nikmax = rho.normalizedProb( amp.getSpinDensity() );

                if ( nikmax > maxfoundprob ) {
                    maxfoundprob = nikmax;
                    EvtGenReport( EVTGEN_NOTICE, "EvtGen" )
                        << "\n maxfoundprob ( s =" << s << ",  t = " << t_for_s
                        << " ) = " << maxfoundprob << "\n ijk =" << ijk
                        << std::endl;
                }

                delete scalar_part;
                //          delete root_part;
                delete gamm;
                delete lep1;
                delete lep2;

            }    // for(ijk=0; ijk<=max_ijk; ijk++)
        }        // i_list - variable loop
    }            // if(res_swch==0||res_swch==1)
    else {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "\n\n In the function EvtbsTollVectorAmpNew::CalcMaxProb(...)"
            << "\n Unexpected value of the variable res_swch !!!"
            << "\n res_swch = " << res_swch << std::endl;
        ::abort();
    }

    if ( maxfoundprob <= 0.0 ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "\n\n In the function EvtbsTollVectorAmpNew::CalcMaxProb(...)"
            << "\n maxfoundprob = " << maxfoundprob << " <0 or =0!"
            << "\n res_swch     = " << res_swch << std::endl;
        ::abort();
    }

    maxfoundprob *= 1.01;

    //  EvtGenReport(EVTGEN_NOTICE,"EvtGen")
    //         << "\n **********************************************************************"
    //         << "\n The function Evtbs2llGammaAmp::CalcMaxProb(...) passed with arguments:"
    //         << "\n mu =" << mu << " Nf =" << Nf
    //         << " res_swch =" << res_swch
    //         << " ias =" << ias
    //         << "\n CKM_A      = " << CKM_A
    //         << " CKM_lambda = " << CKM_lambda
    //         << "\n CKM_barrho = " << CKM_barrho
    //         << " CKM_bareta = " << CKM_bareta
    //         << "\n The distribution maximum maxfoundprob =" << maxfoundprob
    //         << "\n ijk = " << ijk_atmax
    //         << "\n **********************************************************************"
    //         << std::endl;

    return maxfoundprob;
}

// Triangular function
double Evtbs2llGammaAmp::lambda( double a, double b, double c )
{
    double l;

    l = pow( a, 2.0 ) + pow( b, 2.0 ) + pow( c, 2.0 ) - 2.0 * a * b -
        2.0 * a * c - 2.0 * b * c;

    return l;
}