EvtBsMuMuKK.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/EvtBsMuMuKK.hh"
 
#include "EvtGenBase/EvtCPUtil.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector3R.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtdFunction.hh"
 
const double pi = EvtConst::pi;
const EvtComplex I = EvtComplex( 0.0, 1.0 );
const double sq2 = sqrt( 2.0 );
 
std::string EvtBsMuMuKK::getName() const
{
    return "BS_MUMUKK";
}
 
EvtDecayBase* EvtBsMuMuKK::clone() const
{
    return new EvtBsMuMuKK;
}
 
void EvtBsMuMuKK::init()
{
    // DecFile parameters
    checkNArg( 37 );
 
    // Non-resonant S wave
    m_f_S_NR = getArg( 0 );
    m_delta_S_NR = getArg( 1 );
    m_phis_S_NR = getArg( 2 );
    m_lambda_S_NR_abs = getArg( 3 );
 
    // f0 (S wave)
    m_f_f0 = getArg( 4 );
    m_delta_f0 = getArg( 5 );
    m_phis_f0 = getArg( 6 );
    m_lambda_f0_abs = getArg( 7 );
 
    // phi (P wave)
    m_f_phi = getArg( 8 );
    m_f_phi_0 = getArg( 9 );
    m_delta_phi_0 = getArg( 10 );
    m_phis_phi_0 = getArg( 11 );
    m_lambda_phi_0_abs = getArg( 12 );
    m_f_phi_perp = getArg( 13 );
    m_delta_phi_perp = pi - getArg( 14 );
    m_phis_phi_perp = getArg( 15 );
    m_lambda_phi_perp_abs = getArg( 16 );
    m_delta_phi_par = pi - getArg( 17 );
    m_phis_phi_par = getArg( 18 );
    m_lambda_phi_par_abs = getArg( 19 );
 
    // f2' (D wave)
    m_f_f2p_0 = getArg( 20 );
    m_delta_f2p_0 = getArg( 21 );
    m_phis_f2p_0 = getArg( 22 );
    m_lambda_f2p_0_abs = getArg( 23 );
    m_f_f2p_perp = getArg( 24 );
    m_delta_f2p_perp = pi - getArg( 25 );
    m_phis_f2p_perp = getArg( 26 );
    m_lambda_f2p_perp_abs = getArg( 27 );
    m_delta_f2p_par = pi - getArg( 28 );
    m_phis_f2p_par = getArg( 29 );
    m_lambda_f2p_par_abs = getArg( 30 );
 
    // Time dependence
    m_Gamma = getArg( 31 );
    m_deltaGamma = getArg( 32 );
    m_deltaMs = getArg( 33 );
 
    // mKK window
    m_Mf0 = getArg( 34 );
    m_kin_lower_limit = getArg( 35 );    // the minimum is approx 2.03*MKp
    m_kin_upper_limit = getArg( 36 );
 
    // PDG masses
    m_MBs = EvtPDL::getMass( EvtPDL::getId( "B_s0" ) );
    m_MJpsi = EvtPDL::getMeanMass( EvtPDL::getId( "J/psi" ) );
    m_Mphi = EvtPDL::getMeanMass( EvtPDL::getId( "phi" ) );
    m_Mf2p = EvtPDL::getMeanMass( EvtPDL::getId( "f'_2" ) );
    m_MKp = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
    m_MKm = EvtPDL::getMass( EvtPDL::getId( "K-" ) );
    m_MK0 = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
    m_Mpip = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
    m_Mpi0 = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
    m_Mmu = EvtPDL::getMass( EvtPDL::getId( "mu+" ) );
 
    double MBsSq = m_MBs * m_MBs;
 
    // Amplitudes and other time parameters
    m_A_S_NR = sqrt( m_f_S_NR );
    m_A_f0 = sqrt( m_f_f0 );
 
    m_A_phi_0 = sqrt( m_f_phi_0 * m_f_phi );
    m_A_phi_perp = sqrt( m_f_phi_perp * m_f_phi );
    // Use fabs to make sure subtractions are >= 0, since subtracting 0 from 0 can give -0
    m_A_phi_par = sqrt(
        fabs( m_f_phi - m_A_phi_perp * m_A_phi_perp - m_A_phi_0 * m_A_phi_0 ) );
 
    m_f_f2p = fabs( 1.0 - m_f_S_NR - m_f_f0 - m_f_phi );
    m_A_f2p_0 = sqrt( m_f_f2p_0 * m_f_f2p );
    m_A_f2p_perp = sqrt( m_f_f2p_perp * m_f_f2p );
    m_A_f2p_par = sqrt(
        fabs( m_f_f2p - m_A_f2p_perp * m_A_f2p_perp - m_A_f2p_0 * m_A_f2p_0 ) );
 
    m_ctau = 1.0 / m_Gamma;
    m_Gamma0phi = EvtPDL::getWidth( EvtPDL::getId( "phi" ) );
    m_Gamma0f2p = EvtPDL::getWidth( EvtPDL::getId( "f'_2" ) );
 
    m_kin_middle = 0.5 * ( m_kin_upper_limit + m_kin_lower_limit );
 
    m_int_const_NR = sqrt( Integral( 1.0, 1.0, 0, 1, 1.0, m_kin_lower_limit,
                                     m_kin_upper_limit, 0 ) );
 
    m_int_Flatte_f0 = sqrt( Integral( 1.0, m_Mf0, 0, 1, 1.0, m_kin_lower_limit,
                                      m_kin_upper_limit, 1 ) );
 
    m_p30Kp_mid_CMS = sqrt( ( pow( m_kin_middle, 2 ) - pow( m_MKp + m_MKm, 2 ) ) *
                            ( pow( m_kin_middle, 2 ) - pow( m_MKp - m_MKm, 2 ) ) ) /
                      ( 2.0 * m_kin_middle );
 
    m_p30Kp_ll_CMS =
        sqrt( ( pow( m_kin_lower_limit, 2 ) - pow( m_MKp + m_MKm, 2 ) ) *
              ( pow( m_kin_lower_limit, 2 ) - pow( m_MKp - m_MKm, 2 ) ) ) /
        ( 2.0 * m_kin_lower_limit );
 
    m_p30Kp_phi_CMS = sqrt( ( m_Mphi * m_Mphi - pow( m_MKp + m_MKm, 2 ) ) *
                            ( m_Mphi * m_Mphi - pow( m_MKp - m_MKm, 2 ) ) ) /
                      ( 2.0 * m_Mphi );
 
    m_p30Kp_f2p_CMS = sqrt( ( m_Mf2p * m_Mf2p - pow( m_MKp + m_MKm, 2 ) ) *
                            ( m_Mf2p * m_Mf2p - pow( m_MKp - m_MKm, 2 ) ) ) /
                      ( 2.0 * m_Mf2p );
 
    m_p30Jpsi_mid_CMS = sqrt( ( MBsSq - pow( m_kin_middle + m_MJpsi, 2 ) ) *
                              ( MBsSq - pow( m_kin_middle - m_MJpsi, 2 ) ) ) /
                        ( 2.0 * m_MBs );
 
    m_p30Jpsi_ll_CMS = sqrt( ( MBsSq - pow( m_kin_lower_limit + m_MJpsi, 2 ) ) *
                             ( MBsSq - pow( m_kin_lower_limit - m_MJpsi, 2 ) ) ) /
                       ( 2.0 * m_MBs );
 
    m_p30Jpsi_phi_CMS = sqrt( ( MBsSq - pow( m_Mphi + m_MJpsi, 2 ) ) *
                              ( MBsSq - pow( m_Mphi - m_MJpsi, 2 ) ) ) /
                        ( 2.0 * m_MBs );
 
    m_p30Jpsi_f2p_CMS = sqrt( ( MBsSq - pow( m_Mf2p + m_MJpsi, 2 ) ) *
                              ( MBsSq - pow( m_Mf2p - m_MJpsi, 2 ) ) ) /
                        ( 2.0 * m_MBs );
 
    m_int_BW_phi = sqrt( Integral( m_Gamma0phi, m_Mphi, 1, 0, m_p30Kp_phi_CMS,
                                   m_kin_lower_limit, m_kin_upper_limit, 2 ) );
 
    m_int_BW_f2p = sqrt( Integral( m_Gamma0f2p, m_Mf2p, 2, 1, m_p30Kp_f2p_CMS,
                                   m_kin_lower_limit, m_kin_upper_limit, 2 ) );
 
    // 4 daughters
    checkNDaug( 4 );
 
    // Spin-0 parent
    checkSpinParent( EvtSpinType::SCALAR );    // B_s0 (anti-B_s0)
 
    // Daughters
    checkSpinDaughter( 0, EvtSpinType::DIRAC );     // mu+  (mu-)
    checkSpinDaughter( 1, EvtSpinType::DIRAC );     // mu-  (mu+)
    checkSpinDaughter( 2, EvtSpinType::SCALAR );    // K+   (K-)
    checkSpinDaughter( 3, EvtSpinType::SCALAR );    // K-   (K+)
 
    // B_s0 parent (Parent must be B_s0 or anti-B_s0)
    const EvtId p = getParentId();
    if ( p != EvtPDL::getId( "B_s0" ) && p != EvtPDL::getId( "anti-B_s0" ) ) {
        assert( 0 );
    }
 
    // Daughter types and ordering (should be mu+-, mu-+, K+-, K-+)
    const EvtId d1 = getDaug( 0 );
    const EvtId d2 = getDaug( 1 );
    const EvtId d3 = getDaug( 2 );
    const EvtId d4 = getDaug( 3 );
    if ( !( ( d1 == EvtPDL::getId( "mu+" ) || d1 == EvtPDL::getId( "mu-" ) ) &&
            ( d2 == EvtPDL::getId( "mu-" ) || d2 == EvtPDL::getId( "mu+" ) ) &&
            ( d3 == EvtPDL::getId( "K+" ) || d3 == EvtPDL::getId( "K-" ) ) &&
            ( d4 == EvtPDL::getId( "K-" ) || d4 == EvtPDL::getId( "K+" ) ) ) ) {
        assert( 0 );
    }
}
 
