EvtBTo4piCP.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/EvtBTo4piCP.hh"
 
#include "EvtGenBase/EvtCPUtil.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4R.hh"
 
#include <stdlib.h>
#include <string>
 
EvtComplex EvtAmpA2( const EvtVector4R& p4pi1, const EvtVector4R& p4pi2,
                     const EvtVector4R& p4pi3, const EvtVector4R& p4pi4 )
{
    //added by Lange Jan4,2000
    static const EvtId A2M = EvtPDL::getId( "a_2-" );
    static const EvtId RHO0 = EvtPDL::getId( "rho0" );
 
    EvtVector4R p4a2, p4rho, p4b;
 
    p4rho = p4pi1 + p4pi2;
 
    p4a2 = p4rho + p4pi3;
 
    p4b = p4a2 + p4pi4;
 
    EvtVector4R p4b_a2, p4rho_a2, p4pi1_a2, p4a2_a2;
 
    p4b_a2 = boostTo( p4b, p4a2 );
    p4rho_a2 = boostTo( p4rho, p4a2 );
    p4pi1_a2 = boostTo( p4pi1, p4a2 );
    p4a2_a2 = boostTo( p4a2, p4a2 );
 
    EvtVector4R p4pi1_rho;
 
    p4pi1_rho = boostTo( p4pi1_a2, p4rho_a2 );
 
    EvtVector4R vb, vrho, vpi, t;
 
    vb = p4b_a2 / p4b_a2.d3mag();
    vrho = p4rho_a2 / p4rho_a2.d3mag();
    vpi = p4pi1_rho / p4pi1_rho.d3mag();
 
    t.set( 1.0, 0.0, 0.0, 0.0 );
 
    //  EvtComplex amp_a1,amp_a2;
    EvtComplex amp_a2;
 
    //  double bwm_a1=EvtPDL::getMeanMass(A1M);
    //  double gamma_a1=EvtPDL::getWidth(A1M);
    double bwm_a2 = EvtPDL::getMeanMass( A2M );
    double gamma_a2 = EvtPDL::getWidth( A2M );
    double bwm_rho = EvtPDL::getMeanMass( RHO0 );
    double gamma_rho = EvtPDL::getWidth( RHO0 );
 
    amp_a2 =
        ( sqrt( gamma_a2 / EvtConst::twoPi ) /
          ( ( p4a2 ).mass() - bwm_a2 - EvtComplex( 0.0, 0.5 * gamma_a2 ) ) ) *
        ( sqrt( gamma_rho / EvtConst::twoPi ) /
          ( ( p4rho ).mass() - bwm_rho - EvtComplex( 0.0, 0.5 * gamma_rho ) ) );
 
    return amp_a2 *
           ( vb.get( 1 ) * vrho.get( 1 ) + vb.get( 2 ) * vrho.get( 2 ) +
             vb.get( 3 ) * vrho.get( 3 ) ) *
           ( vpi.get( 1 ) *
                 ( vb.get( 2 ) * vrho.get( 3 ) - vb.get( 3 ) * vrho.get( 2 ) ) +
             vpi.get( 2 ) *
                 ( vb.get( 3 ) * vrho.get( 1 ) - vb.get( 1 ) * vrho.get( 3 ) ) +
             vpi.get( 3 ) *
                 ( vb.get( 1 ) * vrho.get( 2 ) - vb.get( 2 ) * vrho.get( 1 ) ) );
}
 
