EvtVubHybrid.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/EvtVubHybrid.hh"
 
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4R.hh"
 
#include "EvtGenModels/EvtPFermi.hh"
#include "EvtGenModels/EvtVubdGamma.hh"
 
#include <fstream>
#include <iostream>
#include <stdlib.h>
#include <string>
using std::cout;
using std::endl;
using std::ifstream;
 
std::string EvtVubHybrid::getName() const
{
    return "VUBHYBRID";
}
 
EvtDecayBase* EvtVubHybrid::clone() const
{
    return new EvtVubHybrid;
}
 
void EvtVubHybrid::init()
{
    // check that there are at least 3 arguments
    if ( getNArg() < EvtVubHybrid::nParameters ) {
        EvtGenReport( EVTGEN_ERROR, "EvtVubHybrid" )
            << "EvtVub generator expected "
            << "at least " << EvtVubHybrid::nParameters
            << " arguments but found: " << getNArg()
            << "\nWill terminate execution!" << endl;
        ::abort();
    } else if ( getNArg() == EvtVubHybrid::nParameters ) {
        EvtGenReport( EVTGEN_WARNING, "EvtVubHybrid" )
            << "EvtVub: generate B -> Xu l nu events "
            << "without using the hybrid reweighting." << endl;
        m_noHybrid = true;
    } else if ( getNArg() < EvtVubHybrid::nParameters + EvtVubHybrid::nVariables ) {
        EvtGenReport( EVTGEN_ERROR, "EvtVubHybrid" )
            << "EvtVub could not read number of bins for "
            << "all variables used in the reweighting\n"
            << "Will terminate execution!" << endl;
        ::abort();
    }
 
    // check that there are 3 daughters
    checkNDaug( 3 );
 
    // read minimum required parameters from decay.dec
    m_mb = getArg( 0 );
    m_a = getArg( 1 );
    m_alphas = getArg( 2 );
 
    // the maximum dGamma*p2 value depends on alpha_s only:
    const double dGMax0 = 3.;
    m_dGMax = 0.21344 + 8.905 * m_alphas;
    if ( m_dGMax < dGMax0 )
        m_dGMax = dGMax0;
 
    // for the Fermi Motion we need a B-Meson mass - but it's not critical
    // to get an exact value; in order to stay in the phase space for
    // B+- and B0 use the smaller mass
 
    static const double mB0 = EvtPDL::getMaxMass( EvtPDL::getId( "B0" ) );
    static const double mBP = EvtPDL::getMaxMass( EvtPDL::getId( "B+" ) );
    static const double mB = ( mB0 < mBP ? mB0 : mBP );
 
    const double xlow = -m_mb;
    const double xhigh = mB - m_mb;
    const int aSize = 10000;
 
    EvtPFermi pFermi( m_a, mB, m_mb );
    // pf is the cumulative distribution normalized to 1.
    m_pf.resize( aSize );
    for ( int i = 0; i < aSize; i++ ) {
        double kplus = xlow + (double)( i + 0.5 ) / ( (double)aSize ) *
                                  ( xhigh - xlow );
        if ( i == 0 )
            m_pf[i] = pFermi.getFPFermi( kplus );
        else
            m_pf[i] = m_pf[i - 1] + pFermi.getFPFermi( kplus );
    }
    for ( size_t index = 0; index < m_pf.size(); index++ ) {
        m_pf[index] /= m_pf[m_pf.size() - 1];
    }
 
    m_dGamma = std::make_unique<EvtVubdGamma>( m_alphas );
 
    if ( m_noHybrid )
        return;    // Without hybrid weighting, nothing else to do
 
    m_bins_mX.resize( abs( (int)getArg( 3 ) ) );
    m_bins_q2.resize( abs( (int)getArg( 4 ) ) );
    m_bins_El.resize( abs( (int)getArg( 5 ) ) );
 
    int nextArg = EvtVubHybrid::nParameters + EvtVubHybrid::nVariables;
 
    m_nbins = m_bins_mX.size() * m_bins_q2.size() *
              m_bins_El.size();    // Binning of weight table
 
    int expectArgs = nextArg + m_bins_mX.size() + m_bins_q2.size() +
                     m_bins_El.size() + m_nbins;
 
    if ( getNArg() < expectArgs ) {
        EvtGenReport( EVTGEN_ERROR, "EvtVubHybrid" )
            << " finds " << getNArg() << " arguments, expected " << expectArgs
            << ".  Something is wrong with the tables of weights or thresholds."
            << "\nWill terminate execution!" << endl;
        ::abort();
    }
 
    // read bin boundaries from decay.dec
    for ( auto& b : m_bins_mX )
        b = getArg( nextArg++ );
    m_masscut = m_bins_mX[0];
 
    for ( auto& b : m_bins_q2 )
        b = getArg( nextArg++ );
    for ( auto& b : m_bins_El )
        b = getArg( nextArg++ );
 
    // read in weights (and rescale to range 0..1)
    readWeights( nextArg );
}
 
