EvtVubNLO Class Reference

#include <EvtVubNLO.hh>

Inheritance diagram for EvtVubNLO:

Public Member Functions
	EvtVubNLO ()
virtual	~EvtVubNLO ()
void	getName (std::string &name)
EvtDecayBase *	clone ()
void	initProbMax ()
void	init ()
void	decay (EvtParticle *p)
Public Member Functions inherited from EvtDecayIncoherent
void	makeDecay (EvtParticle *p)
virtual	~EvtDecayIncoherent ()
void	setDaughterSpinDensity (int daughter)
int	isDaughterSpinDensitySet (int daughter)
Public Member Functions inherited from EvtDecayBase
virtual std::string	commandName ()
virtual void	command (std::string cmd)
double	getProbMax (double prob)
double	resetProbMax (double prob)
	EvtDecayBase ()
virtual	~EvtDecayBase ()
virtual bool	matchingDecay (const EvtDecayBase &other) const
EvtId	getParentId ()
double	getBranchingFraction ()
void	disableCheckQ ()
void	checkQ ()
int	getNDaug ()
EvtId *	getDaugs ()
EvtId	getDaug (int i)
int	getNArg ()
int	getPHOTOS ()
void	setPHOTOS ()
void	setVerbose ()
void	setSummary ()
double *	getArgs ()
std::string *	getArgsStr ()
double	getArg (int j)
std::string	getArgStr (int j)
std::string	getModelName ()
int	getDSum ()
int	summary ()
int	verbose ()
void	saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
void	printSummary ()
void	setProbMax (double prbmx)
void	noProbMax ()
void	checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1)
void	checkNDaug (int d1, int d2=-1)
void	checkSpinParent (EvtSpinType::spintype sp)
void	checkSpinDaughter (int d1, EvtSpinType::spintype sp)
virtual int	nRealDaughters ()

Additional Inherited Members
Static Public Member Functions inherited from EvtDecayBase
static void	findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10])
static void	findMass (EvtParticle *p)
static double	findMaxMass (EvtParticle *p)
Protected Member Functions inherited from EvtDecayBase
bool	daugsDecayedByParentModel ()
Protected Attributes inherited from EvtDecayBase
bool	_daugsDecayedByParentModel

Detailed Description

Definition at line 32 of file EvtVubNLO.hh.

Constructor & Destructor Documentation

EvtVubNLO()

EvtVubNLO::EvtVubNLO ( )

inline

Definition at line 35 of file EvtVubNLO.hh.

{}

Referenced by clone().

~EvtVubNLO()

EvtVubNLO::~EvtVubNLO ( )

virtual

Definition at line 40 of file EvtVubNLO.cc.

                      {
  delete[] _masses;
  delete[] _weights;
  cout << " max pdf : " << _gmax << endl;
  cout << " efficiency : " << (float)_ngood / (float)_ntot << endl;
}

Member Function Documentation

clone()

EvtDecayBase * EvtVubNLO::clone ( )

virtual

Implements EvtDecayBase.

Definition at line 53 of file EvtVubNLO.cc.

{ return new EvtVubNLO; }

decay()

void EvtVubNLO::decay ( EvtParticle * p )

virtual

Implements EvtDecayBase.

Definition at line 159 of file EvtVubNLO.cc.

                                      {
 
  int j;
  // B+ -> u-bar specflav l+ nu
 
  EvtParticle *xuhad, *lepton, *neutrino;
  EvtVector4R  p4;
 
  double pp, pm, pl, ml, El( 0.0 ), Eh( 0.0 ), sh( 0.0 );
 
  p->initializePhaseSpace( getNDaug(), getDaugs() );
 
  xuhad    = p->getDaug( 0 );
  lepton   = p->getDaug( 1 );
  neutrino = p->getDaug( 2 );
 
  _mB = p->mass();
  ml  = lepton->mass();
 
  bool tryit = true;
 
  while ( tryit )
  {
    // pm=(E_H+P_H)
    pm = EvtRandom::Flat( 0., 1 );
    pm = pow( pm, 1. / 3. ) * _mB;
    // pl=mB-2*El
    pl = EvtRandom::Flat( 0., 1 );
    pl = sqrt( pl ) * pm;
    // pp=(E_H-P_H)
    pp = EvtRandom::Flat( 0., pl );
 
    _ntot++;
 
    El = ( _mB - pl ) / 2.;
    Eh = ( pp + pm ) / 2;
    sh = pp * pm;
 
    double pdf( 0. );
    if ( pp < pl && El > ml && sh > _masses[0] * _masses[0] &&
         _mB * _mB + sh - 2 * _mB * Eh > ml * ml )
    {
      double xran = EvtRandom::Flat( 0, _dGMax );
      pdf         = tripleDiff( pp, pl, pm ); // triple differential distribution
      //      cout <<" P+,P-,Pl,Pdf= "<<pp <<" "<<pm<<" "<<pl<<" "<<pdf<<endl;
      if ( pdf > _dGMax )
      {
        report( ERROR, "EvtGen" )
            << "EvtVubNLO pdf above maximum: " << pdf << " P+,P-,Pl,Pdf= " << pp << " " << pm
            << " " << pl << " " << pdf << endl;
        //::abort();
      }
      if ( pdf >= xran ) tryit = false;
 
      if ( pdf > _gmax ) _gmax = pdf;
    }
    else
    {
      //      cout <<" EvtVubNLO incorrect kinematics  sh= "<<sh<<"EH "<<Eh<<endl;
    }
 
    // reweight the Mx distribution
    if ( !tryit && _nbins > 0 )
    {
      _ngood++;
      double xran1 = EvtRandom::Flat();
      double m     = sqrt( sh );
      j            = 0;
      while ( j < _nbins && m > _masses[j] ) j++;
      double w = _weights[j - 1];
      if ( w < xran1 ) tryit = true; // through away this candidate
    }
  }
 
  //  cout <<" max prob "<<gmax<<" " << pp<<" "<<y<<" "<<x<<endl;
 
  // 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 );
 
  // 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;
  //  if(mW2<0.1){
  //  cout <<" low Q2! "<<pp<<" "<<epp<<" "<<x<<" "<<y<<endl;
  //}
  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;
 
  // 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[8] = { 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;
}

