EvtVubHybrid Class Reference

#include <EvtVubHybrid.hh>

Inheritance diagram for EvtVubHybrid:

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 ()
	EvtVubHybrid ()
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)
double	getWeight (double mX, double q2, double El)
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 ()
void	readWeights (int startArg=0)
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	~EvtVubHybrid ()
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 Types
enum	{ nParameters = 3, nVariables = 3 }
Private Member Functions
double	findPFermi ()
Private Attributes
double	_a
double	_alphas
double *	_bins_El
double *	_bins_mX
double *	_bins_q2
EvtVubdGamma *	_dGamma
double	_dGMax
double	_masscut
double	_mb
int	_nbins
int	_nbins_El
int	_nbins_mX
int	_nbins_q2
bool	_noHybrid
std::vector< double >	_pf
bool	_storeQplus
double *	_weights

---

Member Enumeration Documentation

anonymous enum [private]

Enumeration values: nParameters  nVariables  { nParameters = 3, nVariables = 3 };

---

Constructor & Destructor Documentation

EvtVubHybrid::EvtVubHybrid (   )

: _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) {}

EvtVubHybrid::~EvtVubHybrid (   )  [virtual]

{ delete _dGamma; delete [] _bins_mX; delete [] _bins_q2; delete [] _bins_El; 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 * EvtVubHybrid::clone (   )  [virtual]

Implements EvtDecayBase. { return new EvtVubHybrid; }

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 EvtVubHybrid::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; 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 ; }

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 EvtVubHybrid::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 EvtVubHybrid::getName (  std::string &  name  )  [virtual]

Implements EvtDecayBase. { model_name="VUBHYBRID"; }

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

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

{ 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]; }

void EvtVubHybrid::init (   )  [virtual]

Reimplemented from EvtDecayBase. { // 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); }

void EvtVubHybrid::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; }

void EvtVubHybrid::readWeights (  int  startArg = 0  )

{ _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; } }

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 EvtVubHybrid::_a [private]

double EvtVubHybrid::_alphas [private]

double* EvtVubHybrid::_bins_El [private]

double* EvtVubHybrid::_bins_mX [private]

double* EvtVubHybrid::_bins_q2 [private]

bool EvtDecayBase::_daugsDecayedByParentModel [protected, inherited]

EvtVubdGamma* EvtVubHybrid::_dGamma [private]

double EvtVubHybrid::_dGMax [private]

double EvtVubHybrid::_masscut [private]

double EvtVubHybrid::_mb [private]

int EvtVubHybrid::_nbins [private]

int EvtVubHybrid::_nbins_El [private]

int EvtVubHybrid::_nbins_mX [private]

int EvtVubHybrid::_nbins_q2 [private]

bool EvtVubHybrid::_noHybrid [private]

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

bool EvtVubHybrid::_storeQplus [private]

double* EvtVubHybrid::_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/EvtVubHybrid.hh
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/src/EvtGen/EvtGenModels/EvtVubHybrid.cc

---