EvtSSDCP.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/EvtSSDCP.hh"
 
#include "EvtGenBase/EvtCPUtil.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtIncoherentMixing.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"
 
#include <stdlib.h>
#include <string>
using std::endl;
 
std::string EvtSSDCP::getName() const
{
    return "SSD_CP";
}
 
EvtDecayBase* EvtSSDCP::clone() const
{
    return new EvtSSDCP;
}
 
void EvtSSDCP::init()
{
    // check that there are 8 or 12 or 14 arguments
 
    checkNArg( 14, 12, 8 );
    checkNDaug( 2 );
 
    EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) );
    EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) );
 
    // Check it is a B0 or B0s
    if ( ( getParentId() != EvtPDL::getId( "B0" ) ) &&
         ( getParentId() != EvtPDL::getId( "anti-B0" ) ) &&
         ( getParentId() != EvtPDL::getId( "B_s0" ) ) &&
         ( getParentId() != EvtPDL::getId( "anti-B_s0" ) ) ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "EvtSSDCP only decays B0 and B0s" << std::endl;
        ::abort();
    }
 
    if ( ( !( d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR ) ) ||
         ( !( ( d2type == EvtSpinType::SCALAR ) ||
              ( d2type == EvtSpinType::VECTOR ) ||
              ( d2type == EvtSpinType::TENSOR ) ) ) ||
         ( !( ( d1type == EvtSpinType::SCALAR ) ||
              ( d1type == EvtSpinType::VECTOR ) ||
              ( d1type == EvtSpinType::TENSOR ) ) ) ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "EvtSSDCP generator expected "
            << "one of the daugters to be a scalar, the other either scalar, vector, or tensor, found:"
            << EvtPDL::name( getDaug( 0 ) ).c_str() << " and "
            << EvtPDL::name( getDaug( 1 ) ).c_str() << endl;
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "Will terminate execution!" << endl;
        ::abort();
    }
 
    m_dm = getArg( 0 ) / EvtConst::c;    //units of 1/mm
 
    m_dgog = getArg( 1 );
 
    m_qoverp = getArg( 2 ) * EvtComplex( cos( getArg( 3 ) ), sin( getArg( 3 ) ) );
    m_poverq = 1.0 / m_qoverp;
 
    m_A_f = getArg( 4 ) * EvtComplex( cos( getArg( 5 ) ), sin( getArg( 5 ) ) );
 
    m_Abar_f = getArg( 6 ) * EvtComplex( cos( getArg( 7 ) ), sin( getArg( 7 ) ) );
 
    if ( getNArg() >= 12 ) {
        m_eigenstate = false;
        m_A_fbar = getArg( 8 ) *
                   EvtComplex( cos( getArg( 9 ) ), sin( getArg( 9 ) ) );
        m_Abar_fbar = getArg( 10 ) *
                      EvtComplex( cos( getArg( 11 ) ), sin( getArg( 11 ) ) );
    } else {
        //I'm somewhat confused about this. For a CP eigenstate set the
        //amplitudes to the same. For a non CP eigenstate CPT invariance
        //is enforced. (ryd)
        if ( ( getDaug( 0 ) == EvtPDL::chargeConj( getDaug( 0 ) ) &&
               getDaug( 1 ) == EvtPDL::chargeConj( getDaug( 1 ) ) ) ||
             ( getDaug( 0 ) == EvtPDL::chargeConj( getDaug( 1 ) ) &&
               getDaug( 1 ) == EvtPDL::chargeConj( getDaug( 0 ) ) ) ) {
            m_eigenstate = true;
        } else {
            m_eigenstate = false;
            m_A_fbar = conj( m_Abar_f );
            m_Abar_fbar = conj( m_A_f );
        }
    }
 
    //FS: new check for z
    if ( getNArg() == 14 ) {    //FS Set m_z parameter if provided else set it 0
        m_z = EvtComplex( getArg( 12 ), getArg( 13 ) );
    } else {
        m_z = EvtComplex( 0.0, 0.0 );
    }
 
    // FS substituted next 2 lines...
 
    //
    //  m_gamma=EvtPDL::getctau(EvtPDL::getId("B0"));  //units of 1/mm
    //m_dgamma=m_gamma*0.5*m_dgog;
    //
    // ...with:
 
    if ( ( getParentId() == EvtPDL::getId( "B0" ) ) ||
         ( getParentId() == EvtPDL::getId( "anti-B0" ) ) ) {
        m_gamma = 1. / EvtPDL::getctau( EvtPDL::getId( "B0" ) );    //gamma/c (1/mm)
    } else {
        m_gamma = 1. / EvtPDL::getctau( EvtPDL::getId( "B_s0" ) );
    }
 
    m_dgamma = m_gamma * m_dgog;    //dgamma/c (1/mm)
 
    if ( verbose() ) {
        EvtGenReport( EVTGEN_INFO, "EvtGen" )
            << "SSD_CP will generate CP/CPT violation:" << endl
            << endl
            << "    " << EvtPDL::name( getParentId() ).c_str() << " --> "
            << EvtPDL::name( getDaug( 0 ) ).c_str() << " + "
            << EvtPDL::name( getDaug( 1 ) ).c_str() << endl
            << endl
            << "using parameters:" << endl
            << endl
            << "  delta(m)  = " << m_dm << " hbar/ps" << endl
            << "dGamma      = " << m_dgamma << " ps-1" << endl
            << "       q/p  = " << m_qoverp << endl
            << "        z  = " << m_z << endl
            << "       tau  = " << 1. / m_gamma << " ps" << endl;
    }
}
 
