EvtBtoKD3P.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/EvtBtoKD3P.hh"
 
#include "EvtGenBase/EvtCyclic3.hh"
#include "EvtGenBase/EvtDalitzPoint.hh"
#include "EvtGenBase/EvtDecayTable.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
 
#include "EvtGenModels/EvtPto3P.hh"
 
#include <assert.h>
using std::endl;
 
//------------------------------------------------------------------
EvtDecayBase* EvtBtoKD3P::clone() const
{
    return new EvtBtoKD3P();
}
 
//------------------------------------------------------------------
std::string EvtBtoKD3P::getName() const
{
    return "BTOKD3P";
}
 
//------------------------------------------------------------------
void EvtBtoKD3P::init()
{
    checkNArg( 2 );     // r, phase
    checkNDaug( 3 );    // K, D0(allowed), D0(suppressed).
                        // The last two daughters are really one particle
 
    // check that the mother and all daughters are scalars:
    checkSpinParent( EvtSpinType::SCALAR );
    checkSpinDaughter( 0, EvtSpinType::SCALAR );
    checkSpinDaughter( 1, EvtSpinType::SCALAR );
    checkSpinDaughter( 2, EvtSpinType::SCALAR );
 
    // Check that the B dtr types are K D D:
 
    // get the parameters:
    m_r = getArg( 0 );
    double phase = getArg( 1 );
    m_exp = EvtComplex( cos( phase ), sin( phase ) );
}
 
//------------------------------------------------------------------
void EvtBtoKD3P::initProbMax()
{
    setProbMax( 1 );    // this is later changed in decay()
}
 
//------------------------------------------------------------------
void EvtBtoKD3P::decay( EvtParticle* p )
{
    // tell the subclass that we decay the daughter:
    m_daugsDecayedByParentModel = true;
 
    // the K is the 1st daughter of the B EvtParticle.
    // The decay mode of the allowed D (the one produced in b->c decay) is 2nd
    // The decay mode of the suppressed D (the one produced in b->u decay) is 3rd
    const int KIND = 0;
    const int D1IND = 1;
    const int D2IND = 2;
 
    // generate kinematics of daughters (K and D):
    EvtId tempDaug[2] = { getDaug( KIND ), getDaug( D1IND ) };
    p->initializePhaseSpace( 2, tempDaug );
 
    // Get the D daughter particle and the decay models of the allowed
    // and suppressed D modes:
    EvtParticle* theD = p->getDaug( D1IND );
    EvtPto3P* model1 = dynamic_cast<EvtPto3P*>(
        EvtDecayTable::getInstance().getDecayFunc( theD ) );
 
    // For the suppressed mode, re-initialize theD as the suppressed D alias.
    // First set the id, then re-initialize (since it matches the expected id)
    theD->setId( getDaug( D2IND ) );
    theD->init( getDaug( D2IND ), theD->getP4() );
    EvtPto3P* model2 = dynamic_cast<EvtPto3P*>(
        EvtDecayTable::getInstance().getDecayFunc( theD ) );
 
    // on the first call:
    if ( false == m_decayedOnce ) {
        m_decayedOnce = true;
 
        // store the D decay model pointers:
        m_model1 = model1;
        m_model2 = model2;
 
        // check the decay models of the first 2 daughters and that they
        // have the same final states:
        std::string name1 = model1->getName();
        std::string name2 = model2->getName();
 
        if ( name1 != "PTO3P" ) {
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "D daughters of EvtBtoKD3P decay must decay via the \"PTO3P\" model"
                << endl
                << "    but found to decay via " << name1.c_str() << " or "
                << name2.c_str() << ". Will terminate execution!" << endl;
            assert( 0 );
        }
 
        const EvtId* daugs1 = model1->getDaugs();
        const EvtId* daugs2 = model2->getDaugs();
 
        bool idMatch = true;
        int d;
        for ( d = 0; d < 2; ++d ) {
            if ( daugs1[d] != daugs2[d] ) {
                idMatch = false;
            }
        }
        if ( false == idMatch ) {
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "D daughters of EvtBtoKD3P decay must decay to the same final state"
                << endl
                << "   particles in the same order (not CP-conjugate order),"
                << endl
                << "   but they were found to decay to" << endl;
            for ( d = 0; d < model1->getNDaug(); ++d ) {
                EvtGenReport( EVTGEN_ERROR, "" )
                    << "   " << EvtPDL::name( daugs1[d] ).c_str() << " ";
            }
            EvtGenReport( EVTGEN_ERROR, "" ) << endl;
            for ( d = 0; d < model1->getNDaug(); ++d ) {
                EvtGenReport( EVTGEN_ERROR, "" )
                    << "   " << EvtPDL::name( daugs2[d] ).c_str() << " ";
            }
            EvtGenReport( EVTGEN_ERROR, "" )
                << endl
                << ". Will terminate execution!" << endl;
            assert( 0 );
        }
 
        // estimate the probmax. Need to know the probmax's of the 2
        // models for this:
        setProbMax(
            model1->getProbMax( 0 ) + m_r * m_r * model2->getProbMax( 0 ) +
            2 * m_r * sqrt( model1->getProbMax( 0 ) * model2->getProbMax( 0 ) ) );
 
    }    // end of things to do on the first call
 
    // make sure the models haven't changed since the first call:
    if ( m_model1 != model1 || m_model2 != model2 ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "D daughters of EvtBtoKD3P decay should have only 1 decay modes, "
            << endl
            << "    but a new decay mode was found after the first call" << endl
            << "    Will terminate execution!" << endl;
        assert( 0 );
    }
 
    // Reset the D id to the 1st D
    theD->setId( getDaug( D1IND ) );
    theD->init( getDaug( D1IND ), theD->getP4() );
 
    // get the cover function for each of the models and add them up.
    // They are summed with coefficients 1 because we are willing to
    // take a small inefficiency (~50%) in order to ensure that the
    // cover function is large enough without getting into complications
    // associated with the smallness of m_r:
    EvtPdfSum<EvtDalitzPoint>* pc1 = model1->getPC();
    EvtPdfSum<EvtDalitzPoint>* pc2 = model2->getPC();
    EvtPdfSum<EvtDalitzPoint> pc;
    pc.addTerm( 1.0, *pc1 );
    pc.addTerm( 1.0, *pc2 );
 
    // from this combined cover function, generate the Dalitz point:
    EvtDalitzPoint x = pc.randomPoint();
 
    // get the aptitude for each of the models on this point and add them up:
    EvtComplex amp1 = model1->amplNonCP( x );
    EvtComplex amp2 = model2->amplNonCP( x );
    EvtComplex amp = amp1 + amp2 * m_r * m_exp;
 
    // get the value of the cover function for this point and set the
    // relative amplitude for this decay:
 
    double comp = sqrt( pc.evaluate( x ) );
    vertex( amp / comp );
 
    // Make the daughters of theD:
    bool massTreeOK = theD->generateMassTree();
    if ( massTreeOK == false ) {
        return;
    }
 
    // Now generate the p4's of the daughters of theD:
    std::vector<EvtVector4R> v = model2->initDaughters( x );
 
    if ( v.size() != theD->getNDaug() ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "Number of daughters " << theD->getNDaug() << " != "
            << "Momentum vector size " << v.size() << endl
            << "     Terminating execution." << endl;
        assert( 0 );
    }
 
    // Apply the new p4's to the daughters:
    for ( unsigned int i = 0; i < theD->getNDaug(); ++i ) {
        theD->getDaug( i )->init( model2->getDaugs()[i], v[i] );
    }
}