// Get ProbMax
void EvtBsMuMuKK::initProbMax()
{
    const EvtComplex term11 = sqrt( m_p30Jpsi_f2p_CMS * m_p30Kp_f2p_CMS );
 
    const EvtComplex term12 =
        X_J( 2, m_p30Kp_f2p_CMS, 0 ) * X_J( 1, m_p30Jpsi_f2p_CMS, 1 ) *
        m_p30Kp_f2p_CMS * m_p30Kp_f2p_CMS * m_p30Jpsi_f2p_CMS *
        ( m_A_f2p_0 + 0.3 * m_A_f2p_perp + 0.3 * m_A_f2p_par );
 
    const EvtComplex term13 = m_f_f2p *
                              Breit_Wigner( m_Gamma0f2p, m_Mf2p, m_Mf2p, 2,
                                            m_p30Kp_f2p_CMS, m_p30Kp_f2p_CMS ) /
                              m_int_BW_f2p;
 
    const EvtComplex term21 = sqrt( m_p30Jpsi_phi_CMS * m_p30Kp_phi_CMS );
 
    const EvtComplex term22 = X_J( 1, m_p30Kp_phi_CMS, 0 ) * m_p30Kp_phi_CMS *
                              ( 0.65 * m_A_phi_0 + 0.6 * m_A_phi_perp +
                                0.6 * m_A_phi_par );
 
    const EvtComplex term23 = m_f_phi *
                              Breit_Wigner( m_Gamma0phi, m_Mphi, m_Mphi, 1,
                                            m_p30Kp_phi_CMS, m_p30Kp_phi_CMS ) /
                              m_int_BW_phi;
 
    const EvtComplex term31 = sqrt( m_p30Jpsi_ll_CMS * m_p30Kp_ll_CMS );
 
    const EvtComplex term32 = X_J( 1, m_p30Jpsi_ll_CMS, 1 ) * m_p30Jpsi_ll_CMS;
 
    const EvtComplex term33 = m_f_f0 * Flatte( m_Mf0, m_kin_lower_limit ) /
                              m_int_Flatte_f0;
 
    const EvtComplex term41 = sqrt( m_p30Jpsi_mid_CMS * m_p30Kp_mid_CMS );
 
    const EvtComplex term42 = X_J( 1, m_p30Jpsi_mid_CMS, 1 ) * m_p30Jpsi_mid_CMS;
 
    const EvtComplex term43 = 1.2 * m_f_S_NR / m_int_const_NR;
 
    const EvtComplex hm = term11 * term12 * term13 + term21 * term22 * term23 +
                          term31 * term32 * term33 + term41 * term42 * term43;
 
    // Increase by 10%
    setProbMax( 0.5 * abs2( hm ) * 1.1 );
}
 
