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

#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtDiracParticle.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtGammaMatrix.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtRaritaSchwinger.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtSemiLeptonicFF.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"

#include <stdlib.h>

using std::endl;

void EvtSemiLeptonicBaryonAmp::CalcAmp( EvtParticle* parent, EvtAmp& amp,
                                        EvtSemiLeptonicFF* FormFactors )
{
    static EvtId EM = EvtPDL::getId( "e-" );
    static EvtId MUM = EvtPDL::getId( "mu-" );
    static EvtId TAUM = EvtPDL::getId( "tau-" );
    static EvtId EP = EvtPDL::getId( "e+" );
    static EvtId MUP = EvtPDL::getId( "mu+" );
    static EvtId TAUP = EvtPDL::getId( "tau+" );

    //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 f1v, f1a, f2v, f2a;
    double m_meson = parent->getDaug( 0 )->mass();

    FormFactors->getbaryonff( parent->getId(), parent->getDaug( 0 )->getId(),
                              q2, m_meson, &f1v, &f1a, &f2v, &f2a );

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

    EvtVector4C temp_00_term1;
    EvtVector4C temp_00_term2;

    EvtVector4C temp_01_term1;
    EvtVector4C temp_01_term2;

    EvtVector4C temp_10_term1;
    EvtVector4C temp_10_term2;

    EvtVector4C temp_11_term1;
    EvtVector4C temp_11_term2;

    EvtDiracSpinor p0 = parent->sp( 0 );
    EvtDiracSpinor p1 = parent->sp( 1 );

    EvtDiracSpinor d0 = parent->getDaug( 0 )->spParent( 0 );
    EvtDiracSpinor d1 = parent->getDaug( 0 )->spParent( 1 );

    temp_00_term1.set( 0, f1v * ( d0 * ( EvtGammaMatrix::g0() * p0 ) ) );
    temp_00_term2.set(
        0,
        f1a * ( d0 * ( ( EvtGammaMatrix::g0() * EvtGammaMatrix::g5() ) * p0 ) ) );
    temp_01_term1.set( 0, f1v * ( d0 * ( EvtGammaMatrix::g0() * p1 ) ) );
    temp_01_term2.set(
        0,
        f1a * ( d0 * ( ( EvtGammaMatrix::g0() * EvtGammaMatrix::g5() ) * p1 ) ) );
    temp_10_term1.set( 0, f1v * ( d1 * ( EvtGammaMatrix::g0() * p0 ) ) );
    temp_10_term2.set(
        0,
        f1a * ( d1 * ( ( EvtGammaMatrix::g0() * EvtGammaMatrix::g5() ) * p0 ) ) );
    temp_11_term1.set( 0, f1v * ( d1 * ( EvtGammaMatrix::g0() * p1 ) ) );
    temp_11_term2.set(
        0,
        f1a * ( d1 * ( ( EvtGammaMatrix::g0() * EvtGammaMatrix::g5() ) * p1 ) ) );

    temp_00_term1.set( 1, f1v * ( d0 * ( EvtGammaMatrix::g1() * p0 ) ) );
    temp_00_term2.set(
        1,
        f1a * ( d0 * ( ( EvtGammaMatrix::g1() * EvtGammaMatrix::g5() ) * p0 ) ) );
    temp_01_term1.set( 1, f1v * ( d0 * ( EvtGammaMatrix::g1() * p1 ) ) );
    temp_01_term2.set(
        1,
        f1a * ( d0 * ( ( EvtGammaMatrix::g1() * EvtGammaMatrix::g5() ) * p1 ) ) );
    temp_10_term1.set( 1, f1v * ( d1 * ( EvtGammaMatrix::g1() * p0 ) ) );
    temp_10_term2.set(
        1,
        f1a * ( d1 * ( ( EvtGammaMatrix::g1() * EvtGammaMatrix::g5() ) * p0 ) ) );
    temp_11_term1.set( 1, f1v * ( d1 * ( EvtGammaMatrix::g1() * p1 ) ) );
    temp_11_term2.set(
        1,
        f1a * ( d1 * ( ( EvtGammaMatrix::g1() * EvtGammaMatrix::g5() ) * p1 ) ) );

    temp_00_term1.set( 2, f1v * ( d0 * ( EvtGammaMatrix::g2() * p0 ) ) );
    temp_00_term2.set(
        2,
        f1a * ( d0 * ( ( EvtGammaMatrix::g2() * EvtGammaMatrix::g5() ) * p0 ) ) );
    temp_01_term1.set( 2, f1v * ( d0 * ( EvtGammaMatrix::g2() * p1 ) ) );
    temp_01_term2.set(
        2,
        f1a * ( d0 * ( ( EvtGammaMatrix::g2() * EvtGammaMatrix::g5() ) * p1 ) ) );
    temp_10_term1.set( 2, f1v * ( d1 * ( EvtGammaMatrix::g2() * p0 ) ) );
    temp_10_term2.set(
        2,
        f1a * ( d1 * ( ( EvtGammaMatrix::g2() * EvtGammaMatrix::g5() ) * p0 ) ) );
    temp_11_term1.set( 2, f1v * ( d1 * ( EvtGammaMatrix::g2() * p1 ) ) );
    temp_11_term2.set(
        2,
        f1a * ( d1 * ( ( EvtGammaMatrix::g2() * EvtGammaMatrix::g5() ) * p1 ) ) );

    temp_00_term1.set( 3, f1v * ( d0 * ( EvtGammaMatrix::g3() * p0 ) ) );
    temp_00_term2.set(
        3,
        f1a * ( d0 * ( ( EvtGammaMatrix::g3() * EvtGammaMatrix::g5() ) * p0 ) ) );
    temp_01_term1.set( 3, f1v * ( d0 * ( EvtGammaMatrix::g3() * p1 ) ) );
    temp_01_term2.set(
        3,
        f1a * ( d0 * ( ( EvtGammaMatrix::g3() * EvtGammaMatrix::g5() ) * p1 ) ) );
    temp_10_term1.set( 3, f1v * ( d1 * ( EvtGammaMatrix::g3() * p0 ) ) );
    temp_10_term2.set(
        3,
        f1a * ( d1 * ( ( EvtGammaMatrix::g3() * EvtGammaMatrix::g5() ) * p0 ) ) );
    temp_11_term1.set( 3, f1v * ( d1 * ( EvtGammaMatrix::g3() * p1 ) ) );
    temp_11_term2.set(
        3,
        f1a * ( d1 * ( ( EvtGammaMatrix::g3() * EvtGammaMatrix::g5() ) * p1 ) ) );

    EvtVector4C l1, l2;

    EvtId l_num = parent->getDaug( 1 )->getId();
    if ( l_num == EM || l_num == MUM || l_num == TAUM ) {
        l1 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 0 ),
                                 parent->getDaug( 2 )->spParentNeutrino() );
        l2 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 1 ),
                                 parent->getDaug( 2 )->spParentNeutrino() );
    } else {
        if ( l_num == EP || l_num == MUP || l_num == TAUP ) {
            l1 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
                                     parent->getDaug( 1 )->spParent( 0 ) );
            l2 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
                                     parent->getDaug( 1 )->spParent( 1 ) );
        } else {
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "Wrong lepton number" << endl;
        }
    }

    amp.vertex( 0, 0, 0, l1.cont( temp_00_term1 + temp_00_term2 ) );
    amp.vertex( 0, 0, 1, l2.cont( temp_00_term1 + temp_00_term2 ) );

    amp.vertex( 0, 1, 0, l1.cont( temp_01_term1 + temp_01_term2 ) );
    amp.vertex( 0, 1, 1, l2.cont( temp_01_term1 + temp_01_term2 ) );

    amp.vertex( 1, 0, 0, l1.cont( temp_10_term1 + temp_10_term2 ) );
    amp.vertex( 1, 0, 1, l2.cont( temp_10_term1 + temp_10_term2 ) );

    amp.vertex( 1, 1, 0, l1.cont( temp_11_term1 + temp_11_term2 ) );
    amp.vertex( 1, 1, 1, l2.cont( temp_11_term1 + temp_11_term2 ) );

    return;
}

