EvtVubHybrid Class Reference

#include <EvtVubHybrid.hh>

Inheritance diagram for EvtVubHybrid:

Public Member Functions
	EvtVubHybrid ()
virtual	~EvtVubHybrid ()
void	getName (std::string &name)
EvtDecayBase *	clone ()
void	initProbMax ()
void	init ()
void	decay (EvtParticle *p)
void	readWeights (int startArg=0)
double	getWeight (double mX, double q2, double El)
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 44 of file EvtVubHybrid.hh.

Constructor & Destructor Documentation

EvtVubHybrid()

EvtVubHybrid::EvtVubHybrid

(

)

Definition at line 51 of file EvtVubHybrid.cc.

  : _noHybrid(false), _storeQplus(true),
    _mb(4.62), _a(2.27), _alphas(0.22), _dGMax(3.),
    _nbins_mX(0), _nbins_q2(0), _nbins_El(0), _nbins(0),
    _masscut(0.28), _bins_mX(0), _bins_q2(0), _bins_El(0),
    _weights(0), _dGamma(0)
{}

Referenced by clone().

~EvtVubHybrid()

EvtVubHybrid::~EvtVubHybrid ( )

virtual

Definition at line 59 of file EvtVubHybrid.cc.

                            {
  delete _dGamma;
  delete [] _bins_mX;
  delete [] _bins_q2;
  delete [] _bins_El;
  delete [] _weights;
 
}

Member Function Documentation

clone()

EvtDecayBase * EvtVubHybrid::clone ( )

virtual

Implements EvtDecayBase.

Definition at line 74 of file EvtVubHybrid.cc.

                                 {
 
  return new EvtVubHybrid;
 
}

decay()

void EvtVubHybrid::decay ( EvtParticle *  p )

virtual

Implements EvtDecayBase.

