EvtTauHadnu.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/EvtTauHadnu.hh"
 
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtIdSet.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4C.hh"
 
#include <iostream>
#include <stdlib.h>
#include <string>
 
using namespace std;
 
std::string EvtTauHadnu::getName() const
{
    return "TAUHADNU";
}
 
EvtDecayBase* EvtTauHadnu::clone() const
{
    return new EvtTauHadnu;
}
 
void EvtTauHadnu::init()
{
    // check that there are 0 arguments
 
    checkSpinParent( EvtSpinType::DIRAC );
 
    //the last daughter should be a neutrino
    checkSpinDaughter( getNDaug() - 1, EvtSpinType::NEUTRINO );
 
    int i;
    for ( i = 0; i < ( getNDaug() - 1 ); i++ ) {
        checkSpinDaughter( i, EvtSpinType::SCALAR );
    }
 
    bool validndaug = false;
 
    if ( getNDaug() == 4 ) {
        //pipinu
        validndaug = true;
        checkNArg( 7 );
        m_beta = getArg( 0 );
        m_mRho = getArg( 1 );
        m_gammaRho = getArg( 2 );
        m_mRhopr = getArg( 3 );
        m_gammaRhopr = getArg( 4 );
        m_mA1 = getArg( 5 );
        m_gammaA1 = getArg( 6 );
    }
    if ( getNDaug() == 3 ) {
        //pipinu
        validndaug = true;
        checkNArg( 5 );
        m_beta = getArg( 0 );
        m_mRho = getArg( 1 );
        m_gammaRho = getArg( 2 );
        m_mRhopr = getArg( 3 );
        m_gammaRhopr = getArg( 4 );
    }
    if ( getNDaug() == 2 ) {
        //pipinu
        validndaug = true;
        checkNArg( 0 );
    }
 
    if ( !validndaug ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "Have not yet implemented this final state in TAUHADNUKS model"
            << endl;
        EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Ndaug=" << getNDaug() << endl;
        int id;
        for ( id = 0; id < ( getNDaug() - 1 ); id++ )
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "Daug " << id << " " << EvtPDL::name( getDaug( id ) ).c_str()
                << endl;
    }
}
 
void EvtTauHadnu::initProbMax()
{
    if ( getNDaug() == 2 )
        setProbMax( 90.0 );
    if ( getNDaug() == 3 )
        setProbMax( 2500.0 );
    if ( getNDaug() == 4 )
        setProbMax( 30000.0 );
}
 