void EvtVubHybrid::initProbMax()
{
    noProbMax();
}
 
void EvtVubHybrid::decay( EvtParticle* p )
{
    int j;
    // B+ -> u-bar specflav l+ nu
 
    EvtParticle *xuhad( nullptr ), *lepton( nullptr ), *neutrino( nullptr );
    EvtVector4R p4;
    // R. Faccini 21/02/03
    // move the reweighting up , before also shooting the fermi distribution
    double x, z, p2;
    double sh = 0.0;
    double mB, ml, xlow, xhigh, qplus;
    double El = 0.0;
    double Eh = 0.0;
    double kplus;
    double q2, mX;
 
    const double lp2epsilon = -10;
    bool rew( true );
    while ( rew ) {
        p->initializePhaseSpace( getNDaug(), getDaugs() );
 
        xuhad = p->getDaug( 0 );
        lepton = p->getDaug( 1 );
        neutrino = p->getDaug( 2 );
 
        mB = p->mass();
        ml = lepton->mass();
 
        xlow = -m_mb;
        xhigh = mB - m_mb;
 
        // Fermi motion does not need to be computed inside the
        // tryit loop as m_b in Gamma0 does not need to be replaced by (m_b+kplus).
        // The difference however should be of the Order (lambda/m_b)^2 which is
        // beyond the considered orders in the paper anyway ...
 
        // for alpha_S = 0 and a mass cut on X_u not all values of kplus are
        // possible. The maximum value is mB/2-m_mb + sqrt(mB^2/4-m_masscut^2)
        kplus = 2 * xhigh;
 
        while ( kplus >= xhigh || kplus <= xlow ||
                ( m_alphas == 0 &&
                  kplus >= mB / 2 - m_mb +
                               sqrt( mB * mB / 4 - m_masscut * m_masscut ) ) ) {
            kplus = findPFermi();    //_pFermi->shoot();
            kplus = xlow + kplus * ( xhigh - xlow );
        }
        qplus = mB - m_mb - kplus;
        if ( ( mB - qplus ) / 2. <= ml )
            continue;
 
        int tryit = 1;
        while ( tryit ) {
            x = EvtRandom::Flat();
            z = EvtRandom::Flat( 0, 2 );
            p2 = EvtRandom::Flat();
            p2 = pow( 10, lp2epsilon * p2 );
 
            El = x * ( mB - qplus ) / 2;
            if ( El > ml && El < mB / 2 ) {
                Eh = z * ( mB - qplus ) / 2 + qplus;
                if ( Eh > 0 && Eh < mB ) {
                    sh = p2 * pow( mB - qplus, 2 ) +
                         2 * qplus * ( Eh - qplus ) + qplus * qplus;
                    if ( sh > m_masscut * m_masscut &&
                         mB * mB + sh - 2 * mB * Eh > ml * ml ) {
                        double xran = EvtRandom::Flat();
 
                        double y = m_dGamma->getdGdxdzdp( x, z, p2 ) / m_dGMax *
                                   p2;
 
                        if ( y > 1 )
                            EvtGenReport( EVTGEN_WARNING, "EvtVubHybrid" )
                                << "EvtVubHybrid decay probability > 1 found: "
                                << y << endl;
                        if ( y >= xran )
                            tryit = 0;
                    }
                }
            }
        }
 
        // compute all kinematic variables needed for reweighting (J. Dingfelder)
        mX = sqrt( sh );
        q2 = mB * mB + sh - 2 * mB * Eh;
 
        // Reweighting in bins of mX, q2, El (J. Dingfelder)
        if ( !m_weights.empty() ) {
            double xran1 = EvtRandom::Flat();
            double w = 1.0;
            if ( !m_noHybrid )
                w = getWeight( mX, q2, El );
            if ( w >= xran1 )
                rew = false;
        } else {
            rew = false;
        }
    }
 
    // o.k. we have the three kineamtic variables
    // now calculate a flat cos Theta_H [-1,1] distribution of the
    // hadron flight direction w.r.t the B flight direction
    // because the B is a scalar and should decay isotropic.
    // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
    // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
    // W flight direction.
 
    double ctH = EvtRandom::Flat( -1, 1 );
    double phH = EvtRandom::Flat( 0, 2 * M_PI );
    double phL = EvtRandom::Flat( 0, 2 * M_PI );
 
    // now compute the four vectors in the B Meson restframe
 
    double ptmp, sttmp;
    // calculate the hadron 4 vector in the B Meson restframe
 
    sttmp = sqrt( 1 - ctH * ctH );
    ptmp = sqrt( Eh * Eh - sh );
    double pHB[4] = { Eh, ptmp * sttmp * cos( phH ), ptmp * sttmp * sin( phH ),
                      ptmp * ctH };
    p4.set( pHB[0], pHB[1], pHB[2], pHB[3] );
    xuhad->init( getDaug( 0 ), p4 );
 
    if ( m_storeQplus ) {
        // cludge to store the hidden parameter q+ with the decay;
        // the lifetime of the Xu is abused for this purpose.
        // tau = 1 ps corresponds to ctau = 0.3 mm -> in order to
        // stay well below BaBars sensitivity we take q+/(10000 GeV) which
        // goes up to 0.0005 in the most extreme cases as ctau in mm.
        // To extract q+ back from the StdHepTrk its necessary to get
        // delta_ctau = Xu->anyDaughter->getVertexTime()-Xu->getVertexTime()
        // where these pseudo calls refere to the StdHep time stored at
        // the production vertex in the lab for each particle. The boost
        // has to be reversed and the result is:
        //
        // q+ = delta_ctau * 10000 GeV/mm * Mass_Xu/Energy_Xu
        //
        xuhad->setLifetime( qplus / 10000. );
    }
 
    // calculate the W 4 vector in the B Meson restrframe
 
    double apWB = ptmp;
    double pWB[4] = { mB - Eh, -pHB[1], -pHB[2], -pHB[3] };
 
    // first go in the W restframe and calculate the lepton and
    // the neutrino in the W frame
 
    double mW2 = mB * mB + sh - 2 * mB * Eh;
    double beta = ptmp / pWB[0];
    double gamma = pWB[0] / sqrt( mW2 );
 
    double pLW[4];
 
    ptmp = ( mW2 - ml * ml ) / 2 / sqrt( mW2 );
    pLW[0] = sqrt( ml * ml + ptmp * ptmp );
 
    double ctL = ( El - gamma * pLW[0] ) / beta / gamma / ptmp;
    if ( ctL < -1 )
        ctL = -1;
    if ( ctL > 1 )
        ctL = 1;
    sttmp = sqrt( 1 - ctL * ctL );
 
    // eX' = eZ x eW
    double xW[3] = { -pWB[2], pWB[1], 0 };
    // eZ' = eW
    double zW[3] = { pWB[1] / apWB, pWB[2] / apWB, pWB[3] / apWB };
 
    double lx = sqrt( xW[0] * xW[0] + xW[1] * xW[1] );
    for ( j = 0; j < 2; j++ )
        xW[j] /= lx;
 
    // eY' = eZ' x eX'
    double yW[3] = { -pWB[1] * pWB[3], -pWB[2] * pWB[3],
                     pWB[1] * pWB[1] + pWB[2] * pWB[2] };
    double ly = sqrt( yW[0] * yW[0] + yW[1] * yW[1] + yW[2] * yW[2] );
    for ( j = 0; j < 3; j++ )
        yW[j] /= ly;
 
    // p_lep = |p_lep| * (  sin(Theta) * cos(Phi) * eX'
    //                    + sin(Theta) * sin(Phi) * eY'
    //                    + cos(Theta) *            eZ')
    for ( j = 0; j < 3; j++ )
        pLW[j + 1] = sttmp * cos( phL ) * ptmp * xW[j] +
                     sttmp * sin( phL ) * ptmp * yW[j] + ctL * ptmp * zW[j];
 
    double apLW = ptmp;
    // calculate the neutrino 4 vector in the W restframe
    //double pNW[4] = {sqrt(mW2)-pLW[0],-pLW[1],-pLW[2],-pLW[3]};
 
    // boost them back in the B Meson restframe
 
    double appLB = beta * gamma * pLW[0] + gamma * ctL * apLW;
 
    ptmp = sqrt( El * El - ml * ml );
    double ctLL = appLB / ptmp;
 
    if ( ctLL > 1 )
        ctLL = 1;
    if ( ctLL < -1 )
        ctLL = -1;
 
    double pLB[4] = { El, 0, 0, 0 };
    double pNB[4] = { pWB[0] - El, 0, 0, 0 };
 
    for ( j = 1; j < 4; j++ ) {
        pLB[j] = pLW[j] + ( ctLL * ptmp - ctL * apLW ) / apWB * pWB[j];
        pNB[j] = pWB[j] - pLB[j];
    }
 
    p4.set( pLB[0], pLB[1], pLB[2], pLB[3] );
    lepton->init( getDaug( 1 ), p4 );
 
    p4.set( pNB[0], pNB[1], pNB[2], pNB[3] );
    neutrino->init( getDaug( 2 ), p4 );
 
    return;
}
 
