EvtVub Class Reference

#include <EvtVub.hh>

Inheritance diagram for EvtVub:

List of all members.

Public Member Functions
void	checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1)
void	checkNDaug (int d1, int d2=-1)
void	checkQ ()
void	checkSpinDaughter (int d1, EvtSpinType::spintype sp)
void	checkSpinParent (EvtSpinType::spintype sp)
EvtDecayBase *	clone ()
virtual void	command (std::string cmd)
virtual std::string	commandName ()
void	decay (EvtParticle *p)
void	disableCheckQ ()
	EvtVub ()
double	getArg (int j)
double *	getArgs ()
std::string *	getArgsStr ()
std::string	getArgStr (int j)
double	getBranchingFraction ()
EvtId	getDaug (int i)
EvtId *	getDaugs ()
int	getDSum ()
std::string	getModelName ()
void	getName (std::string &name)
int	getNArg ()
int	getNDaug ()
EvtId	getParentId ()
int	getPHOTOS ()
double	getProbMax (double prob)
void	init ()
void	initProbMax ()
int	isDaughterSpinDensitySet (int daughter)
void	makeDecay (EvtParticle *p)
virtual bool	matchingDecay (const EvtDecayBase &other) const
void	noProbMax ()
virtual int	nRealDaughters ()
void	printSummary ()
double	resetProbMax (double prob)
void	saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
void	setDaughterSpinDensity (int daughter)
void	setPHOTOS ()
void	setProbMax (double prbmx)
void	setSummary ()
void	setVerbose ()
int	summary ()
int	verbose ()
virtual	~EvtVub ()
Static Public Member Functions
void	findMass (EvtParticle *p)
void	findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10])
double	findMaxMass (EvtParticle *p)
Protected Member Functions
bool	daugsDecayedByParentModel ()
Protected Attributes
bool	_daugsDecayedByParentModel
Private Member Functions
double	findPFermi ()
Private Attributes
double	_a
double	_alphas
EvtVubdGamma *	_dGamma
double	_dGMax
double *	_masses
double	_mb
int	_nbins
std::vector< double >	_pf
int	_storeQplus
double *	_weights

---

Constructor & Destructor Documentation

EvtVub::EvtVub (   )  [inline]

{}

EvtVub::~EvtVub (   )  [virtual]

{ delete _dGamma; delete [] _masses; delete [] _weights; }

---

Member Function Documentation

void EvtDecayBase::checkNArg (  int  a1, int  a2 = -1, int  a3 = -1, int  a4 = -1 )  [inherited]

