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 void	getName (std::string &name)=0
virtual void	decay (EvtParticle *p)=0
virtual void	makeDecay (EvtParticle *p)=0
virtual EvtDecayBase *	clone ()=0
virtual void	init ()
virtual void	initProbMax ()
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 36 of file EvtVubNLO.hh.

{}

Referenced by clone().

~EvtVubNLO()

EvtVubNLO::~EvtVubNLO ( )

virtual

Definition at line 39 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 57 of file EvtVubNLO.cc.

                              {
 
  return new EvtVubNLO;
 
}

decay()

void EvtVubNLO::decay ( EvtParticle *  p )

virtual

Implements EvtDecayBase.

Definition at line 173 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 64 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 167 of file EvtVubNLO.cc.

                           {
  
  noProbMax();
  
}

---

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

• EvtVubNLO.hh
• EvtVubNLO.cc