double EvtSemiLeptonicBaryonAmp::CalcMaxProb( EvtId parent, EvtId baryon,
                                              EvtId lepton, EvtId nudaug,
                                              EvtSemiLeptonicFF* FormFactors,
                                              EvtComplex r00, EvtComplex r01,
                                              EvtComplex r10, EvtComplex r11 )
{
    //This routine takes the arguements parent, baryon, and lepton
    //number, and a form factor model, and returns a maximum
    //probability for this semileptonic form factor model.  A
    //brute force method is used.  The 2D cos theta lepton and
    //q2 phase space is probed.

    //Start by declaring a particle at rest.

    //It only makes sense to have a scalar parent.  For now.
    //This should be generalized later.

    //  EvtScalarParticle *scalar_part;
    //  scalar_part=new EvtScalarParticle;

    EvtDiracParticle* dirac_part;
    EvtParticle* root_part;
    dirac_part = new EvtDiracParticle;

    //cludge to avoid generating random numbers!
    //  scalar_part->noLifeTime();
    dirac_part->noLifeTime();

    EvtVector4R p_init;

    p_init.set( EvtPDL::getMass( parent ), 0.0, 0.0, 0.0 );
    //  scalar_part->init(parent,p_init);
    //  root_part=(EvtParticle *)scalar_part;
    //  root_part->set_type(EvtSpinType::SCALAR);

    dirac_part->init( parent, p_init );
    root_part = (EvtParticle*)dirac_part;
    root_part->setDiagonalSpinDensity();

    EvtParticle *daughter, *lep, *trino;

    EvtAmp amp;

    EvtId listdaug[3];
    listdaug[0] = baryon;
    listdaug[1] = lepton;
    listdaug[2] = nudaug;

    amp.init( parent, 3, listdaug );

    root_part->makeDaughters( 3, listdaug );
    daughter = root_part->getDaug( 0 );
    lep = root_part->getDaug( 1 );
    trino = root_part->getDaug( 2 );

    //cludge to avoid generating random numbers!
    daughter->noLifeTime();
    lep->noLifeTime();
    trino->noLifeTime();

    //Initial particle is unpolarized, well it is a scalar so it is
    //trivial
    EvtSpinDensity rho;
    rho.setDiag( root_part->getSpinStates() );

    double mass[3];

    double m = root_part->mass();

    EvtVector4R p4baryon, p4lepton, p4nu, p4w;
    double q2max;

    double q2, elepton, plepton;
    int i, j;
    double erho, prho, costl;

    double maxfoundprob = 0.0;
    double prob = -10.0;
    int massiter;

    for ( massiter = 0; massiter < 3; massiter++ ) {
        mass[0] = EvtPDL::getMass( baryon );
        mass[1] = EvtPDL::getMass( lepton );
        mass[2] = EvtPDL::getMass( nudaug );
        if ( massiter == 1 ) {
            mass[0] = EvtPDL::getMinMass( baryon );
        }
        if ( massiter == 2 ) {
            mass[0] = EvtPDL::getMaxMass( baryon );
        }

        q2max = ( m - mass[0] ) * ( m - mass[0] );

        //loop over q2

        for ( i = 0; i < 25; i++ ) {
            q2 = ( ( i + 0.5 ) * q2max ) / 25.0;

            erho = ( m * m + mass[0] * mass[0] - q2 ) / ( 2.0 * m );

            prho = sqrt( erho * erho - mass[0] * mass[0] );

            p4baryon.set( erho, 0.0, 0.0, -1.0 * prho );
            p4w.set( m - erho, 0.0, 0.0, prho );

            //This is in the W rest frame
            elepton = ( q2 + mass[1] * mass[1] ) / ( 2.0 * sqrt( q2 ) );
            plepton = sqrt( elepton * elepton - mass[1] * mass[1] );

            double probctl[3];

            for ( j = 0; j < 3; j++ ) {
                costl = 0.99 * ( j - 1.0 );

                //These are in the W rest frame. Need to boost out into
                //the B frame.
                p4lepton.set( elepton, 0.0, plepton * sqrt( 1.0 - costl * costl ),
                              plepton * costl );
                p4nu.set( plepton, 0.0,
                          -1.0 * plepton * sqrt( 1.0 - costl * costl ),
                          -1.0 * plepton * costl );

                EvtVector4R boost( ( m - erho ), 0.0, 0.0, 1.0 * prho );
                p4lepton = boostTo( p4lepton, boost );
                p4nu = boostTo( p4nu, boost );

                //Now initialize the daughters...

                daughter->init( baryon, p4baryon );
                lep->init( lepton, p4lepton );
                trino->init( nudaug, p4nu );

                CalcAmp( root_part, amp, FormFactors, r00, r01, r10, r11 );

                //Now find the probability at this q2 and cos theta lepton point
                //and compare to maxfoundprob.

                //Do a little magic to get the probability!!
                prob = rho.normalizedProb( amp.getSpinDensity() );

                probctl[j] = prob;
            }

            //probclt contains prob at ctl=-1,0,1.
            //prob=a+b*ctl+c*ctl^2

            double a = probctl[1];
            double b = 0.5 * ( probctl[2] - probctl[0] );
            double c = 0.5 * ( probctl[2] + probctl[0] ) - probctl[1];

            prob = probctl[0];
            if ( probctl[1] > prob )
                prob = probctl[1];
            if ( probctl[2] > prob )
                prob = probctl[2];

            if ( fabs( c ) > 1e-20 ) {
                double ctlx = -0.5 * b / c;
                if ( fabs( ctlx ) < 1.0 ) {
                    double probtmp = a + b * ctlx + c * ctlx * ctlx;
                    if ( probtmp > prob )
                        prob = probtmp;
                }
            }

            //EvtGenReport(EVTGEN_DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
            //                      << probctl[0]<<" "
            //                      << probctl[1]<<" "
            //                      << probctl[2]<<std::endl;

            if ( prob > maxfoundprob ) {
                maxfoundprob = prob;
            }
        }
        if ( EvtPDL::getWidth( baryon ) <= 0.0 ) {
            //if the particle is narrow dont bother with changing the mass.
            massiter = 4;
        }
    }
    root_part->deleteTree();

    maxfoundprob *= 1.1;
    return maxfoundprob;
}
void EvtSemiLeptonicBaryonAmp::CalcAmp( EvtParticle* parent, EvtAmp& amp,
                                        EvtSemiLeptonicFF* FormFactors,
                                        EvtComplex r00, EvtComplex r01,
                                        EvtComplex r10, EvtComplex r11 )
{
    //  Leptons
    static EvtId EM = EvtPDL::getId( "e-" );
    static EvtId MUM = EvtPDL::getId( "mu-" );
    static EvtId TAUM = EvtPDL::getId( "tau-" );
    //  Anti-Leptons
    static EvtId EP = EvtPDL::getId( "e+" );
    static EvtId MUP = EvtPDL::getId( "mu+" );
    static EvtId TAUP = EvtPDL::getId( "tau+" );

    //  Baryons
    static EvtId LAMCP = EvtPDL::getId( "Lambda_c+" );
    static EvtId LAMC1P = EvtPDL::getId( "Lambda_c(2593)+" );
    static EvtId LAMC2P = EvtPDL::getId( "Lambda_c(2625)+" );
    static EvtId LAMB = EvtPDL::getId( "Lambda_b0" );

    // Anti-Baryons
    static EvtId LAMCM = EvtPDL::getId( "anti-Lambda_c-" );
    static EvtId LAMC1M = EvtPDL::getId( "anti-Lambda_c(2593)-" );
    static EvtId LAMC2M = EvtPDL::getId( "anti-Lambda_c(2625)-" );
    static EvtId LAMBB = EvtPDL::getId( "anti-Lambda_b0" );

    // Set the spin density matrix of the parent baryon
    EvtSpinDensity rho;
    rho.setDim( 2 );
    rho.set( 0, 0, r00 );
    rho.set( 0, 1, r01 );

    rho.set( 1, 0, r10 );
    rho.set( 1, 1, r11 );

    EvtVector4R vector4P = parent->getP4Lab();
    double pmag = vector4P.d3mag();
    double cosTheta = vector4P.get( 3 ) / pmag;

    double theta = acos( cosTheta );
    double phi = atan2( vector4P.get( 2 ), vector4P.get( 1 ) );

    parent->setSpinDensityForwardHelicityBasis( rho, phi, theta, 0.0 );
    //parent->setSpinDensityForward(rho);

    // Set the four momentum of the parent baryon in it's rest frame
    EvtVector4R p4b;
    p4b.set( parent->mass(), 0.0, 0.0, 0.0 );

    // Get the four momentum of the daughter baryon in the parent's rest frame
    EvtVector4R p4daught = parent->getDaug( 0 )->getP4();

    // Add the lepton and neutrino 4 momenta to find q (q^2)
    EvtVector4R q = parent->getDaug( 1 )->getP4() + parent->getDaug( 2 )->getP4();

    double q2 = q.mass2();

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

    double baryonmass = parent->getDaug( 0 )->mass();

    // Handle spin-1/2 daughter baryon Dirac spinor cases
    if ( EvtPDL::getSpinType( parent->getDaug( 0 )->getId() ) ==
         EvtSpinType::DIRAC ) {
        // Set the form factors
        double f1, f2, f3, g1, g2, g3;
        FormFactors->getdiracff( par_num, bar_num, q2, baryonmass, &f1, &f2,
                                 &f3, &g1, &g2, &g3 );

        const double form_fact[6] = {f1, f2, f3, g1, g2, g3};

        EvtVector4C b11, b12, b21, b22, l1, l2;

        //  Lepton Current
        if ( l_num == EM || l_num == MUM || l_num == TAUM ) {
            l1 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 0 ),
                                     parent->getDaug( 2 )->spParentNeutrino() );
            l2 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 1 ),
                                     parent->getDaug( 2 )->spParentNeutrino() );

        } else if ( l_num == EP || l_num == MUP || l_num == TAUP ) {
            l1 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
                                     parent->getDaug( 1 )->spParent( 0 ) );
            l2 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
                                     parent->getDaug( 1 )->spParent( 1 ) );

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

        // Baryon current

        // Flag for particle/anti-particle parent, daughter with same/opp. parity
        // pflag = 0 => particle, same parity parent, daughter
        // pflag = 1 => particle, opp. parity parent, daughter
        // pflag = 2 => anti-particle, same parity parent, daughter
        // pflag = 3 => anti-particle, opp. parity parent, daughter

        int pflag = 0;

        // Handle 1/2+ -> 1/2+ first
        if ( ( par_num == LAMB && bar_num == LAMCP ) ||
             ( par_num == LAMBB && bar_num == LAMCM ) ) {
            // Set particle/anti-particle flag
            if ( bar_num == LAMCP )
                pflag = 0;
            else if ( bar_num == LAMCM )
                pflag = 2;

            b11 = EvtBaryonVACurrent( parent->getDaug( 0 )->spParent( 0 ),
                                      parent->sp( 0 ), p4b, p4daught, form_fact,
                                      pflag );
            b21 = EvtBaryonVACurrent( parent->getDaug( 0 )->spParent( 0 ),
                                      parent->sp( 1 ), p4b, p4daught, form_fact,
                                      pflag );
            b12 = EvtBaryonVACurrent( parent->getDaug( 0 )->spParent( 1 ),
                                      parent->sp( 0 ), p4b, p4daught, form_fact,
                                      pflag );
            b22 = EvtBaryonVACurrent( parent->getDaug( 0 )->spParent( 1 ),
                                      parent->sp( 1 ), p4b, p4daught, form_fact,
                                      pflag );
        }

        // Handle 1/2+ -> 1/2- second
        else if ( ( par_num == LAMB && bar_num == LAMC1P ) ||
                  ( par_num == LAMBB && bar_num == LAMC1M ) ) {
            // Set particle/anti-particle flag
            if ( bar_num == LAMC1P )
                pflag = 1;
            else if ( bar_num == LAMC1M )
                pflag = 3;

            b11 = EvtBaryonVACurrent( ( parent->getDaug( 0 )->spParent( 0 ) ),
                                      ( EvtGammaMatrix::g5() * parent->sp( 0 ) ),
                                      p4b, p4daught, form_fact, pflag );
            b21 = EvtBaryonVACurrent( ( parent->getDaug( 0 )->spParent( 0 ) ),
                                      ( EvtGammaMatrix::g5() * parent->sp( 1 ) ),
                                      p4b, p4daught, form_fact, pflag );
            b12 = EvtBaryonVACurrent( ( parent->getDaug( 0 )->spParent( 1 ) ),
                                      ( EvtGammaMatrix::g5() * parent->sp( 0 ) ),
                                      p4b, p4daught, form_fact, pflag );
            b22 = EvtBaryonVACurrent( ( parent->getDaug( 0 )->spParent( 1 ) ),
                                      ( EvtGammaMatrix::g5() * parent->sp( 1 ) ),
                                      p4b, p4daught, form_fact, pflag );

        }

        else {
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "Dirac semilep. baryon current "
                << "not implemented for this decay sequence." << std::endl;
            ::abort();
        }

        amp.vertex( 0, 0, 0, l1 * b11 );
        amp.vertex( 0, 0, 1, l2 * b11 );

        amp.vertex( 1, 0, 0, l1 * b21 );
        amp.vertex( 1, 0, 1, l2 * b21 );

        amp.vertex( 0, 1, 0, l1 * b12 );
        amp.vertex( 0, 1, 1, l2 * b12 );

        amp.vertex( 1, 1, 0, l1 * b22 );
        amp.vertex( 1, 1, 1, l2 * b22 );

    }

    // Need special handling for the spin-3/2 daughter baryon
    // Rarita-Schwinger spinor cases
    else if ( EvtPDL::getSpinType( parent->getDaug( 0 )->getId() ) ==
              EvtSpinType::RARITASCHWINGER ) {
        // Set the form factors
        double f1, f2, f3, f4, g1, g2, g3, g4;
        FormFactors->getraritaff( par_num, bar_num, q2, baryonmass, &f1, &f2,
                                  &f3, &f4, &g1, &g2, &g3, &g4 );

        const double form_fact[8] = {f1, f2, f3, f4, g1, g2, g3, g4};

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

        EvtVector4C b11, b12, b21, b22, b13, b23, b14, b24, l1, l2;

        //  Lepton Current
        if ( l_num == EM || l_num == MUM || l_num == TAUM ) {
            //  Lepton Current
            l1 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 0 ),
                                     parent->getDaug( 2 )->spParentNeutrino() );
            l2 = EvtLeptonVACurrent( parent->getDaug( 1 )->spParent( 1 ),
                                     parent->getDaug( 2 )->spParentNeutrino() );
        } else if ( l_num == EP || l_num == MUP || l_num == TAUP ) {
            l1 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
                                     parent->getDaug( 1 )->spParent( 0 ) );
            l2 = EvtLeptonVACurrent( parent->getDaug( 2 )->spParentNeutrino(),
                                     parent->getDaug( 1 )->spParent( 1 ) );

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

        //  Baryon Current
        // Declare particle, anti-particle flag, same/opp. parity
        // pflag = 0 => particle
        // pflag = 1 => anti-particle
        int pflag = 0;

        // Handle cases of 1/2+ -> 3/2-
        if ( par_num == LAMB && bar_num == LAMC2P ) {
            // Set flag for particle case
            pflag = 0;
        } else if ( par_num == LAMBB && bar_num == LAMC2M ) {
            // Set flag for anti-particle case
            pflag = 1;
        } else {
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "Rarita-Schwinger semilep. baryon current "
                << "not implemented for this decay sequence." << std::endl;
            ::abort();
        }

        // Baryon current
        b11 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 0 ),
                                        parent->sp( 0 ), p4b, p4daught,
                                        form_fact, pflag );
        b21 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 0 ),
                                        parent->sp( 1 ), p4b, p4daught,
                                        form_fact, pflag );

        b12 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 1 ),
                                        parent->sp( 0 ), p4b, p4daught,
                                        form_fact, pflag );
        b22 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 1 ),
                                        parent->sp( 1 ), p4b, p4daught,
                                        form_fact, pflag );

        b13 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 2 ),
                                        parent->sp( 0 ), p4b, p4daught,
                                        form_fact, pflag );
        b23 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 2 ),
                                        parent->sp( 1 ), p4b, p4daught,
                                        form_fact, pflag );

        b14 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 3 ),
                                        parent->sp( 0 ), p4b, p4daught,
                                        form_fact, pflag );
        b24 = EvtBaryonVARaritaCurrent( parent->getDaug( 0 )->spRSParent( 3 ),
                                        parent->sp( 1 ), p4b, p4daught,
                                        form_fact, pflag );

        amp.vertex( 0, 0, 0, l1 * b11 );
        amp.vertex( 0, 0, 1, l2 * b11 );

        amp.vertex( 1, 0, 0, l1 * b21 );
        amp.vertex( 1, 0, 1, l2 * b21 );

        amp.vertex( 0, 1, 0, l1 * b12 );
        amp.vertex( 0, 1, 1, l2 * b12 );

        amp.vertex( 1, 1, 0, l1 * b22 );
        amp.vertex( 1, 1, 1, l2 * b22 );

        amp.vertex( 0, 2, 0, l1 * b13 );
        amp.vertex( 0, 2, 1, l2 * b13 );

        amp.vertex( 1, 2, 0, l1 * b23 );
        amp.vertex( 1, 2, 1, l2 * b23 );

        amp.vertex( 0, 3, 0, l1 * b14 );
        amp.vertex( 0, 3, 1, l2 * b14 );

        amp.vertex( 1, 3, 0, l1 * b24 );
        amp.vertex( 1, 3, 1, l2 * b24 );
    }
}