EvtComplex EvtAmpA1( const EvtVector4R& p4pi1, const EvtVector4R& p4pi2,
                     const EvtVector4R& p4pi3, const EvtVector4R& p4pi4 )
{
    //added by Lange Jan4,2000
    static const EvtId A1M = EvtPDL::getId( "a_1-" );
    static const EvtId RHO0 = EvtPDL::getId( "rho0" );
 
    EvtVector4R p4a1, p4rho, p4b;
 
    p4rho = p4pi1 + p4pi2;
 
    p4a1 = p4rho + p4pi3;
 
    p4b = p4a1 + p4pi4;
 
    EvtVector4R p4b_a1, p4rho_a1, p4pi1_a1, p4a1_a1;
 
    p4b_a1 = boostTo( p4b, p4a1 );
    p4rho_a1 = boostTo( p4rho, p4a1 );
    p4pi1_a1 = boostTo( p4pi1, p4a1 );
    p4a1_a1 = boostTo( p4a1, p4a1 );
 
    EvtVector4R p4pi1_rho;
 
    p4pi1_rho = boostTo( p4pi1_a1, p4rho_a1 );
 
    EvtVector4R vb, vrho, vpi, t;
 
    vb = p4b_a1 / p4b_a1.d3mag();
    vrho = p4rho_a1 / p4rho_a1.d3mag();
    vpi = p4pi1_rho / p4pi1_rho.d3mag();
 
    t.set( 1.0, 0.0, 0.0, 0.0 );
 
    EvtComplex amp_a1;
 
    double bwm_a1 = EvtPDL::getMeanMass( A1M );
    double gamma_a1 = EvtPDL::getWidth( A1M );
    //  double bwm_a2=EvtPDL::getMeanMass(A2M);
    //  double gamma_a2=EvtPDL::getWidth(A2M);
    double bwm_rho = EvtPDL::getMeanMass( RHO0 );
    double gamma_rho = EvtPDL::getWidth( RHO0 );
 
    amp_a1 =
        ( sqrt( gamma_a1 / EvtConst::twoPi ) /
          ( ( p4a1 ).mass() - bwm_a1 - EvtComplex( 0.0, 0.5 * gamma_a1 ) ) ) *
        ( sqrt( gamma_rho / EvtConst::twoPi ) /
          ( ( p4rho ).mass() - bwm_rho - EvtComplex( 0.0, 0.5 * gamma_rho ) ) );
 
    return amp_a1 * ( vb.get( 1 ) * vpi.get( 1 ) + vb.get( 2 ) * vpi.get( 2 ) +
                      vb.get( 3 ) * vpi.get( 3 ) );
}
 
std::string EvtBTo4piCP::getName() const
{
    return "BTO4PI_CP";
}
 
EvtBTo4piCP* EvtBTo4piCP::clone() const
{
    return new EvtBTo4piCP;
}
 
void EvtBTo4piCP::initProbMax()
{
    setProbMax( 50.0 );
}
 
void EvtBTo4piCP::init()
{
    // check that there are 18 arguments
    checkNArg( 18 );
    checkNDaug( 4 );
 
    checkSpinParent( EvtSpinType::SCALAR );
 
    checkSpinDaughter( 0, EvtSpinType::SCALAR );
    checkSpinDaughter( 1, EvtSpinType::SCALAR );
    checkSpinDaughter( 2, EvtSpinType::SCALAR );
    checkSpinDaughter( 3, EvtSpinType::SCALAR );
}
 