void EvtSSDCP::initProbMax()
{
    double theProbMax = abs( m_A_f ) * abs( m_A_f ) +
                        abs( m_Abar_f ) * abs( m_Abar_f ) +
                        abs( m_A_fbar ) * abs( m_A_fbar ) +
                        abs( m_Abar_fbar ) * abs( m_Abar_fbar );
 
    if ( m_eigenstate )
        theProbMax *= 2;
 
    EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) );
    EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) );
    if ( d1type == EvtSpinType::TENSOR || d2type == EvtSpinType::TENSOR )
        theProbMax *= 10;
 
    setProbMax( theProbMax );
}
 
void EvtSSDCP::decay( EvtParticle* p )
{
    static const EvtId B0 = EvtPDL::getId( "B0" );
    static const EvtId B0B = EvtPDL::getId( "anti-B0" );
 
    static const EvtId B0s = EvtPDL::getId( "B_s0" );
    static const EvtId B0Bs = EvtPDL::getId( "anti-B_s0" );
 
    double t;
    EvtId other_b;
    EvtId daugs[2];
 
    int flip = 0;
    if ( !m_eigenstate ) {
        if ( EvtRandom::Flat( 0.0, 1.0 ) < 0.5 )
            flip = 1;
    }
 
    if ( !flip ) {
        daugs[0] = getDaug( 0 );
        daugs[1] = getDaug( 1 );
    } else {
        daugs[0] = EvtPDL::chargeConj( getDaug( 0 ) );
        daugs[1] = EvtPDL::chargeConj( getDaug( 1 ) );
    }
 
    EvtParticle* d;
    p->initializePhaseSpace( 2, daugs );
 
    EvtComplex amp;
 
    EvtCPUtil::getInstance()->OtherB( p, t, other_b, 0.5 );    // t is c*Dt (mm)
    //  EvtIncoherentMixing::OtherB( p , t , other_b , 0.5 ) ;
 
    //if (flip) t=-t;
 
    //FS We assume DGamma=GammaLow-GammaHeavy and Dm=mHeavy-mLow
    EvtComplex expH = exp( -EvtComplex( -0.25 * m_dgamma * t, 0.5 * m_dm * t ) );
    EvtComplex expL = exp( EvtComplex( -0.25 * m_dgamma * t, 0.5 * m_dm * t ) );
    //FS Definition of gp and gm
    EvtComplex gp = 0.5 * ( expL + expH );
    EvtComplex gm = 0.5 * ( expL - expH );
    //FS Calculation os sqrt(1-z^2)
    EvtComplex sqz = sqrt( abs( 1 - m_z * m_z ) ) *
                     exp( EvtComplex( 0, arg( 1 - m_z * m_z ) / 2 ) );
 
    //EvtComplex BB=0.5*(expL+expH);                  // <B0|B0(t)>
    //EvtComplex barBB=m_qoverp*0.5*(expL-expH);       // <B0bar|B0(t)>
    //EvtComplex BbarB=m_poverq*0.5*(expL-expH);       // <B0|B0bar(t)>
    //EvtComplex barBbarB=BB;                         // <B0bar|B0bar(t)>
    //  FS redefinition of these guys... (See BAD #188 eq.35 for ref.)
    //  q/p is taken as in the BaBar Phys. Book (opposite sign wrt ref.)
    EvtComplex BB = gp + m_z * gm;             // <B0|B0(t)>
    EvtComplex barBB = sqz * m_qoverp * gm;    // <B0bar|B0(t)>
    EvtComplex BbarB = sqz * m_poverq * gm;    // <B0|B0bar(t)>
    EvtComplex barBbarB = gp - m_z * gm;       // <B0bar|B0bar(t)>
 
    if ( !flip ) {
        if ( other_b == B0B || other_b == B0Bs ) {
            //at t=0 we have a B0
            //EvtGenReport(EVTGEN_INFO,"EvtGen") << "B0B"<<endl;
            amp = BB * m_A_f + barBB * m_Abar_f;
            //std::cout << "noflip B0B tag:"<<amp<<std::endl;
            //amp=0.0;
        }
        if ( other_b == B0 || other_b == B0s ) {
            //EvtGenReport(EVTGEN_INFO,"EvtGen") << "B0"<<endl;
            amp = BbarB * m_A_f + barBbarB * m_Abar_f;
        }
    } else {
        if ( other_b == B0 || other_b == B0s ) {
            amp = BbarB * m_A_fbar + barBbarB * m_Abar_fbar;
            //std::cout << "flip B0 tag:"<<amp<<std::endl;
            //amp=0.0;
        }
        if ( other_b == B0B || other_b == B0Bs ) {
            amp = BB * m_A_fbar + barBB * m_Abar_fbar;
        }
    }
 
    EvtVector4R p4_parent = p->getP4Restframe();
 
    EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) );
 
    EvtVector4R momv;
    EvtVector4R moms;
 
    if ( d2type == EvtSpinType::SCALAR ) {
        d2type = EvtPDL::getSpinType( getDaug( 0 ) );
        d = p->getDaug( 0 );
        momv = d->getP4();
        moms = p->getDaug( 1 )->getP4();
    } else {
        d = p->getDaug( 1 );
        momv = d->getP4();
        moms = p->getDaug( 0 )->getP4();
    }
 
    if ( d2type == EvtSpinType::SCALAR ) {
        vertex( amp );
    }
 
    if ( d2type == EvtSpinType::VECTOR ) {
        double norm = momv.mass() / ( momv.d3mag() * p->mass() );
 
        //std::cout << amp << " " << norm << " " << p4_parent << d->getP4()<< std::endl;
        //    std::cout << EvtPDL::name(d->getId()) << " " << EvtPDL::name(p->getDaug(0)->getId()) <<
        //  " 1and2 " << EvtPDL::name(p->getDaug(1)->getId()) << std::endl;
        //std::cout << d->eps(0) << std::endl;
        //std::cout << d->epsParent(0) << std::endl;
        vertex( 0, amp * norm * p4_parent * ( d->epsParent( 0 ) ) );
        vertex( 1, amp * norm * p4_parent * ( d->epsParent( 1 ) ) );
        vertex( 2, amp * norm * p4_parent * ( d->epsParent( 2 ) ) );
    }
 
    if ( d2type == EvtSpinType::TENSOR ) {
        double norm = d->mass() * d->mass() /
                      ( p4_parent.mass() * d->getP4().d3mag() *
                        d->getP4().d3mag() );
 
        vertex( 0, amp * norm * d->epsTensorParent( 0 ).cont1( p4_parent ) *
                       p4_parent );
        vertex( 1, amp * norm * d->epsTensorParent( 1 ).cont1( p4_parent ) *
                       p4_parent );
        vertex( 2, amp * norm * d->epsTensorParent( 2 ).cont1( p4_parent ) *
                       p4_parent );
        vertex( 3, amp * norm * d->epsTensorParent( 3 ).cont1( p4_parent ) *
                       p4_parent );
        vertex( 4, amp * norm * d->epsTensorParent( 4 ).cont1( p4_parent ) *
                       p4_parent );
    }
 
    return;
}
 
std::string EvtSSDCP::getParamName( int i )
{
    switch ( i ) {
        case 0:
            return "deltaM";
        case 1:
            return "deltaGammaOverGamma";
        case 2:
            return "qOverP";
        case 3:
            return "qOverPPhase";
        case 4:
            return "Af";
        case 5:
            return "AfPhase";
        case 6:
            return "Abarf";
        case 7:
            return "AbarfPhase";
        case 8:
            return "Afbar";
        case 9:
            return "AfbarPhase";
        case 10:
            return "Abarfbar";
        case 11:
            return "AbarfbarPhase";
        case 12:
            return "Z";
        case 13:
            return "ZPhase";
        default:
            return "";
    }
}
 
std::string EvtSSDCP::getParamDefault( int i )
{
    switch ( i ) {
        case 3:
            return "0.0";
        case 4:
            return "1.0";
        case 5:
            return "0.0";
        case 6:
            return "1.0";
        case 7:
            return "0.0";
        default:
            return "";
    }
}