EvtVector4C EvtSemiLeptonicBaryonAmp::EvtBaryonVACurrent(
    const EvtDiracSpinor& Bf, const EvtDiracSpinor& Bi, EvtVector4R parent,
    EvtVector4R daught, const double* ff, int pflag )
{
    // flag == 0 => particle, same parity
    // flag == 1 => particle, opposite parity
    // flag == 2 => anti-particle, same parity
    // flag == 3 => anti-particle, opposite parity

    // particle
    EvtComplex cv = EvtComplex( 1.0, 0. );
    EvtComplex ca = EvtComplex( 1.0, 0. );

    EvtComplex cg0 = EvtComplex( 1.0, 0. );
    EvtComplex cg5 = EvtComplex( 1.0, 0. );

    // antiparticle- same parity parent & daughter
    if ( pflag == 2 ) {
        cv = EvtComplex( -1.0, 0. );
        ca = EvtComplex( 1.0, 0. );

        cg0 = EvtComplex( 1.0, 0.0 );
        cg5 = EvtComplex( 0.0, -1.0 );
    }
    // antiparticle- opposite parity parent & daughter
    else if ( pflag == 3 ) {
        cv = EvtComplex( 1.0, 0. );
        ca = EvtComplex( -1.0, 0. );

        cg0 = EvtComplex( 0.0, -1.0 );
        cg5 = EvtComplex( 1.0, 0.0 );
    }

    EvtVector4C t[6];

    // Term 1 = \bar{u}(p',s')*(F_1(q^2)*\gamma_{mu})*u(p,s)
    t[0] = cv * EvtLeptonVCurrent( Bf, Bi );

    // Term 2 = \bar{u}(p',s')*(F_2(q^2)*(p_{mu}/m_{\Lambda_Q}))*u(p,s)
    t[1] = cg0 * EvtLeptonSCurrent( Bf, Bi ) * ( parent / parent.mass() );

    // Term 3 = \bar{u}(p',s')*(F_3(q^2)*(p'_{mu}/m_{\Lambda_q}))*u(p,s)
    t[2] = cg0 * EvtLeptonSCurrent( Bf, Bi ) * ( daught / daught.mass() );

    // Term 4 = \bar{u}(p',s')*(G_1(q^2)*\gamma_{mu}*\gamma_5)*u(p,s)
    t[3] = ca * EvtLeptonACurrent( Bf, Bi );

    // Term 5 =  \bar{u}(p',s')*(G_2(q^2)*(p_{mu}/m_{\Lambda_Q})*\gamma_5)*u(p,s)
    t[4] = cg5 * EvtLeptonPCurrent( Bf, Bi ) * ( parent / parent.mass() );

    // Term 6 = \bar{u}(p',s')*(G_3(q^2)*(p'_{mu}/m_{\Lambda_q})*\gamma_5)*u(p,s)
    t[5] = cg5 * EvtLeptonPCurrent( Bf, Bi ) * ( daught / daught.mass() );

    // Sum the individual terms
    EvtVector4C current = ( ff[0] * t[0] + ff[1] * t[1] + ff[2] * t[2] -
                            ff[3] * t[3] - ff[4] * t[4] - ff[5] * t[5] );

    return current;
}