double EvtVubHybrid::findPFermi()
{
    double ranNum = EvtRandom::Flat();
    double oOverBins = 1.0 / ( float( m_pf.size() ) );
    int nBinsBelow = 0;    // largest k such that I[k] is known to be <= rand
    int nBinsAbove = m_pf.size();    // largest k such that I[k] is known to be >  rand
    int middle;
 
    while ( nBinsAbove > nBinsBelow + 1 ) {
        middle = ( nBinsAbove + nBinsBelow + 1 ) >> 1;
        if ( ranNum >= m_pf[middle] ) {
            nBinsBelow = middle;
        } else {
            nBinsAbove = middle;
        }
    }
 
    double bSize = m_pf[nBinsAbove] - m_pf[nBinsBelow];
    // binMeasure is always aProbFunc[nBinsBelow],
 
    if ( bSize == 0 ) {
        // rand lies right in a bin of measure 0.  Simply return the center
        // of the range of that bin.  (Any value between k/N and (k+1)/N is
        // equally good, in this rare case.)
        return ( nBinsBelow + .5 ) * oOverBins;
    }
 
    double bFract = ( ranNum - m_pf[nBinsBelow] ) / bSize;
    return ( nBinsBelow + bFract ) * oOverBins;
}
 
double EvtVubHybrid::getWeight( double mX, double q2, double El )
{
    int ibin_mX = -1;
    int ibin_q2 = -1;
    int ibin_El = -1;
 
    for ( unsigned i = 0; i < m_bins_mX.size(); i++ ) {
        if ( mX >= m_bins_mX[i] )
            ibin_mX = i;
    }
    for ( unsigned i = 0; i < m_bins_q2.size(); i++ ) {
        if ( q2 >= m_bins_q2[i] )
            ibin_q2 = i;
    }
    for ( unsigned i = 0; i < m_bins_El.size(); i++ ) {
        if ( El >= m_bins_El[i] )
            ibin_El = i;
    }
    int ibin = ibin_mX + ibin_q2 * m_bins_mX.size() +
               ibin_El * m_bins_mX.size() * m_bins_q2.size();
 
    if ( ( ibin_mX < 0 ) || ( ibin_q2 < 0 ) || ( ibin_El < 0 ) ) {
        EvtGenReport( EVTGEN_ERROR, "EvtVubHybrid" )
            << "Cannot determine hybrid weight "
            << "for this event "
            << "-> assign weight = 0" << endl;
        return 0.0;
    }
 
    return m_weights[ibin];
}
 
void EvtVubHybrid::readWeights( int startArg )
{
    m_weights.resize( m_nbins );
 
    double maxw = 0.0;
    for ( auto& w : m_weights ) {
        w = getArg( startArg++ );
        if ( w > maxw )
            maxw = w;
    }
 
    if ( maxw == 0 ) {
        EvtGenReport( EVTGEN_ERROR, "EvtVubHybrid" )
            << "EvtVub generator expected at least one "
            << " weight > 0, but found none! "
            << "Will terminate execution!" << endl;
        ::abort();
    }
 
    // rescale weights (to be in range 0..1)
    for ( auto& w : m_weights )
        w /= maxw;
}