{ if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) { report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl; report(ERROR,"EvtGen") << a1<<endl;; if ( a2>-1) { report(ERROR,"EvtGen") << " or " << a2<<endl; } if ( a3>-1) { report(ERROR,"EvtGen") << " or " << a3<<endl; } if ( a4>-1) { report(ERROR,"EvtGen") << " or " << a4<<endl; } report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

void EvtDecayBase::checkNDaug (  int  d1, int  d2 = -1 )  [inherited]

{ if ( _ndaug != d1 && _ndaug != d2 ) { report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "; report(ERROR,"EvtGen") << d1; if ( d2>-1) { report(ERROR,"EvtGen") << " or " << d2; } report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

void EvtDecayBase::checkQ (   )  [inherited]

{ int i; int q=0; int qpar; //If there are no daughters (jetset etc) then we do not //want to do this test. Why? Because sometimes the parent //will have a nonzero charge. if ( _ndaug != 0) { for(i=0; i<_ndaug; i++ ) { q += EvtPDL::chg3(_daug[i]); } qpar = EvtPDL::chg3(_parent); if ( q != qpar ) { report(ERROR,"EvtGen") <<_modelname.c_str()<< " generator expected " << " charge to be conserved, found:"<<endl; report(ERROR,"EvtGen") << "Parent charge of "<<(qpar/3)<<endl; report(ERROR,"EvtGen") << "Sum of daughter charge of "<<(q/3)<<endl; report(ERROR,"EvtGen") << "The parent is "<< EvtPDL::name(_parent).c_str()<<endl; for(i=0; i<_ndaug; i++ ) { report(ERROR,"EvtGen") << "Daughter "<< EvtPDL::name(_daug[i]).c_str()<<endl; } report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } } }

void EvtDecayBase::checkSpinDaughter (  int  d1, EvtSpinType::spintype  sp )  [inherited]

{ EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1)); if ( parenttype != sp ) { report(ERROR,"EvtGen") << _modelname.c_str() << " did not get the correct daughter spin d=" << d1 << endl; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

void EvtDecayBase::checkSpinParent (  EvtSpinType::spintype  sp  )  [inherited]

{ EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId()); if ( parenttype != sp ) { report(ERROR,"EvtGen") << _modelname.c_str() << " did not get the correct parent spin\n"; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

EvtDecayBase * EvtVub::clone (   )  [virtual]

Implements EvtDecayBase. { return new EvtVub; }

void EvtDecayBase::command (  std::string  cmd  )  [virtual, inherited]

Reimplemented in EvtJetSet, EvtLunda, EvtLundCharm, EvtPythia, and EvtTauola. { report(ERROR,"EvtGen") << "Should never call EvtDecayBase::command"<<endl; ::abort(); }

std::string EvtDecayBase::commandName (   )  [virtual, inherited]

Reimplemented in EvtJetSet, EvtLunda, EvtLundCharm, EvtPythia, and EvtTauola. { return std::string(""); }

bool EvtDecayBase::daugsDecayedByParentModel (   )  [inline, protected, inherited]

{return _daugsDecayedByParentModel;}

void EvtVub::decay (  EvtParticle *  p  )  [virtual]

Implements EvtDecayBase. { 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; 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-_masses[0]^2) kplus = 2*xhigh; while( kplus >= xhigh || kplus <= xlow || (_alphas == 0 && kplus >= mB/2-_mb + sqrt(mB*mB/4-_masses[0]*_masses[0]))) { 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 > _masses[0]*_masses[0] && 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,"EvtGen")<<"EvtVub decay probability > 1 found: " << y << endl; if ( y >= xran ) tryit = 0; } } } } // reweight the Mx distribution if(_nbins>0){ 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 ) 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 ; }

void EvtDecayBase::disableCheckQ (   )  [inline, inherited]

{_chkCharge=0;};

void EvtDecayBase::findMass (  EvtParticle *  p  )  [static, inherited]

{ //Need to also check that this mass does not screw //up the parent //This code assumes that for the ith daughter, 0..i-1 //already have a mass double maxOkMass=findMaxMass(p); int count=0; double mass; bool massOk=false; int i; while (!massOk) { count++; if ( count > 10000 ) { report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl; report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n"; if ( p->getParent() ) { if ( p->getParent()->getParent() ) { p->getParent()->getParent()->printTree(); report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl; report(INFO,"EvtGen") << p->getParent()->mass() <<endl; } else{ p->getParent()->printTree(); report(INFO,"EvtGen") << p->getParent()->mass() <<endl; } } else p->printTree(); report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl; if ( p->getNDaug() ) { for (i=0; i<p->getNDaug(); i++) { report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" "; } report(INFO,"EvtGen") << endl; } if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) { report(INFO,"EvtGen") << "taking a default value\n"; p->setMass(maxOkMass); return; } assert(0); } mass = EvtPDL::getMass(p->getId()); //Just need to check that this mass is > than //the mass of all daughters double massSum=0.; if ( p->getNDaug() ) { for (i=0; i<p->getNDaug(); i++) { massSum+= p->getDaug(i)->mass(); } } //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters if (p->getNDaug()<2) massOk=true; if ( p->getParent() ) { if ( p->getParent()->getNDaug()==1 ) massOk=true; } if ( !massOk ) { if (massSum < mass) massOk=true; if ( mass> maxOkMass) massOk=false; } } p->setMass(mass); }

void EvtDecayBase::findMasses (  EvtParticle *  p, int  ndaugs, EvtId  daugs[10], double  masses[10] )  [static, inherited]

{ int i; double mass_sum; int count=0; if (!( p->firstornot() )) { for (i = 0; i < ndaugs; i++ ) { masses[i] = p->getDaug(i)->mass(); } //for } //if else { p->setFirstOrNot(); // if only one daughter do it if (ndaugs==1) { masses[0]=p->mass(); return; } //until we get a combo whose masses are less than _parent mass. do { mass_sum = 0.0; for (i = 0; i < ndaugs; i++ ) { masses[i] = EvtPDL::getMass(daugs[i]); mass_sum = mass_sum + masses[i]; } count++; if(count==10000) { report(ERROR,"EvtGen") <<"Decaying particle:"<< EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl; report(ERROR,"EvtGen") <<"To the following daugthers"<<endl; for (i = 0; i < ndaugs; i++ ) { report(ERROR,"EvtGen") << EvtPDL::name(daugs[i]).c_str() << endl; } report(ERROR,"EvtGen") << "Has been rejected "<<count << " times, will now take minimal masses " << " of daugthers"<<endl; mass_sum=0.; for (i = 0; i < ndaugs; i++ ) { masses[i] = EvtPDL::getMinMass(daugs[i]); mass_sum = mass_sum + masses[i]; } if (mass_sum > p->mass()){ report(ERROR,"EvtGen") << "Parent mass="<<p->mass() << "to light for daugthers."<<endl << "Will throw the event away."<<endl; //dont terminate - start over on the event. EvtStatus::setRejectFlag(); mass_sum=0.; // ::abort(); } } } while ( mass_sum > p->mass()); } //else return; }

double EvtDecayBase::findMaxMass (  EvtParticle *  p  )  [static, inherited]

{ double maxOkMass=EvtPDL::getMaxMass(p->getId()); //protect against vphotons if ( maxOkMass < 0.0000000001 ) return 10000000.; //and against already determined masses if ( p->hasValidP4() ) maxOkMass=p->mass(); EvtParticle *par=p->getParent(); if ( par ) { double maxParMass=findMaxMass(par); int i; double minDaugMass=0.; for(i=0;i<par->getNDaug();i++){ EvtParticle *dau=par->getDaug(i); if ( dau!=p) { // it might already have a mass if ( dau->isInitialized() || dau->hasValidP4() ) minDaugMass+=dau->mass(); else //give it a bit of phase space minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId()); } } if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass; } return maxOkMass; }

double EvtVub::findPFermi (   )  [private]

{ double ranNum=EvtRandom::Flat(); double oOverBins= 1.0/(float(_pf.size())); int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand int nBinsAbove = _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 >= _pf[middle]) { nBinsBelow = middle; } else { nBinsAbove = middle; } } double bSize = _pf[nBinsAbove] - _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; } HepDouble bFract = (ranNum - _pf[nBinsBelow]) / bSize; return (nBinsBelow + bFract) * oOverBins; }

double EvtDecayBase::getArg (  int  j  )  [inherited]

{ // Verify string const char* str = _args[j].c_str(); int i = 0; while(str[i]!=0){ if (isalpha(str[i]) && str[i]!='e') { report(INFO,"EvtGen") << "String " << str << " is not a number" << endl; assert(0); } i++; } char** tc=0; return strtod(_args[j].c_str(),tc); }

double * EvtDecayBase::getArgs (   )  [inherited]

{ if ( _argsD ) return _argsD; //The user has asked for a list of doubles - the arguments //better all be doubles... if ( _narg==0 ) return _argsD; _argsD = new double[_narg]; int i; char * tc; for(i=0;i<_narg;i++) { _argsD[i] = strtod(_args[i].c_str(),&tc); } return _argsD; }

std::string* EvtDecayBase::getArgsStr (   )  [inline, inherited]

{return _args;}

std::string EvtDecayBase::getArgStr (  int  j  )  [inline, inherited]

{return _args[j];}

double EvtDecayBase::getBranchingFraction (   )  [inline, inherited]

{return _brfr;}

EvtId EvtDecayBase::getDaug (  int  i  )  [inline, inherited]

{return _daug[i];}

EvtId* EvtDecayBase::getDaugs (   )  [inline, inherited]

{return _daug;}

int EvtDecayBase::getDSum (   )  [inline, inherited]

{return _dsum; }

std::string EvtDecayBase::getModelName (   )  [inline, inherited]

{return _modelname; }

void EvtVub::getName (  std::string &  name  )  [virtual]

Implements EvtDecayBase. { model_name="VUB"; }

int EvtDecayBase::getNArg (   )  [inline, inherited]

{return _narg;}

int EvtDecayBase::getNDaug (   )  [inline, inherited]

{return _ndaug;}

EvtId EvtDecayBase::getParentId (   )  [inline, inherited]

{return _parent;}

int EvtDecayBase::getPHOTOS (   )  [inline, inherited]

{return _photos;}

double EvtDecayBase::getProbMax (  double  prob  )  [inherited]

{ int i; //diagnostics sum_prob+=prob; if (prob>max_prob) max_prob=prob; if ( defaultprobmax && ntimes_prob<=500 ) { //We are building up probmax with this iteration ntimes_prob += 1; if ( prob > probmax ) { probmax = prob;} if (ntimes_prob==500) { probmax*=1.2; } return 1000000.0*prob; } if ( prob> probmax*1.0001) { report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")"; report(INFO,"") << "("<<_modelname.c_str()<<") "; report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; for(i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; } report(INFO,"") << endl; if (defaultprobmax) probmax = prob; } ntimes_prob += 1; return probmax; } //getProbMax

void EvtVub::init (   )  [virtual]

Reimplemented from EvtDecayBase. { // check that there are at least 6 arguments if (getNArg()<6) { report(ERROR,"EvtGen") << "EvtVub generator expected " << " at least 6 arguments (mb,a,alpha_s,Nbins,m1,w1,...) but found: " <<getNArg()<<endl; report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } _mb = getArg(0); _a = getArg(1); _alphas = getArg(2); _nbins = abs((int)getArg(3)); _storeQplus = (getArg(3)<0?1:0); _masses = new double[_nbins]; _weights = new double[_nbins]; if (getNArg()-4 != 2*_nbins) { report(ERROR,"EvtGen") << "EvtVub generator expected " << _nbins << " masses and weights but found: " <<(getNArg()-4)/2 <<endl; report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } int i,j = 4; double maxw = 0.; for (i=0;i<_nbins;i++) { _masses[i] = getArg(j++); if (i>0 && _masses[i] <= _masses[i-1]) { report(ERROR,"EvtGen") << "EvtVub generator expected " << " mass bins in ascending order!" << "Will terminate execution!"<<endl; ::abort(); } _weights[i] = getArg(j++); if (_weights[i] < 0) { report(ERROR,"EvtGen") << "EvtVub 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") << "EvtVub 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: 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 EvtId BP=EvtPDL::getId("B+"); EvtId B0=EvtPDL::getId("B0"); double mB0 = EvtPDL::getMaxMass(B0); double mBP = EvtPDL::getMaxMass(BP); 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(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 (i=0; i<_pf.size(); i++) _pf[i]/=_pf[_pf.size()-1]; // static EvtHepRandomEngine myEngine; // _pFermi = new RandGeneral(myEngine,pf,aSize,0); _dGamma = new EvtVubdGamma(_alphas); // check that there are 3 daughters checkNDaug(3); }

void EvtVub::initProbMax (   )  [virtual]

Reimplemented from EvtDecayBase. { noProbMax(); }

int EvtDecayIncoherent::isDaughterSpinDensitySet (  int  daughter  )  [inline, inherited]

{return spinDensitySet[daughter];}

void EvtDecayIncoherent::makeDecay (  EvtParticle *  p  )  [virtual, inherited]

Implements EvtDecayBase. { int i; //initialize this the hard way.. //Lange June 26, 2000 for (i=0; i<MAX_DAUG; i++ ) { spinDensitySet[i]=0;} _daugsDecayedByParentModel=false; decay(p); p->setDecayProb(1.); EvtSpinDensity rho; rho.SetDiag(p->getSpinStates()); p->setSpinDensityBackward(rho); if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) { EvtRadCorr::doRadCorr(p); } //Now decay the daughters. if ( !daugsDecayedByParentModel()) { for(i=0;i<p->getNDaug();i++){ //Need to set the spin density of the daughters to be //diagonal. rho.SetDiag(p->getDaug(i)->getSpinStates()); //if (p->getDaug(i)->getNDaug()==0){ //only do this if the user has not already set the //spin density matrix herself. //Lange June 26, 2000 if ( isDaughterSpinDensitySet(i)==0 ) { p->getDaug(i)->setSpinDensityForward(rho); } else{ //report(INFO,"EvtGen") << "spinDensitymatrix already set!!!\n"; EvtSpinDensity temp=p->getDaug(i)->getSpinDensityForward(); // report(INFO,"EvtGen") <<temp<<endl; } //Now decay the daughter. Really! p->getDaug(i)->decay(); //} } } }

bool EvtDecayBase::matchingDecay (  const EvtDecayBase &  other  )  const [virtual, inherited]

{ if ( _ndaug != other._ndaug) return false; if ( _parent != other._parent) return false; std::vector<int> useDs; for ( unsigned int i=0; i<_ndaug; i++) useDs.push_back(0); for ( unsigned int i=0; i<_ndaug; i++) { bool foundIt=false; for ( unsigned int j=0; j<_ndaug; j++) { if ( useDs[j] == 1 ) continue; if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) { foundIt=true; useDs[j]=1; break; } } if ( foundIt==false) return false; } for ( unsigned int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false; return true; }

void EvtDecayBase::noProbMax (   )  [inherited]

{ defaultprobmax=0; }

virtual int EvtDecayBase::nRealDaughters (   )  [inline, virtual, inherited]

Reimplemented in EvtBtoKD3P, and EvtVSSBMixCPT. { return _ndaug;}

void EvtDecayBase::printSummary (   )  [inherited]

{ int i; if (ntimes_prob>0) { report(INFO,"EvtGen") << "Calls="<<ntimes_prob<<" eff:"<< sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax; report(INFO,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : "; } report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; for(i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; } report(INFO,"") << " ("<<_modelname.c_str()<<"):"<< endl; }

double EvtDecayBase::resetProbMax (  double  prob  )  [inherited]

{ report(INFO,"EvtGen") << "Reseting prob max\n"; report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")"; report(INFO,"") << "("<<_modelname.c_str()<<")"; report(INFO,"") << EvtPDL::getStdHep(_parent)<<"->"; for( int i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::getStdHep(_daug[i]) << " "; } report(INFO,"") << endl; probmax = 0.0; defaultprobmax = 0; ntimes_prob = 0; return prob; }

void EvtDecayBase::saveDecayInfo (  EvtId  ipar, int  ndaug, EvtId *  daug, int  narg, std::vector< std::string > &  args, std::string  name, double  brfr )  [inherited]

{ int i; _brfr=brfr; _ndaug=ndaug; _narg=narg; _parent=ipar; _dsum=0; if (_ndaug>0) { _daug=new EvtId [_ndaug]; for(i=0;i<_ndaug;i++){ _daug[i]=daug[i]; _dsum+=daug[i].getAlias(); } } else{ _daug=0; } if (_narg>0) { _args=new std::string[_narg+1]; for(i=0;i<_narg;i++){ _args[i]=args[i]; } } else{ _args = 0; } _modelname=name; this->init(); this->initProbMax(); if (_chkCharge){ this->checkQ(); } if (defaultprobmax){ report(INFO,"EvtGen") << "No default probmax for "; report(INFO,"") << "("<<_modelname.c_str()<<") "; report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; for(i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; } report(INFO,"") << endl; report(INFO,"") << "This is fine for development, but must be provided for production."<<endl; report(INFO,"EvtGen") << "Never fear though - the decay will use the \n"; report(INFO,"EvtGen") << "500 iterations to build up a good probmax \n"; report(INFO,"EvtGen") << "before accepting a decay. "<<endl; } }

void EvtDecayIncoherent::setDaughterSpinDensity (  int  daughter  )  [inline, inherited]

{ spinDensitySet[daughter]=1; return;}

void EvtDecayBase::setPHOTOS (   )  [inline, inherited]

{_photos=1;}

void EvtDecayBase::setProbMax (  double  prbmx  )  [inherited]

{ defaultprobmax=0; probmax=prbmx; }

void EvtDecayBase::setSummary (   )  [inline, inherited]

{_summary=1;}

void EvtDecayBase::setVerbose (   )  [inline, inherited]

{_verbose=1;}

int EvtDecayBase::summary (   )  [inline, inherited]

{return _summary; }

int EvtDecayBase::verbose (   )  [inline, inherited]

{return _verbose; }

---

Member Data Documentation

double EvtVub::_a [private]

double EvtVub::_alphas [private]

bool EvtDecayBase::_daugsDecayedByParentModel [protected, inherited]

EvtVubdGamma* EvtVub::_dGamma [private]

double EvtVub::_dGMax [private]

double* EvtVub::_masses [private]

double EvtVub::_mb [private]

int EvtVub::_nbins [private]

std::vector<double> EvtVub::_pf [private]

int EvtVub::_storeQplus [private]

double* EvtVub::_weights [private]

---

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

• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/src/EvtGen/EvtGenModels/EvtVub.hh
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/src/EvtGen/EvtGenModels/EvtVub.cc

---