EvtVector4C EvtSemiLeptonicBaryonAmp::EvtBaryonVARaritaCurrent(
    const EvtRaritaSchwinger& Bf, const EvtDiracSpinor& Bi, EvtVector4R parent,
    EvtVector4R daught, const double* ff, int pflag )
{
    // flag == 0 => particle
    // flag == 1 => anti-particle

    // particle
    EvtComplex cv = EvtComplex( 1.0, 0. );
    EvtComplex ca = EvtComplex( 1.0, 0. );

    EvtComplex cg0 = EvtComplex( 1.0, 0. );
    EvtComplex cg5 = EvtComplex( 1.0, 0. );

    // antiparticle
    if ( pflag == 1 ) {
        cv = EvtComplex( -1.0, 0. );
        ca = EvtComplex( 1.0, 0. );

        cg0 = EvtComplex( 1.0, 0.0 );
        cg5 = EvtComplex( 0.0, -1.0 );
    }

    EvtVector4C t[8];
    EvtTensor4C id;
    id.setdiag( 1.0, 1.0, 1.0, 1.0 );

    EvtDiracSpinor tmp;
    for ( int i = 0; i < 4; i++ ) {
        tmp.set_spinor( i, Bf.getVector( i ) * parent );
    }

    EvtVector4C v1, v2;
    for ( int i = 0; i < 4; i++ ) {
        v1.set( i, EvtLeptonSCurrent( Bf.getSpinor( i ), Bi ) );
        v2.set( i, EvtLeptonPCurrent( Bf.getSpinor( i ), Bi ) );
    }

    // Term 1 = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(F_1(q^2)*\gamma_{mu})*u(p,s)
    t[0] = ( cv / parent.mass() ) * EvtLeptonVCurrent( tmp, Bi );

    // Term 2
    // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(F_2(q^2)*(p_{mu}/m_{\Lambda_Q}))*u(p,s)
    t[1] = ( ( cg0 / parent.mass() ) * EvtLeptonSCurrent( tmp, Bi ) ) *
           ( parent / parent.mass() );

    // Term 3
    // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(F_3(q^2)*(p'_{mu}/m_{\Lambda_q}))*u(p,s)
    t[2] = ( ( cg0 / parent.mass() ) * EvtLeptonSCurrent( tmp, Bi ) ) *
           ( daught / daught.mass() );

    // Term 4 = \bar{u}^{\alpha}(p',s')*(F_4(q^2)*g_{\alpha,\mu})*u(p,s)
    t[3] = cg0 * ( id.cont2( v1 ) );

    // Term 5
    // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})*(G_1(q^2)*\gamma_{mu}*\gamma_5)*u(p,s)
    t[4] = ( ca / parent.mass() ) * EvtLeptonACurrent( tmp, Bi );

    // Term 6
    // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})
    //      *(G_2(q^2)*(p_{mu}/m_{\Lambda_Q})*\gamma_5)*u(p,s)
    t[5] = ( ( cg5 / parent.mass() ) * EvtLeptonPCurrent( tmp, Bi ) ) *
           ( parent / parent.mass() );

    // Term 7
    // = \bar{u}^{\alpha}(p',s')*(p_{\alpha}/m_{\Lambda_Q})
    //      *(G_3(q^2)*(p'_{mu}/m_{\Lambda_q})*\gamma_5)*u(p,s)
    t[6] = ( ( cg5 / parent.mass() ) * EvtLeptonPCurrent( tmp, Bi ) ) *
           ( daught / daught.mass() );

    // Term 8 = \bar{u}^{\alpha}(p',s')*(G_4(q^2)*g_{\alpha,\mu}*\gamma_5))*u(p,s)
    t[7] = cg5 * ( id.cont2( v2 ) );

    // Sum the individual terms
    EvtVector4C current = ( ff[0] * t[0] + ff[1] * t[1] + ff[2] * t[2] +
                            ff[3] * t[3] - ff[4] * t[4] - ff[5] * t[5] -
                            ff[6] * t[6] - ff[7] * t[7] );

    return current;
}