void EvtTauHadnu::decay( EvtParticle* p )
{
    static const EvtId TAUM = EvtPDL::getId( "tau-" );
 
    EvtIdSet thePis{ "pi+", "pi-", "pi0" };
    EvtIdSet theKs{ "K+", "K-" };
 
    p->initializePhaseSpace( getNDaug(), getDaugs() );
 
    EvtParticle* nut;
    nut = p->getDaug( getNDaug() - 1 );
    p->getDaug( 0 )->getP4();
 
    //get the leptonic current
    EvtVector4C tau1, tau2;
 
    if ( p->getId() == TAUM ) {
        tau1 = EvtLeptonVACurrent( nut->spParentNeutrino(), p->sp( 0 ) );
        tau2 = EvtLeptonVACurrent( nut->spParentNeutrino(), p->sp( 1 ) );
    } else {
        tau1 = EvtLeptonVACurrent( p->sp( 0 ), nut->spParentNeutrino() );
        tau2 = EvtLeptonVACurrent( p->sp( 1 ), nut->spParentNeutrino() );
    }
 
    EvtVector4C hadCurr;
    bool foundHadCurr = false;
    if ( getNDaug() == 2 ) {
        hadCurr = p->getDaug( 0 )->getP4();
        foundHadCurr = true;
    }
    if ( getNDaug() == 3 ) {
        //pi pi0 nu with rho and rhopr resonance
        if ( thePis.contains( getDaug( 0 ) ) && thePis.contains( getDaug( 1 ) ) ) {
            EvtVector4R q1 = p->getDaug( 0 )->getP4();
            EvtVector4R q2 = p->getDaug( 1 )->getP4();
            double m1 = q1.mass();
            double m2 = q2.mass();
 
            hadCurr = Fpi( ( q1 + q2 ).mass2(), m1, m2 ) * ( q1 - q2 );
 
            foundHadCurr = true;
        }
    }
    if ( getNDaug() == 4 ) {
        if ( thePis.contains( getDaug( 0 ) ) && thePis.contains( getDaug( 1 ) ) &&
             thePis.contains( getDaug( 2 ) ) ) {
            //figure out which is the different charged pi
            //want it to be q3
 
            int diffPi( 0 ), samePi1( 0 ), samePi2( 0 );
            if ( getDaug( 0 ) == getDaug( 1 ) ) {
                diffPi = 2;
                samePi1 = 0;
                samePi2 = 1;
            }
            if ( getDaug( 0 ) == getDaug( 2 ) ) {
                diffPi = 1;
                samePi1 = 0;
                samePi2 = 2;
            }
            if ( getDaug( 1 ) == getDaug( 2 ) ) {
                diffPi = 0;
                samePi1 = 1;
                samePi2 = 2;
            }
 
            EvtVector4R q1 = p->getDaug( samePi1 )->getP4();
            EvtVector4R q2 = p->getDaug( samePi2 )->getP4();
            EvtVector4R q3 = p->getDaug( diffPi )->getP4();
 
            double m1 = q1.mass();
            double m2 = q2.mass();
            double m3 = q3.mass();
 
            EvtVector4R Q = q1 + q2 + q3;
            double Q2 = Q.mass2();
            double mA12 = m_mA1 * m_mA1;
 
            double gammaA1X = m_gammaA1 * gFunc( Q2, samePi1 ) /
                              gFunc( mA12, samePi1 );
 
            EvtComplex denBW_A1( mA12 - Q2, -1. * m_mA1 * gammaA1X );
            EvtComplex BW_A1 = mA12 / denBW_A1;
 
            hadCurr = BW_A1 *
                      ( ( ( q1 - q3 ) - ( Q * ( Q * ( q1 - q3 ) ) / Q2 ) ) *
                            Fpi( ( q1 + q3 ).mass2(), m1, m3 ) +
                        ( ( q2 - q3 ) - ( Q * ( Q * ( q2 - q3 ) ) / Q2 ) ) *
                            Fpi( ( q2 + q3 ).mass2(), m2, m3 ) );
 
            foundHadCurr = true;
        }
    }
 
    if ( !foundHadCurr ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "Have not yet implemented this final state in TAUHADNUKS model"
            << endl;
        EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Ndaug=" << getNDaug() << endl;
        int id;
        for ( id = 0; id < ( getNDaug() - 1 ); id++ )
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "Daug " << id << " " << EvtPDL::name( getDaug( id ) ).c_str()
                << endl;
    }
 
    vertex( 0, tau1 * hadCurr );
    vertex( 1, tau2 * hadCurr );
 
    return;
}
 
double EvtTauHadnu::gFunc( double Q2, int dupD )
{
    double mpi = EvtPDL::getMeanMass( getDaug( dupD ) );
    double mpi2 = pow( mpi, 2. );
    if ( Q2 < pow( m_mRho + mpi, 2. ) ) {
        double arg = Q2 - 9. * mpi2;
        return 4.1 * pow( arg, 3. ) * ( 1. - 3.3 * arg + 5.8 * pow( arg, 2. ) );
    } else
        return Q2 * ( 1.623 + 10.38 / Q2 - 9.32 / pow( Q2, 2. ) +
                      0.65 / pow( Q2, 3. ) );
}
 
EvtComplex EvtTauHadnu::Fpi( double s, double xm1, double xm2 )
{
    EvtComplex BW_rho = BW( s, m_mRho, m_gammaRho, xm1, xm2 );
    EvtComplex BW_rhopr = BW( s, m_mRhopr, m_gammaRhopr, xm1, xm2 );
 
    return ( BW_rho + m_beta * BW_rhopr ) / ( 1. + m_beta );
}
 
EvtComplex EvtTauHadnu::BW( double s, double m, double gamma, double xm1,
                            double xm2 )
{
    double m2 = pow( m, 2. );
 
    if ( s > pow( xm1 + xm2, 2. ) ) {
        double qs = sqrt( fabs( ( s - pow( xm1 + xm2, 2. ) ) *
                                ( s - pow( xm1 - xm2, 2. ) ) ) ) /
                    sqrt( s );
        double qm = sqrt( fabs( ( m2 - pow( xm1 + xm2, 2. ) ) *
                                ( m2 - pow( xm1 - xm2, 2. ) ) ) ) /
                    m;
 
        gamma *= m2 / s * pow( qs / qm, 3. );
    } else
        gamma = 0.;
 
    EvtComplex denBW( m2 - s, -1. * sqrt( s ) * gamma );
 
    return m2 / denBW;
}