Definition at line 187 of file EvtVubHybrid.cc.

                                        {
 
  int j;
  // B+ -> u-bar specflav l+ nu
  
  EvtParticle *xuhad, *lepton, *neutrino;
  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 = -_mb;
    xhigh = mB-_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-_mb + sqrt(mB^2/4-_masscut^2)
    kplus = 2*xhigh;
    
    while( kplus >= xhigh || kplus <= xlow 
       || (_alphas == 0 && kplus >= mB/2-_mb 
           + sqrt(mB*mB/4-_masscut*_masscut))) {
      kplus = findPFermi(); //_pFermi->shoot();
      kplus = xlow + kplus*(xhigh-xlow);
    }
    qplus = mB-_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 > _masscut*_masscut
           && mB*mB + sh - 2*mB*Eh > ml*ml) {
        
        double xran = EvtRandom::Flat();
        
        double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2;
        
        if ( y > 1 ) report(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 (_nbins>0) {
      double xran1 = EvtRandom::Flat();
      double w = 1.0;
      if (!_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 (_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 ;
}

getName()

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

virtual

Implements EvtDecayBase.

Definition at line 68 of file EvtVubHybrid.cc.

                                               {
 
  model_name="VUBHYBRID";     
 
}

getWeight()

double EvtVubHybrid::getWeight	(	double	mX,
		double	q2,
		double	El
	)

Definition at line 435 of file EvtVubHybrid.cc.

                                                              {
 
  int ibin_mX = -1;
  int ibin_q2 = -1;
  int ibin_El = -1;
 
  for (int i = 0; i < _nbins_mX; i++) {
    if (mX >= _bins_mX[i]) ibin_mX = i;
  }
  for (int i = 0; i < _nbins_q2; i++) {
    if (q2 >= _bins_q2[i]) ibin_q2 = i;
  }
  for (int i = 0; i < _nbins_El; i++) {
    if (El >= _bins_El[i]) ibin_El = i;
  }
  int ibin = ibin_mX + ibin_q2*_nbins_mX + ibin_El*_nbins_mX*_nbins_q2;
 
  if ( (ibin_mX < 0) || (ibin_q2 < 0) || (ibin_El < 0) ) {
    report(ERROR,"EvtVubHybrid") << "Cannot determine hybrid weight "
                                 << "for this event " 
                 << "-> assign weight = 0" << endl;
    return 0.0;
  }
 
  return _weights[ibin];
}

Referenced by decay().

init()

void EvtVubHybrid::init ( )

virtual

Reimplemented from EvtDecayBase.

Definition at line 80 of file EvtVubHybrid.cc.

                       {
 
  // check that there are at least 3 arguments
  if (getNArg() < EvtVubHybrid::nParameters) {
    report(ERROR,"EvtVubHybrid") << "EvtVub generator expected "
                 << "at least " << EvtVubHybrid::nParameters
                 << " arguments but found: " << getNArg()
                 << "\nWill terminate execution!"<<endl;
    ::abort(); 
  } else if (getNArg() == EvtVubHybrid::nParameters) {
    report(WARNING,"EvtVubHybrid") << "EvtVub: generate B -> Xu l nu events " 
                   << "without using the hybrid reweighting." 
                   << endl;
    _noHybrid = true;
  } else if (getNArg() < EvtVubHybrid::nParameters+EvtVubHybrid::nVariables) {
     report(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
  _mb     = getArg(0);
  _a      = getArg(1);
  _alphas = getArg(2);
 
  // the maximum dGamma*p2 value depends on alpha_s only:
  const double dGMax0 = 3.;
  _dGMax = 0.21344+8.905*_alphas;
  if ( _dGMax < dGMax0 ) _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 double mB0 = EvtPDL::getMaxMass(EvtPDL::getId("B0"));
  static double mBP = EvtPDL::getMaxMass(EvtPDL::getId("B+"));
  static double mB = (mB0<mBP?mB0:mBP);
  
  const double xlow = -_mb;
  const double xhigh = mB-_mb;
  const int aSize = 10000;
 
  EvtPFermi pFermi(_a,mB,_mb);
  // pf is the cumulative distribution normalized to 1.
  _pf.resize(aSize);
  for(int i=0;i<aSize;i++){
    double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
    if ( i== 0 )
      _pf[i] = pFermi.getFPFermi(kplus);
    else
      _pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus);
  }
  for (int i=0; i<_pf.size(); i++) _pf[i]/=_pf[_pf.size()-1];
 
  _dGamma = new EvtVubdGamma(_alphas);
  
  if (_noHybrid) return;    // Without hybrid weighting, nothing else to do
 
  _nbins_mX = abs((int)getArg(3));
  _nbins_q2 = abs((int)getArg(4));
  _nbins_El = abs((int)getArg(5));
 
  int nextArg = EvtVubHybrid::nParameters + EvtVubHybrid::nVariables;
 
  _nbins = _nbins_mX*_nbins_q2*_nbins_El;   // Binning of weight table
 
  int expectArgs = nextArg + _nbins_mX +_nbins_q2 + _nbins_El + _nbins;
 
  if (getNArg() < expectArgs) {
    report(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
  int i;
 
  _bins_mX  = new double[_nbins_mX];
  for (i = 0; i < _nbins_mX; i++,nextArg++) {
    _bins_mX[i] = getArg(nextArg);
  }
  _masscut = _bins_mX[0];
 
  _bins_q2  = new double[_nbins_q2];
  for (i = 0; i < _nbins_q2; i++,nextArg++) {
    _bins_q2[i] = getArg(nextArg);    
  }
 
  _bins_El  = new double[_nbins_El];
  for (i = 0; i < _nbins_El; i++,nextArg++) {
    _bins_El[i] = getArg(nextArg);    
  }
    
  // read in weights (and rescale to range 0..1)
  readWeights(nextArg); 
}

initProbMax()

void EvtVubHybrid::initProbMax ( )

virtual

Reimplemented from EvtDecayBase.

Definition at line 183 of file EvtVubHybrid.cc.

                              {
  noProbMax();
}

readWeights()

void EvtVubHybrid::readWeights

(

int

startArg = 0

)

Definition at line 463 of file EvtVubHybrid.cc.

                                           {
  _weights  = new double[_nbins];
 
  double maxw = 0.0;
  for (int i = 0; i < _nbins; i++, startArg++) {
    _weights[i] = getArg(startArg);
    if (_weights[i] > maxw) maxw = _weights[i];
  }
 
  if (maxw == 0) {
    report(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 (int i = 0; i < _nbins; i++) {
    _weights[i] /= maxw;
  }
}

Referenced by init().

---

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

• EvtVubHybrid.hh
• EvtVubHybrid.cc