void EvtBTo4piCP::decay( EvtParticle* p )
{
    //added by Lange Jan4,2000
    static const EvtId B0 = EvtPDL::getId( "B0" );
    static const EvtId B0B = EvtPDL::getId( "anti-B0" );
 
    double t;
    EvtId other_b;
 
    EvtCPUtil::getInstance()->OtherB( p, t, other_b, 0.5 );
 
    p->initializePhaseSpace( getNDaug(), getDaugs() );
    EvtVector4R mom1 = p->getDaug( 0 )->getP4();
    EvtVector4R mom2 = p->getDaug( 1 )->getP4();
    EvtVector4R mom3 = p->getDaug( 2 )->getP4();
    EvtVector4R mom4 = p->getDaug( 3 )->getP4();
 
    //  double alpha=getArg(0);
    //double dm=getArg(1);
 
    EvtComplex amp;
 
    EvtComplex A, Abar;
 
    EvtComplex A_a1p, Abar_a1p, A_a2p, Abar_a2p;
    EvtComplex A_a1m, Abar_a1m, A_a2m, Abar_a2m;
 
    A_a1p = EvtComplex( getArg( 2 ) * cos( getArg( 3 ) ),
                        getArg( 2 ) * sin( getArg( 3 ) ) );
    Abar_a1p = EvtComplex( getArg( 4 ) * cos( getArg( 5 ) ),
                           getArg( 4 ) * sin( getArg( 5 ) ) );
 
    A_a2p = EvtComplex( getArg( 6 ) * cos( getArg( 7 ) ),
                        getArg( 6 ) * sin( getArg( 7 ) ) );
    Abar_a2p = EvtComplex( getArg( 8 ) * cos( getArg( 9 ) ),
                           getArg( 8 ) * sin( getArg( 9 ) ) );
 
    A_a1m = EvtComplex( getArg( 10 ) * cos( getArg( 11 ) ),
                        getArg( 10 ) * sin( getArg( 11 ) ) );
    Abar_a1m = EvtComplex( getArg( 12 ) * cos( getArg( 13 ) ),
                           getArg( 12 ) * sin( getArg( 13 ) ) );
 
    A_a2m = EvtComplex( getArg( 14 ) * cos( getArg( 15 ) ),
                        getArg( 14 ) * sin( getArg( 15 ) ) );
    Abar_a2m = EvtComplex( getArg( 16 ) * cos( getArg( 17 ) ),
                           getArg( 16 ) * sin( getArg( 17 ) ) );
 
    EvtComplex a2p_amp = EvtAmpA2( mom1, mom2, mom3, mom4 ) +
                         EvtAmpA2( mom1, mom4, mom3, mom2 ) +
                         EvtAmpA2( mom3, mom2, mom1, mom4 ) +
                         EvtAmpA2( mom3, mom4, mom1, mom2 );
 
    EvtComplex a2m_amp = EvtAmpA2( mom2, mom3, mom4, mom1 ) +
                         EvtAmpA2( mom2, mom1, mom4, mom3 ) +
                         EvtAmpA2( mom4, mom3, mom2, mom1 ) +
                         EvtAmpA2( mom4, mom1, mom2, mom3 );
 
    EvtComplex a1p_amp = EvtAmpA1( mom1, mom2, mom3, mom4 ) +
                         EvtAmpA1( mom1, mom4, mom3, mom2 ) +
                         EvtAmpA1( mom3, mom2, mom1, mom4 ) +
                         EvtAmpA1( mom3, mom4, mom1, mom2 );
 
    EvtComplex a1m_amp = EvtAmpA1( mom2, mom3, mom4, mom1 ) +
                         EvtAmpA1( mom2, mom1, mom4, mom3 ) +
                         EvtAmpA1( mom4, mom3, mom2, mom1 ) +
                         EvtAmpA1( mom4, mom1, mom2, mom3 );
 
    A = A_a2p * a2p_amp + A_a1p * a1p_amp + A_a2m * a2m_amp + A_a1m * a1m_amp;
    Abar = Abar_a2p * a2p_amp + Abar_a1p * a1p_amp + Abar_a2m * a2m_amp +
           Abar_a1m * a1m_amp;
 
    if ( other_b == B0B ) {
        amp = A * cos( getArg( 1 ) * t / ( 2 * EvtConst::c ) ) +
              EvtComplex( cos( -2.0 * getArg( 0 ) ), sin( -2.0 * getArg( 0 ) ) ) *
                  getArg( 2 ) * EvtComplex( 0.0, 1.0 ) * Abar *
                  sin( getArg( 1 ) * t / ( 2 * EvtConst::c ) );
    }
    if ( other_b == B0 ) {
        amp = A * EvtComplex( cos( 2.0 * getArg( 0 ) ), sin( 2.0 * getArg( 0 ) ) ) *
                  EvtComplex( 0.0, 1.0 ) *
                  sin( getArg( 1 ) * t / ( 2 * EvtConst::c ) ) +
              getArg( 2 ) * Abar * cos( getArg( 1 ) * t / ( 2 * EvtConst::c ) );
    }
 
    vertex( amp );
 
    return;
}