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

#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtConst.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/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"

#include "EvtGenModels/EvtbTosllAmp.hh"
#include "EvtGenModels/EvtbTosllFF.hh"

void EvtbTosllScalarAmp::CalcAmp( EvtParticle* parent, EvtAmp& amp,
                                  EvtbTosllFF* formFactors )
{
    //Add the lepton and neutrino 4 momenta to find q2

    EvtVector4R q = parent->getDaug( 1 )->getP4() + parent->getDaug( 2 )->getP4();
    double q2 = ( q.mass2() );

    double fp( 0. ), f0( 0. ), ft( 0. );
    double mesonmass = parent->getDaug( 0 )->mass();
    double parentmass = parent->mass();

    formFactors->getScalarFF( parent->getId(), parent->getDaug( 0 )->getId(),
                              q2, mesonmass, fp, f0, ft );

    EvtId daught = parent->getDaug( 0 )->getId();
    bool btod = false;
    bool nnlo = true;
    if ( daught == EvtPDL::getId( std::string( "pi+" ) ) ||
         daught == EvtPDL::getId( std::string( "pi-" ) ) ||
         daught == EvtPDL::getId( std::string( "pi0" ) ) ||
         daught == EvtPDL::getId( std::string( "eta" ) ) ||
         daught == EvtPDL::getId( std::string( "eta'" ) ) )
        btod = true;

    EvtVector4R p4b;
    p4b.set( parent->mass(), 0.0, 0.0, 0.0 );

    EvtVector4C l11, l12;
    EvtVector4C l21, l22;

    EvtVector4C a11, a12;
    EvtVector4C a21, a22;

    EvtId l_num = parent->getDaug( 1 )->getId();

    EvtVector4C T1, T2;

    EvtVector4R phat = p4b / parentmass;
    EvtVector4R qhat = q / parentmass;

    EvtComplex c7eff = EvtbTosllAmp::GetC7Eff( q2, nnlo );
    EvtComplex c9eff = EvtbTosllAmp::GetC9Eff( q2, nnlo, btod );
    EvtComplex c10eff = EvtbTosllAmp::GetC10Eff( q2, nnlo );

    //double mbhat=1;
    double mbhat = 4.4 / ( parentmass );
    //double mkhat = 0.15;
    double mkhat = mesonmass / ( parentmass );
    double shat = q2 / ( parentmass * parentmass );

    double fm = ( f0 - fp ) * ( 1 - mkhat * mkhat ) / shat;

    EvtComplex aprime;
    aprime = c9eff * fp + 2.0 * mbhat * c7eff * ft / ( 1 + mkhat );
    EvtComplex bprime;
    bprime = c9eff * fm - 2 * mbhat * c7eff * ft * ( 1 - mkhat ) / shat;

    EvtComplex cprime;
    cprime = c10eff * fp;
    EvtComplex dprime;
    dprime = c10eff * fm;

    static EvtIdSet leptons( "e-", "mu-", "tau-" );
    static EvtIdSet antileptons( "e+", "mu+", "tau+" );

    if ( leptons.contains( l_num ) ) {
        T1 = aprime * phat + bprime * qhat;
        T2 = cprime * phat + dprime * qhat;

        l11 = EvtLeptonVCurrent( parent->getDaug( 1 )->spParent( 0 ),
                                 parent->getDaug( 2 )->spParent( 0 ) );
        l21 = EvtLeptonVCurrent( parent->getDaug( 1 )->spParent( 1 ),
                                 parent->getDaug( 2 )->spParent( 0 ) );
        l12 = EvtLeptonVCurrent( parent->getDaug( 1 )->spParent( 0 ),
                                 parent->getDaug( 2 )->spParent( 1 ) );
        l22 = EvtLeptonVCurrent( parent->getDaug( 1 )->spParent( 1 ),
                                 parent->getDaug( 2 )->spParent( 1 ) );
        a11 = EvtLeptonACurrent( parent->getDaug( 1 )->spParent( 0 ),
                                 parent->getDaug( 2 )->spParent( 0 ) );
        a21 = EvtLeptonACurrent( parent->getDaug( 1 )->spParent( 1 ),
                                 parent->getDaug( 2 )->spParent( 0 ) );
        a12 = EvtLeptonACurrent( parent->getDaug( 1 )->spParent( 0 ),
                                 parent->getDaug( 2 )->spParent( 1 ) );
        a22 = EvtLeptonACurrent( parent->getDaug( 1 )->spParent( 1 ),
                                 parent->getDaug( 2 )->spParent( 1 ) );
    } else {
        if ( antileptons.contains( l_num ) ) {
            T1 = aprime * phat + bprime * qhat;
            T2 = cprime * phat + dprime * qhat;

            l11 = EvtLeptonVCurrent( parent->getDaug( 1 )->spParent( 1 ),
                                     parent->getDaug( 2 )->spParent( 1 ) );
            l21 = EvtLeptonVCurrent( parent->getDaug( 1 )->spParent( 0 ),
                                     parent->getDaug( 2 )->spParent( 1 ) );
            l12 = EvtLeptonVCurrent( parent->getDaug( 1 )->spParent( 1 ),
                                     parent->getDaug( 2 )->spParent( 0 ) );
            l22 = EvtLeptonVCurrent( parent->getDaug( 1 )->spParent( 0 ),
                                     parent->getDaug( 2 )->spParent( 0 ) );

            a11 = EvtLeptonACurrent( parent->getDaug( 1 )->spParent( 1 ),
                                     parent->getDaug( 2 )->spParent( 1 ) );
            a21 = EvtLeptonACurrent( parent->getDaug( 1 )->spParent( 0 ),
                                     parent->getDaug( 2 )->spParent( 1 ) );
            a12 = EvtLeptonACurrent( parent->getDaug( 1 )->spParent( 1 ),
                                     parent->getDaug( 2 )->spParent( 0 ) );
            a22 = EvtLeptonACurrent( parent->getDaug( 1 )->spParent( 0 ),
                                     parent->getDaug( 2 )->spParent( 0 ) );

        } else {
            EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Wrong lepton number\n";
        }
    }

    amp.vertex( 0, 0, l11 * T1 + a11 * T2 );
    amp.vertex( 0, 1, l12 * T1 + a12 * T2 );
    amp.vertex( 1, 0, l21 * T1 + a21 * T2 );
    amp.vertex( 1, 1, l22 * T1 + a22 * T2 );
}