EvtTVP.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/EvtTVP.hh"
 
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtSpinType.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"
 
#include <cmath>
 
std::string EvtTVP::getName() const
{
    return "TVP";
}
 
EvtDecayBase* EvtTVP::clone() const
{
    return new EvtTVP;
}
 
void EvtTVP::decay( EvtParticle* root )
{
    if ( getNDaug() == 2 ) {
        decay_2body( root );
    } else if ( getNDaug() == 3 ) {
        decay_3body( root );
    }
}
 
void EvtTVP::init()
{
    checkSpinParent( EvtSpinType::TENSOR );
 
    if ( getNDaug() == 2 ) {    // chi -> gamma psi radiative mode
        checkNArg( 0 );
        checkSpinDaughter( 0, EvtSpinType::PHOTON );
        checkSpinDaughter( 1, EvtSpinType::VECTOR );
    } else if ( getNDaug() == 3 ) {    // chi -> psi lepton lepton
        checkNDaug( 3 );
        checkSpinDaughter( 0, EvtSpinType::VECTOR );
        checkSpinDaughter( 1, EvtSpinType::DIRAC );
        checkSpinDaughter( 2, EvtSpinType::DIRAC );
        checkNArg( 1 );
        m_delta = getArg( 0 );
    }
}
 
void EvtTVP::initProbMax()
{
    if ( getNDaug() == 2 ) {
        const EvtId parId = getParentId();
        if ( parId == EvtPDL::getId( "chi_b2" ) ) {
            setProbMax( 15.0 );
        } else {
            setProbMax( 2.5 );
        }
 
    } else if ( getNDaug() == 3 ) {
        double dSq = m_delta * m_delta;
        double denom = dSq - 0.2;
        double ratio( 1.0 );
        if ( fabs( denom ) > 1e-10 ) {
            ratio = dSq / denom;
        }
        double ffCor = ratio * ratio;
 
        const EvtId daugId = getDaug( 1 );
        const EvtId parId = getParentId();
 
        if ( daugId == EvtPDL::getId( "mu+" ) ||
             daugId == EvtPDL::getId( "mu-" ) ) {
            if ( parId == EvtPDL::getId( "chi_c2" ) ) {
                setProbMax( ffCor * 85.0 );    // tested on 1e6 events
            } else if ( parId == EvtPDL::getId( "chi_b2" ) ) {
                setProbMax( ffCor * 750.0 );    // tested on 1e6 events
            }
 
        } else if ( daugId == EvtPDL::getId( "e+" ) ||
                    daugId == EvtPDL::getId( "e-" ) ) {
            if ( parId == EvtPDL::getId( "chi_c2" ) ) {
                setProbMax( ffCor * 3.5e3 );    // tested on 1e5 events
            } else if ( parId == EvtPDL::getId( "chi_b2" ) ) {
                setProbMax( ffCor * 2.6e4 );
            }
        }
    }
}
 
void EvtTVP::decay_2body( EvtParticle* root )
{
    root->initializePhaseSpace( getNDaug(), getDaugs() );
 
    // Photon is the first particle and psi is the second
    // to ensure decay file backwards compatibility
    EvtParticle* photon = root->getDaug( 0 );
    EvtParticle* psi = root->getDaug( 1 );
 
    EvtVector4R p = psi->getP4(),    // psi momentum
        k = photon->getP4();         // Photon momentum
 
    for ( int iPsi = 0; iPsi < 3; iPsi++ ) {
        EvtVector4C epsPsi = psi->epsParent( iPsi ).conj();
 
        for ( int iGamma = 0; iGamma < 2; iGamma++ ) {
            EvtVector4C epsGamma = photon->epsParentPhoton( iGamma ).conj();
 
            for ( int iChi = 0; iChi < 5; iChi++ ) {
                EvtTensor4C epsChi = root->epsTensor( iChi );
 
                // Baranov PRD 85,014034 (2012), Eq 11
                // amp = p^mu epsPsi^a epsChi_{a b} [k_mu epsGamma_b  - k_b epsGamma_mu]
                EvtVector4C eee = epsChi.cont1( epsPsi );
                EvtVector4C vvv = ( p * k ) * eee - ( k * eee ) * p;
                EvtComplex amp = vvv * epsGamma;
                vertex( iChi, iGamma, iPsi, amp );
            }
        }
    }
}
 
void EvtTVP::decay_3body( EvtParticle* root )
{
    root->initializePhaseSpace( getNDaug(), getDaugs() );
    EvtParticle* psi = root->getDaug( 0 );
    EvtParticle* mup = root->getDaug( 1 );
    EvtParticle* mum = root->getDaug( 2 );
 
    EvtVector4R p = psi->getP4(),    // psi momentum
        k1 = mup->getP4(),           // mu+ momentum
        k2 = mum->getP4(),           // mu- momentum
        k = k1 + k2;                 // photon momentum
 
    double kSq = k * k;
 
    // The decay amplitude needs four-vector products. Make sure we have
    // valid values for these, otherwise set the amplitude to zero.
    // We need to set _amp2 (EvtDecayAmp) via the vertex() function call
    // even when the amplitude is zero, otherwise the amplitude from the
    // previous accepted event will be used, potentially leading to biases
 
    // Selection on k^2 to avoid inefficient generation for the electron modes
    bool validAmp( true );
    if ( kSq < 1e-3 ) {
        validAmp = false;
    }
 
    double dSq = m_delta * m_delta;
    double dSqDenom = dSq - kSq;
    if ( fabs( dSqDenom ) < 1e-10 ) {
        validAmp = false;
    }
 
    double factor( 1.0 );
    if ( validAmp ) {
        factor = dSq / ( dSqDenom * kSq );
    }
 
    // Calculate the amplitude terms, looping over the psi and lepton states
    int iPols[4] = { 0, 0, 0, 0 };
 
    for ( int iChi = 0; iChi < 5; iChi++ ) {
        iPols[0] = iChi;
        EvtTensor4C epsChi = root->epsTensor( iChi );
 
        for ( int iPsi = 0; iPsi < 3; iPsi++ ) {
            iPols[1] = iPsi;
            EvtVector4C epsPsi = psi->epsParent( iPsi ).conj();
 
            for ( int iMplus = 0; iMplus < 2; iMplus++ ) {
                iPols[2] = iMplus;
                EvtDiracSpinor spMplus = mup->spParent( iMplus );
 
                for ( int iMminus = 0; iMminus < 2; iMminus++ ) {
                    iPols[3] = iMminus;
                    EvtDiracSpinor spMminus = mum->spParent( iMminus );
                    EvtVector4C epsGamma = EvtLeptonVCurrent( spMplus, spMminus );
 
                    // Based on Baranov PRD 85,014034 (2012), Eq 11
                    // amp = p^mu epsPsi^a epsChi_{a b} [k_mu epsGamma_b  - k_b epsGamma_mu]/k^2
                    EvtVector4C eee = epsChi.cont1( epsPsi );
                    EvtVector4C vvv = ( p * k ) * eee - ( k * eee ) * p;
                    EvtComplex amp( 0.0, 0.0 );
                    if ( validAmp ) {
                        amp = vvv * epsGamma;
                    }
                    amp *= factor;
 
                    // Set the amplitude matrix element using the vertex function
                    vertex( iPols, amp );
                }
            }
        }
    }
}