// Decay function
void EvtBsMuMuKK::decay( EvtParticle* p )
{
    EvtId other_b;
    double time( 0.0 );
    EvtCPUtil::getInstance()->OtherB( p, time, other_b );
    time = -log( EvtRandom::Flat() ) *
           m_ctau;    // This overrules the ctau made in OtherB
 
    if ( EvtCPUtil::getInstance()->isBsMixed( p ) ) {
        p->getParent()->setLifetime( time * EvtConst::c / 1e12 );    // units: mm
    } else {
        p->setLifetime( time * EvtConst::c / 1e12 );    // units: mm
    }
 
    double DGtime = 0.25 * m_deltaGamma * time;
    double DMtime = 0.5 * m_deltaMs * time;
    double mt = exp( -DGtime );
    double pt = exp( +DGtime );
    double cDMt = cos( DMtime );
    double sDMt = sin( DMtime );
    EvtComplex termplus = EvtComplex( cDMt, sDMt );
    EvtComplex terminus = EvtComplex( cDMt, -sDMt );
 
    EvtComplex gplus = 0.5 * ( mt * termplus + pt * terminus );
    EvtComplex gminus = 0.5 * ( mt * termplus - pt * terminus );
 
    EvtId BSB = EvtPDL::getId( "anti-B_s0" );
 
    // Flavour: first assume B_s0, otherwise choose anti-B_s0
    int q( 1 );
    if ( other_b == BSB ) {
        q = -1;
    }
    p->setAttribute( "q", q );
 
    // Amplitudes
    EvtComplex a_S_NR = AmpTime( q, gplus, gminus, m_delta_S_NR,
                                 m_lambda_S_NR_abs, m_A_S_NR, m_phis_S_NR, -1 );
 
    EvtComplex a_f0 = AmpTime( q, gplus, gminus, m_delta_f0, m_lambda_f0_abs,
                               m_A_f0, m_phis_f0, -1 );
 
    EvtComplex a0_phi = AmpTime( q, gplus, gminus, m_delta_phi_0,
                                 m_lambda_phi_0_abs, m_A_phi_0, m_phis_phi_0, 1 );
 
    EvtComplex aperp_phi = AmpTime( q, gplus, gminus, m_delta_phi_perp,
                                    m_lambda_phi_perp_abs, m_A_phi_perp,
                                    m_phis_phi_perp, -1 );
 
    EvtComplex apar_phi = AmpTime( q, gplus, gminus, m_delta_phi_par,
                                   m_lambda_phi_par_abs, m_A_phi_par,
                                   m_phis_phi_par, 1 );
 
    EvtComplex a0_f2p = AmpTime( q, gplus, gminus, m_delta_f2p_0,
                                 m_lambda_f2p_0_abs, m_A_f2p_0, m_phis_f2p_0,
                                 -1 );
 
    EvtComplex aperp_f2p = AmpTime( q, gplus, gminus, m_delta_f2p_perp,
                                    m_lambda_f2p_perp_abs, m_A_f2p_perp,
                                    m_phis_f2p_perp, 1 );
 
    EvtComplex apar_f2p = AmpTime( q, gplus, gminus, m_delta_f2p_par,
                                   m_lambda_f2p_par_abs, m_A_f2p_par,
                                   m_phis_f2p_par, -1 );
 
    // Generate 4-momenta
    double mKK = EvtRandom::Flat( m_kin_lower_limit, m_kin_upper_limit );
    double mass[10] = { m_MJpsi, mKK, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
    double Kmass[10] = { m_MKp, m_MKm, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
    double muMass[10] = { m_Mmu, m_Mmu, 0.0, 0.0, 0.0,
                          0.0,   0.0,   0.0, 0.0, 0.0 };
 
    EvtVector4R mypV[2], mypK[2], mypmu[2];
    EvtGenKine::PhaseSpace( 2, mass, mypV, m_MBs );
    EvtGenKine::PhaseSpace( 2, Kmass, mypK, mKK );
    EvtGenKine::PhaseSpace( 2, muMass, mypmu, m_MJpsi );
 
    EvtVector4R p4mup = boostTo( mypmu[0], mypV[0] );
    EvtVector4R p4mum = boostTo( mypmu[1], mypV[0] );
    EvtVector4R p4Kp = boostTo( mypK[0], mypV[1] );
    EvtVector4R p4Km = boostTo( mypK[1], mypV[1] );
 
    p->makeDaughters( getNDaug(), getDaugs() );
 
    EvtParticle* thisparticle;
    EvtParticle *muplus, *muminus, *Kplus, *Kminus;
 
    // Check particle ID
    for ( int k = 0; k <= 3; k++ ) {
        thisparticle = p->getDaug( k );
        EvtId pId = thisparticle->getId();
 
        if ( pId == EvtPDL::getId( "mu+" ) ) {
            muplus = thisparticle;
            muplus->init( getDaug( k ), p4mup );
 
        } else if ( pId == EvtPDL::getId( "mu-" ) ) {
            muminus = thisparticle;
            muminus->init( getDaug( k ), p4mum );
 
        } else if ( pId == EvtPDL::getId( "K+" ) ) {
            Kplus = thisparticle;
            Kplus->init( getDaug( k ), p4Kp );
 
        } else if ( pId == EvtPDL::getId( "K-" ) ) {
            Kminus = thisparticle;
            Kminus->init( getDaug( k ), p4Km );
        }
    }
 
    EvtVector4R p4KK = p4Kp + p4Km;
    EvtVector4R p4mumu = p4mup + p4mum;
    EvtVector4R p4Bs = p4mumu + p4KK;
 
    double p4KK_mass2 = p4KK.mass2();
    double p4KK_mass = p4KK.mass();
    double p4Bs_mass2 = p4Bs.mass2();
    double p4Bs_mass = p4Bs.mass();
 
    // Kp momentum in the KK CMS
    double p3Kp_KK_CMS = sqrt( ( p4KK_mass2 - pow( m_MKp + m_MKm, 2 ) ) *
                               ( p4KK_mass2 - pow( m_MKp - m_MKm, 2 ) ) ) /
                         ( 2.0 * p4KK_mass );
 
    // J/psi momentum in the KK CMS
    double p3Jpsi_KK_CMS = sqrt( ( p4Bs_mass2 - pow( p4KK_mass + m_MJpsi, 2 ) ) *
                                 ( p4Bs_mass2 - pow( p4KK_mass - m_MJpsi, 2 ) ) ) /
                           ( 2.0 * p4Bs_mass );
 
    // Mass lineshapes
 
    // Non-resonant S wave
    EvtComplex P_NR = 1.0 / m_int_const_NR;
 
    // f0 Flatte
    EvtComplex F_f0 = Flatte( m_Mf0, p4KK_mass ) / m_int_Flatte_f0;
 
    // phi Breit Wigner
    EvtComplex BW_phi = Breit_Wigner( m_Gamma0phi, m_Mphi, p4KK_mass, 1,
                                      m_p30Kp_phi_CMS, p3Kp_KK_CMS ) /
                        m_int_BW_phi;
 
    // f2' Breit Wigner
    EvtComplex BW_f2p = Breit_Wigner( m_Gamma0f2p, m_Mf2p, p4KK_mass, 1,
                                      m_p30Kp_f2p_CMS, p3Kp_KK_CMS ) /
                        m_int_BW_f2p;
 
    // Barrier factors: Always taking the lowest Bs L
    double X_KK_0 = 1.0;
    double X_KK_1 = X_J( 1, p3Kp_KK_CMS, 0 );
    double X_KK_2 = X_J( 2, p3Kp_KK_CMS, 0 );
    double X_NR_Jpsi_1 = X_J( 1, p3Jpsi_KK_CMS, 1 );
    double X_f0_Jpsi_1 = X_J( 1, p3Jpsi_KK_CMS, 1 );
    double X_phi_Jpsi_0 = 1.0;
    double X_f2p_Jpsi_1 = X_J( 1, p3Jpsi_KK_CMS, 1 );
 
    // Birth momentum factors: pow(p3(K+),LR)* pow(p3(J/psi),LB)
    double f_PHSP = sqrt( p3Jpsi_KK_CMS * p3Kp_KK_CMS );
    double f_BMF_NR = p3Jpsi_KK_CMS;
    double f_BMF_f0 = p3Jpsi_KK_CMS;
    double f_BMF_phi = p3Kp_KK_CMS;
    double f_BMF_f2p = p3Kp_KK_CMS * p3Kp_KK_CMS * p3Jpsi_KK_CMS;
 
    // Angular distribution and sum over KK states
    double CosK = EvtDecayAngle( p4Bs, p4KK, p4Kp );
    double CosMu = EvtDecayAngle( p4Bs, p4mumu, p4mup );
    double chi = EvtDecayAngleChi( p4Bs, p4mup, p4mum, p4Kp, p4Km );
 
    // Build helicity amplitudes
 
    // phi
    EvtComplex H0_phi = a0_phi;
    EvtComplex Hp_phi = ( apar_phi + aperp_phi ) / sq2;
    EvtComplex Hm_phi = ( apar_phi - aperp_phi ) / sq2;
 
    // f2p
    EvtComplex H0_f2p = a0_f2p;
    EvtComplex Hp_f2p = ( apar_f2p + aperp_f2p ) / sq2;
    EvtComplex Hm_f2p = ( apar_f2p - aperp_f2p ) / sq2;
 
    // muon polarization +1
    EvtComplex swaveangdist1 = AngularDist( 0, 0, 1, CosK, CosMu, chi );
 
    // KK Spin-0 NR
    EvtComplex mp_hS_NR = a_S_NR * swaveangdist1;
    EvtComplex Amp_p_NR = P_NR * X_KK_0 * X_NR_Jpsi_1 * f_BMF_NR * mp_hS_NR;
 
    // KK Spin-0 f0
    EvtComplex mp_h_f0 = a_f0 * swaveangdist1;
    EvtComplex Amp_p_f0 = F_f0 * X_KK_0 * X_f0_Jpsi_1 * f_BMF_f0 * mp_h_f0;
 
    // KK Spin-1
    EvtComplex mp_h0_phi = H0_phi * AngularDist( 1, 0, 1, CosK, CosMu, chi );
    EvtComplex mp_hp_phi = Hp_phi * AngularDist( 1, 1, 1, CosK, CosMu, chi );
    EvtComplex mp_hm_phi = Hm_phi * AngularDist( 1, -1, 1, CosK, CosMu, chi );
    EvtComplex Amp_p_phi = BW_phi * X_KK_1 * X_phi_Jpsi_0 * f_BMF_phi *
                           ( mp_h0_phi + mp_hp_phi + mp_hm_phi );
 
    // KK Spin-2
    EvtComplex mp_h0_f2p = H0_f2p * AngularDist( 2, 0, 1, CosK, CosMu, chi );
    EvtComplex mp_hp_f2p = Hp_f2p * AngularDist( 2, 1, 1, CosK, CosMu, chi );
    EvtComplex mp_hm_f2p = Hm_f2p * AngularDist( 2, -1, 1, CosK, CosMu, chi );
    EvtComplex Amp_p_f2p = BW_f2p * X_KK_2 * X_f2p_Jpsi_1 * f_BMF_f2p *
                           ( mp_h0_f2p + mp_hp_f2p + mp_hm_f2p );
 
    // muon polarization -1
    EvtComplex swaveangdist2 = AngularDist( 0, 0, -1, CosK, CosMu, chi );
 
    // KK Spin-0 NR
    EvtComplex mm_hS_NR = a_S_NR * swaveangdist2;
    EvtComplex Amp_m_NR = P_NR * X_KK_0 * X_NR_Jpsi_1 * f_BMF_NR * mm_hS_NR;
 
    // KK Spin-0
    EvtComplex mm_h_f0 = a_f0 * swaveangdist2;
    EvtComplex Amp_m_f0 = F_f0 * X_KK_0 * X_f0_Jpsi_1 * f_BMF_f0 * mm_h_f0;
 
    // KK Spin-1
    EvtComplex mm_h0_phi = H0_phi * AngularDist( 1, 0, -1, CosK, CosMu, chi );
    EvtComplex mm_hp_phi = Hp_phi * AngularDist( 1, +1, -1, CosK, CosMu, chi );
    EvtComplex mm_hm_phi = Hm_phi * AngularDist( 1, -1, -1, CosK, CosMu, chi );
    EvtComplex Amp_m_phi = BW_phi * X_KK_1 * X_phi_Jpsi_0 * f_BMF_phi *
                           ( mm_h0_phi + mm_hp_phi + mm_hm_phi );
 
    // KK Spin-2
    EvtComplex mm_h0_f2p = H0_f2p * AngularDist( 2, 0, -1, CosK, CosMu, chi );
    EvtComplex mm_hp_f2p = Hp_f2p * AngularDist( 2, 1, -1, CosK, CosMu, chi );
    EvtComplex mm_hm_f2p = Hm_f2p * AngularDist( 2, -1, -1, CosK, CosMu, chi );
    EvtComplex Amp_m_f2p = BW_f2p * X_KK_2 * X_f2p_Jpsi_1 * f_BMF_f2p *
                           ( mm_h0_f2p + mm_hp_f2p + mm_hm_f2p );
 
    // Total amplitudes
    EvtComplex Amp_tot_plus = f_PHSP *
                              ( Amp_p_NR + Amp_p_f0 + Amp_p_phi + Amp_p_f2p );
    EvtComplex Amp_tot_minus = f_PHSP *
                               ( Amp_m_NR + Amp_m_f0 + Amp_m_phi + Amp_m_f2p );
 
    vertex( 0, 0, 0.0 );
    vertex( 0, 1, Amp_tot_plus );
    vertex( 1, 0, Amp_tot_minus );
    vertex( 1, 1, 0.0 );
}
 
// Rho function
EvtComplex EvtBsMuMuKK::GetRho( const double m0, const double m ) const
{
    double rho_sq = 1.0 - ( 4.0 * m0 * m0 / ( m * m ) );
    EvtComplex rho;
 
    if ( rho_sq > 0.0 ) {
        rho = EvtComplex( sqrt( rho_sq ), 0.0 );
    } else {
        rho = EvtComplex( 0.0, sqrt( -rho_sq ) );
    }
 
    return rho;
}
 
// Flatte function
EvtComplex EvtBsMuMuKK::Flatte( const double m0, const double m ) const
{
    double gpipi = 0.167;
    double gKK = 3.05 * gpipi;
 
    EvtComplex term1 = ( 2.0 * GetRho( m_Mpip, m ) + GetRho( m_Mpi0, m ) ) / 3.0;
    EvtComplex term2 = ( GetRho( m_MKp, m ) + GetRho( m_MK0, m ) ) / 2.0;
 
    EvtComplex w = gpipi * term1 + gKK * term2;
 
    EvtComplex Flatte_0 = 1.0 / ( m0 * m0 - m * m - I * m0 * w );
 
    return Flatte_0;
}
 
// Breit-Wigner function
EvtComplex EvtBsMuMuKK::Breit_Wigner( const double Gamma0, const double m0,
                                      const double m, const int J,
                                      const double q0, const double q ) const
{
    double X_J_q0_sq = pow( X_J( J, q0, 0 ), 2 );
    double X_J_q_sq = pow( X_J( J, q, 0 ), 2 );
 
    double Gamma = Gamma0 * pow( q / q0, 2 * J + 1 ) * ( m0 / m ) *
                   ( X_J_q_sq / X_J_q0_sq );
 
    return 1.0 / ( m0 * m0 - m * m - I * m0 * Gamma );
}
 
// Integral
double EvtBsMuMuKK::Integral( const double Gamma0, const double m0, const int JR,
                              const int JB, const double q0, const double M_KK_ll,
                              const double M_KK_ul, const int fcntype ) const
{
    const int bins = 1000;
    const double bin_width = ( M_KK_ul - M_KK_ll ) / static_cast<double>( bins );
    const double sumMKpKm2 = pow( m_MKp + m_MKm, 2 );
    const double diffMKpKm2 = pow( m_MKp - m_MKm, 2 );
    const double MBs2 = pow( m_MBs, 2 );
 
    EvtComplex integral( 0.0, 0.0 );
 
    for ( int i = 0; i < bins; i++ ) {
        const double M_KK_i = M_KK_ll + static_cast<double>( i ) * bin_width;
        const double M_KK_f = M_KK_ll + static_cast<double>( i + 1 ) * bin_width;
        const double M_KK_i_sq = M_KK_i * M_KK_i;
        const double M_KK_f_sq = M_KK_f * M_KK_f;
 
        const double p3Kp_KK_CMS_i = sqrt( ( M_KK_i_sq - sumMKpKm2 ) *
                                           ( M_KK_i_sq - diffMKpKm2 ) ) /
                                     ( 2.0 * M_KK_i );
        const double p3Kp_KK_CMS_f = sqrt( ( M_KK_f_sq - sumMKpKm2 ) *
                                           ( M_KK_f_sq - diffMKpKm2 ) ) /
                                     ( 2.0 * M_KK_f );
 
        const double p3Jpsi_Bs_CMS_i =
            sqrt( ( MBs2 - pow( M_KK_i + m_MJpsi, 2 ) ) *
                  ( MBs2 - pow( M_KK_i - m_MJpsi, 2 ) ) ) /
            ( 2.0 * m_MBs );
        const double p3Jpsi_Bs_CMS_f =
            sqrt( ( MBs2 - pow( M_KK_f + m_MJpsi, 2 ) ) *
                  ( MBs2 - pow( M_KK_f - m_MJpsi, 2 ) ) ) /
            ( 2.0 * m_MBs );
 
        const double f_PHSP_i = sqrt( p3Kp_KK_CMS_i * p3Jpsi_Bs_CMS_i );
        const double f_PHSP_f = sqrt( p3Kp_KK_CMS_f * p3Jpsi_Bs_CMS_f );
 
        const double f_MBF_KK_i = pow( p3Kp_KK_CMS_i, JR );
        const double f_MBF_KK_f = pow( p3Kp_KK_CMS_f, JR );
 
        const double f_MBF_Bs_i = pow( p3Jpsi_Bs_CMS_i, JB );
        const double f_MBF_Bs_f = pow( p3Jpsi_Bs_CMS_f, JB );
 
        const double X_JR_i = X_J( JR, p3Kp_KK_CMS_i, 0 );
        const double X_JR_f = X_J( JR, p3Kp_KK_CMS_f, 0 );
 
        const double X_JB_i = X_J( JB, p3Jpsi_Bs_CMS_i, 1 );
        const double X_JB_f = X_J( JB, p3Jpsi_Bs_CMS_f, 1 );
 
        EvtComplex fcn_i( 1.0, 0.0 ), fcn_f( 1.0, 0.0 );
 
        if ( fcntype == 1 ) {
            fcn_i = Flatte( m0, M_KK_i );
            fcn_f = Flatte( m0, M_KK_f );
 
        } else if ( fcntype == 2 ) {
            fcn_i = Breit_Wigner( Gamma0, m0, M_KK_i, JR, q0, p3Kp_KK_CMS_i );
            fcn_f = Breit_Wigner( Gamma0, m0, M_KK_f, JR, q0, p3Kp_KK_CMS_f );
        }
 
        const EvtComplex a_i = f_PHSP_i * f_MBF_KK_i * f_MBF_Bs_i * X_JR_i *
                               X_JB_i * fcn_i;
        const EvtComplex a_st_i = conj( a_i );
        const EvtComplex a_f = f_PHSP_f * f_MBF_KK_f * f_MBF_Bs_f * X_JR_f *
                               X_JB_f * fcn_f;
        const EvtComplex a_st_f = conj( a_f );
 
        integral += 0.5 * bin_width * ( a_i * a_st_i + a_f * a_st_f );
    }
 
    return sqrt( abs2( integral ) );
}
 
// Blatt-Weisskopf barrier factors
double EvtBsMuMuKK::X_J( const int J, const double q, const int isB ) const
{
    double r_BW = 1.0;
 
    if ( isB == 0 ) {
        r_BW = 1.5;
    } else if ( isB == 1 ) {
        r_BW = 5.0;
    }
 
    double zsq = pow( r_BW * q, 2 );
 
    double X_J( 1.0 );
 
    if ( J == 1 ) {
        X_J = sqrt( 1.0 / ( 1.0 + zsq ) );
    } else if ( J == 2 ) {
        X_J = sqrt( 1.0 / ( zsq * zsq + 3.0 * zsq + 9.0 ) );
    }
 
    return X_J;
}
 
// EvtGen d matrix: Input is 2J instead of J etc
double EvtBsMuMuKK::Wignerd( const int J, const int l, const int alpha,
                             const double theta ) const
{
    return EvtdFunction::d( 2 * J, 2 * l, 2 * alpha, theta );
}
 
// J spin of KK, l spin proj of J/psi, alpha dimuon spin
EvtComplex EvtBsMuMuKK::AngularDist( const int J, const int l, const int alpha,
                                     const double cK, const double cL,
                                     const double chi ) const
{
    double thetaL = acos( cL );
    double thetaK = acos( cK );
 
    EvtComplex out = 0.5 * sqrt( ( 2 * J + 1 ) / pi ) *
                     exp( EvtComplex( 0, -l * chi ) );
 
    out *= Wignerd( 1, l, alpha, thetaL ) * Wignerd( J, -l, 0, thetaK );
 
    return out;
}
 
// Time-dependent amplitude calculation
EvtComplex EvtBsMuMuKK::AmpTime( const int q, const EvtComplex& gplus,
                                 const EvtComplex& gminus, const double delta,
                                 const double lambda_abs, const double Amp,
                                 const double phis, const int eta ) const
{
    EvtComplex amp_time = Amp * EvtComplex( cos( -delta ), sin( -delta ) );
    double qphis = q * phis;
    amp_time *= ( gplus + eta * pow( lambda_abs, -1.0 * q ) *
                              EvtComplex( cos( qphis ), sin( qphis ) ) * gminus );
 
    if ( q == 1 ) {
        amp_time *= eta;
    }
 
    return amp_time;
}