getName()

void EvtVubNLO::getName ( std::string & name )

virtual

Implements EvtDecayBase.

Definition at line 51 of file EvtVubNLO.cc.

{ model_name = "VUB_NLO"; }

init()

void EvtVubNLO::init ( )

virtual

Reimplemented from EvtDecayBase.

Definition at line 55 of file EvtVubNLO.cc.

                     {
 
  // max pdf
  _gmax  = 0;
  _ntot  = 0;
  _ngood = 0;
  _lbar  = -1000;
  _mupi2 = -1000;
 
  // check that there are at least 6 arguments
  int npar = 8;
  if ( getNArg() < npar )
  {
 
    report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected "
                              << " at least npar arguments  but found: " << getNArg() << endl;
    report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
    ::abort();
  }
  // this is the shape function parameter
  _mb       = getArg( 0 );
  _b        = getArg( 1 );
  _lambdaSF = getArg( 2 );             // shape function lambda is different from lambda
  _mui      = 1.5;                     // GeV (scale)
  _kpar     = getArg( 3 );             // 0
  _idSF     = abs( (int)getArg( 4 ) ); // type of shape function 1: exponential (from Neubert)
  _nbins    = abs( (int)getArg( 5 ) );
  _masses   = new double[_nbins];
  _weights  = new double[_nbins];
 
  // Shape function normalization
  _mB = 5.28; // temporary B meson mass for normalization
 
  std::vector<double> sCoeffs( 11 );
  sCoeffs[3]  = _b;
  sCoeffs[4]  = _mb;
  sCoeffs[5]  = _mB;
  sCoeffs[6]  = _idSF;
  sCoeffs[7]  = lambda_SF();
  sCoeffs[8]  = mu_h();
  sCoeffs[9]  = mu_i();
  sCoeffs[10] = 1.;
  _SFNorm     = SFNorm( sCoeffs ); // SF normalization;
 
  cout << " pdf 0.66, 1.32 , 4.32 " << tripleDiff( 0.66, 1.32, 4.32 ) << endl;
  cout << " pdf 0.23,0.37,3.76 " << tripleDiff( 0.23, 0.37, 3.76 ) << endl;
  cout << " pdf 0.97,4.32,4.42 " << tripleDiff( 0.97, 4.32, 4.42 ) << endl;
  cout << " pdf 0.52,1.02,2.01 " << tripleDiff( 0.52, 1.02, 2.01 ) << endl;
  cout << " pdf 1.35,1.39,2.73 " << tripleDiff( 1.35, 1.39, 2.73 ) << endl;
 
  if ( getNArg() - npar + 2 != 2 * _nbins )
  {
    report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected " << _nbins
                              << " masses and weights but found: " << ( getNArg() - npar ) / 2
                              << endl;
    report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
    ::abort();
  }
  int    i, j = npar - 2;
  double maxw = 0.;
  for ( i = 0; i < _nbins; i++ )
  {
    _masses[i] = getArg( j++ );
    if ( i > 0 && _masses[i] <= _masses[i - 1] )
    {
      report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected "
                                << " mass bins in ascending order!"
                                << "Will terminate execution!" << endl;
      ::abort();
    }
    _weights[i] = getArg( j++ );
    if ( _weights[i] < 0 )
    {
      report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected "
                                << " weights >= 0, but found: " << _weights[i] << endl;
      report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
      ::abort();
    }
    if ( _weights[i] > maxw ) maxw = _weights[i];
  }
  if ( maxw == 0 )
  {
    report( ERROR, "EvtGen" ) << "EvtVubNLO generator expected at least one "
                              << " weight > 0, but found none! "
                              << "Will terminate execution!" << endl;
    ::abort();
  }
  for ( i = 0; i < _nbins; i++ ) _weights[i] /= maxw;
 
  // the maximum dGamma*p2 value depends on alpha_s only:
 
  //  _dGMax = 0.05;
  _dGMax = 150.;
 
  // for the Fermi Motion we need a B-Meso\n 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
 
  // check that there are 3 daughters
  checkNDaug( 3 );
}

initProbMax()

void EvtVubNLO::initProbMax ( )

virtual

Reimplemented from EvtDecayBase.

Definition at line 157 of file EvtVubNLO.cc.

{ noProbMax(); }

---

The documentation for this class was generated from the following files:

• EvtVubNLO.hh
• EvtVubNLO.cc