[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtVub.cxx

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information:
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtVub.cc
//
// Description: Routine to decay a particle according th phase space 
//
// Modification history:
//
//    Sven Menke       January 17, 2001       Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtVub.hh"
#include <string>
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenModels/EvtPFermi.hh"
#include "EvtGenModels/EvtVubdGamma.hh"
#include "EvtGenBase/EvtRandom.hh"
using std::endl;

EvtVub::~EvtVub() {
  if (_dGamma) delete _dGamma;
  if (_masses) delete [] _masses;
  if (_weights) delete [] _weights;
}

std::string EvtVub::getName(){

  return "VUB";     

}

EvtDecayBase* EvtVub::clone(){

  return new EvtVub;

}


void EvtVub::init(){

  // 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 (size_t index=0; index<_pf.size(); index++) {
    _pf[index]/=_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(){

  noProbMax();

}

void EvtVub::decay( EvtParticle *p ){

  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.0,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*EvtConst::pi);
  double phL = EvtRandom::Flat(0,2*EvtConst::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;
    
  // 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 ;
}


double EvtVub::findPFermi() {

  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;
  }
  
  double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
  
  return (nBinsBelow + bFract) * oOverBins;

}
Commit	Line	Data
da0e9ce3	1	//--------------------------------------------------------------------------
	2	//
	3	// Environment:
	4	// This software is part of the EvtGen package developed jointly
	5	// for the BaBar and CLEO collaborations. If you use all or part
	6	// of it, please give an appropriate acknowledgement.
	7	//
	8	// Copyright Information:
	9	// Copyright (C) 1998 Caltech, UCSB
	10	//
	11	// Module: EvtVub.cc
	12	//
	13	// Description: Routine to decay a particle according th phase space
	14	//
	15	// Modification history:
	16	//
	17	// Sven Menke January 17, 2001 Module created
	18	//
	19	//------------------------------------------------------------------------
	20	//
	21	#include "EvtGenBase/EvtPatches.hh"
	22	#include <stdlib.h>
	23	#include "EvtGenBase/EvtParticle.hh"
	24	#include "EvtGenBase/EvtGenKine.hh"
	25	#include "EvtGenBase/EvtPDL.hh"
	26	#include "EvtGenBase/EvtReport.hh"
	27	#include "EvtGenModels/EvtVub.hh"
	28	#include <string>
	29	#include "EvtGenBase/EvtVector4R.hh"
	30	#include "EvtGenModels/EvtPFermi.hh"
	31	#include "EvtGenModels/EvtVubdGamma.hh"
	32	#include "EvtGenBase/EvtRandom.hh"
	33	using std::endl;
	34
	35	EvtVub::~EvtVub() {
	36	if (_dGamma) delete _dGamma;
	37	if (_masses) delete [] _masses;
	38	if (_weights) delete [] _weights;
	39	}
	40
	41	std::string EvtVub::getName(){
	42
	43	return "VUB";
	44
	45	}
	46
	47	EvtDecayBase* EvtVub::clone(){
	48
	49	return new EvtVub;
	50
	51	}
	52
	53
	54	void EvtVub::init(){
	55
	56	// check that there are at least 6 arguments
	57
	58	if (getNArg()<6) {
	59
	60	report(ERROR,"EvtGen") << "EvtVub generator expected "
	61	<< " at least 6 arguments (mb,a,alpha_s,Nbins,m1,w1,...) but found: "
	62	<<getNArg()<<endl;
	63	report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
	64	::abort();
65
66	}
67
68	_mb = getArg(0);
69	_a = getArg(1);
70	_alphas = getArg(2);
71	_nbins = abs((int)getArg(3));
72	_storeQplus = (getArg(3)<0?1:0);
73	_masses = new double[_nbins];
74	_weights = new double[_nbins];
75
76	if (getNArg()-4 != 2*_nbins) {
77	report(ERROR,"EvtGen") << "EvtVub generator expected "
78	<< _nbins << " masses and weights but found: "
79	<<(getNArg()-4)/2 <<endl;
80	report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
81	::abort();
82	}
83	int i,j = 4;
84	double maxw = 0.;
85	for (i=0;i<_nbins;i++) {
86	_masses[i] = getArg(j++);
87	if (i>0 && _masses[i] <= _masses[i-1]) {
88	report(ERROR,"EvtGen") << "EvtVub generator expected "
89	<< " mass bins in ascending order!"
90	<< "Will terminate execution!"<<endl;
91	::abort();
92	}
93	_weights[i] = getArg(j++);
94	if (_weights[i] < 0) {
95	report(ERROR,"EvtGen") << "EvtVub generator expected "
96	<< " weights >= 0, but found: "
97	<<_weights[i] <<endl;
98	report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
99	::abort();
100	}
101	if ( _weights[i] > maxw ) maxw = _weights[i];
102	}
103	if (maxw == 0) {
104	report(ERROR,"EvtGen") << "EvtVub generator expected at least one "
105	<< " weight > 0, but found none! "
106	<< "Will terminate execution!"<<endl;
107	::abort();
108	}
109	for (i=0;i<_nbins;i++) _weights[i]/=maxw;
110
111	// the maximum dGamma*p2 value depends on alpha_s only:
112
113	const double dGMax0 = 3.;
114	_dGMax = 0.21344+8.905*_alphas;
115	if ( _dGMax < dGMax0 ) _dGMax = dGMax0;
116
117	// for the Fermi Motion we need a B-Meson mass - but it's not critical
118	// to get an exact value; in order to stay in the phase space for
119	// B+- and B0 use the smaller mass
120
121	EvtId BP=EvtPDL::getId("B+");
122	EvtId B0=EvtPDL::getId("B0");
123
124	double mB0 = EvtPDL::getMaxMass(B0);
125	double mBP = EvtPDL::getMaxMass(BP);
126
127	double mB = (mB0<mBP?mB0:mBP);
128
129	const double xlow = -_mb;
130	const double xhigh = mB-_mb;
131	const int aSize = 10000;
132
133	EvtPFermi pFermi(_a,mB,_mb);
134	// pf is the cumulative distribution
135	// normalized to 1.
136	_pf.resize(aSize);
137	for(i=0;i<aSize;i++){
138	double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
139	if ( i== 0 )
140	_pf[i] = pFermi.getFPFermi(kplus);
141	else
142	_pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus);
143	}
144	for (size_t index=0; index<_pf.size(); index++) {
145	_pf[index]/=_pf[_pf.size()-1];
146	}
147
148	// static EvtHepRandomEngine myEngine;
149
150	// _pFermi = new RandGeneral(myEngine,pf,aSize,0);
151	_dGamma = new EvtVubdGamma(_alphas);
152
153	// check that there are 3 daughters
154	checkNDaug(3);
155	}
156
157	void EvtVub::initProbMax(){
158
159	noProbMax();
160
161	}
162
163	void EvtVub::decay( EvtParticle *p ){
164
165	int j;
166	// B+ -> u-bar specflav l+ nu
167
168	EvtParticle xuhad, lepton, *neutrino;
169	EvtVector4R p4;
170	// R. Faccini 21/02/03
171	// move the reweighting up , before also shooting the fermi distribution
172	double x,z,p2;
173	double sh=0.0;
174	double mB,ml,xlow,xhigh,qplus;
175	double El=0.0;
176	double Eh=0.0;
177	double kplus;
178	const double lp2epsilon=-10;
179	bool rew(true);
180	while(rew){
181
182	p->initializePhaseSpace(getNDaug(),getDaugs());
183
184	xuhad=p->getDaug(0);
185	lepton=p->getDaug(1);
186	neutrino=p->getDaug(2);
187
188	mB = p->mass();
189	ml = lepton->mass();
190
191	xlow = -_mb;
192	xhigh = mB-_mb;
193
194
195	// Fermi motion does not need to be computed inside the
196	// tryit loop as m_b in Gamma0 does not need to be replaced by (m_b+kplus).
197	// The difference however should be of the Order (lambda/m_b)^2 which is
198	// beyond the considered orders in the paper anyway ...
199
200	// for alpha_S = 0 and a mass cut on X_u not all values of kplus are
201	// possible. The maximum value is mB/2-_mb + sqrt(mB^2/4-_masses[0]^2)
202	kplus = 2*xhigh;
203
204	while( kplus >= xhigh \|\| kplus <= xlow
205	\|\| (_alphas == 0 && kplus >= mB/2-_mb
206	+ sqrt(mBmB/4-_masses[0]_masses[0]))) {
207	kplus = findPFermi(); //_pFermi->shoot();
208	kplus = xlow + kplus*(xhigh-xlow);
209	}
210	qplus = mB-_mb-kplus;
211	if( (mB-qplus)/2.<=ml)continue;
212
213	int tryit = 1;
214	while (tryit) {
215
216	x = EvtRandom::Flat();
217	z = EvtRandom::Flat(0,2);
218	p2=EvtRandom::Flat();
219	p2 = pow(10.0,lp2epsilon*p2);
220
221	El = x*(mB-qplus)/2;
222	if ( El > ml && El < mB/2) {
223
224	Eh = z*(mB-qplus)/2+qplus;
225	if ( Eh > 0 && Eh < mB ) {
226
227	sh = p2pow(mB-qplus,2)+2qplus(Eh-qplus)+qplusqplus;
228	if ( sh > _masses[0]*_masses[0]
229	&& mBmB + sh - 2mBEh > mlml) {
230
231	double xran = EvtRandom::Flat();
232
233	double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2;
234
235	if ( y > 1 ) report(WARNING,"EvtGen")<<"EvtVub decay probability > 1 found: " << y << endl;
236	if ( y >= xran ) tryit = 0;
237	}
238	}
239	}
240	}
241	// reweight the Mx distribution
242	if(_nbins>0){
243	double xran1 = EvtRandom::Flat();
244	double m = sqrt(sh);j=0;
245	while ( j < _nbins && m > _masses[j] ) j++;
246	double w = _weights[j-1];
247	if ( w >= xran1 ) rew = false;
248	} else {
249	rew = false;
250	}
251
252	}
253
254	// o.k. we have the three kineamtic variables
255	// now calculate a flat cos Theta_H [-1,1] distribution of the
256	// hadron flight direction w.r.t the B flight direction
257	// because the B is a scalar and should decay isotropic.
258	// Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
259	// and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
260	// W flight direction.
261
262	double ctH = EvtRandom::Flat(-1,1);
263	double phH = EvtRandom::Flat(0,2*EvtConst::pi);
264	double phL = EvtRandom::Flat(0,2*EvtConst::pi);
265
266	// now compute the four vectors in the B Meson restframe
267
268	double ptmp,sttmp;
269	// calculate the hadron 4 vector in the B Meson restframe
270
271	sttmp = sqrt(1-ctH*ctH);
272	ptmp = sqrt(Eh*Eh-sh);
273	double pHB[4] = {Eh,ptmpsttmpcos(phH),ptmpsttmpsin(phH),ptmp*ctH};
274	p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
275	xuhad->init( getDaug(0), p4);
276
277	if (_storeQplus ) {
278	// cludge to store the hidden parameter q+ with the decay;
279	// the lifetime of the Xu is abused for this purpose.
280	// tau = 1 ps corresponds to ctau = 0.3 mm -> in order to
281	// stay well below BaBars sensitivity we take q+/(10000 GeV) which
282	// goes up to 0.0005 in the most extreme cases as ctau in mm.
283	// To extract q+ back from the StdHepTrk its necessary to get
284	// delta_ctau = Xu->anyDaughter->getVertexTime()-Xu->getVertexTime()
285	// where these pseudo calls refere to the StdHep time stored at
286	// the production vertex in the lab for each particle. The boost
287	// has to be reversed and the result is:
288	//
289	// q+ = delta_ctau * 10000 GeV/mm * Mass_Xu/Energy_Xu
290	//
291	xuhad->setLifetime(qplus/10000.);
292	}
293
294	// calculate the W 4 vector in the B Meson restrframe
295
296	double apWB = ptmp;
297	double pWB[4] = {mB-Eh,-pHB[1],-pHB[2],-pHB[3]};
298
299	// first go in the W restframe and calculate the lepton and
300	// the neutrino in the W frame
301
302	double mW2 = mBmB + sh - 2mB*Eh;
303	double beta = ptmp/pWB[0];
304	double gamma = pWB[0]/sqrt(mW2);
305
306	double pLW[4];
307
308	ptmp = (mW2-ml*ml)/2/sqrt(mW2);
309	pLW[0] = sqrt(mlml + ptmpptmp);
310
311	double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
312	if ( ctL < -1 ) ctL = -1;
313	if ( ctL > 1 ) ctL = 1;
314	sttmp = sqrt(1-ctL*ctL);
315
316	// eX' = eZ x eW
317	double xW[3] = {-pWB[2],pWB[1],0};
318	// eZ' = eW
319	double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
320
321	double lx = sqrt(xW[0]xW[0]+xW[1]xW[1]);
322	for (j=0;j<2;j++)
323	xW[j] /= lx;
324
325	// eY' = eZ' x eX'
326	double yW[3] = {-pWB[1]pWB[3],-pWB[2]pWB[3],pWB[1]pWB[1]+pWB[2]pWB[2]};
327	double ly = sqrt(yW[0]yW[0]+yW[1]yW[1]+yW[2]*yW[2]);
328	for (j=0;j<3;j++)
329	yW[j] /= ly;
330
331	// p_lep = \|p_lep\| * ( sin(Theta) * cos(Phi) * eX'
332	// + sin(Theta) * sin(Phi) * eY'
333	// + cos(Theta) * eZ')
334	for (j=0;j<3;j++)
335	pLW[j+1] = sttmpcos(phL)ptmp*xW[j]
336	+ sttmpsin(phL)ptmp*yW[j]
337	+ ctL ptmpzW[j];
338
339	double apLW = ptmp;
340
341	// boost them back in the B Meson restframe
342	double appLB = betagammapLW[0] + gammactLapLW;
343
344	ptmp = sqrt(ElEl-mlml);
345	double ctLL = appLB/ptmp;
346
347	if ( ctLL > 1 ) ctLL = 1;
348	if ( ctLL < -1 ) ctLL = -1;
349
350	double pLB[4] = {El,0,0,0};
351	double pNB[4] = {pWB[0]-El,0,0,0};
352
353	for (j=1;j<4;j++) {
354	pLB[j] = pLW[j] + (ctLLptmp - ctLapLW)/apWB*pWB[j];
355	pNB[j] = pWB[j] - pLB[j];
356	}
357
358	p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
359	lepton->init( getDaug(1), p4);
360
361	p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
362	neutrino->init( getDaug(2), p4);
363
364	return ;
365	}
366
367
368	double EvtVub::findPFermi() {
369
370	double ranNum=EvtRandom::Flat();
371	double oOverBins= 1.0/(float(_pf.size()));
372	int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand
373	int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand
374	int middle;
375
376	while (nBinsAbove > nBinsBelow+1) {
377	middle = (nBinsAbove + nBinsBelow+1)>>1;
378	if (ranNum >= _pf[middle]) {
379	nBinsBelow = middle;
380	} else {
381	nBinsAbove = middle;
382	}
383	}
384
385	double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
386	// binMeasure is always aProbFunc[nBinsBelow],
387
388	if ( bSize == 0 ) {
389	// rand lies right in a bin of measure 0. Simply return the center
390	// of the range of that bin. (Any value between k/N and (k+1)/N is
391	// equally good, in this rare case.)
392	return (nBinsBelow + .5) * oOverBins;
393	}
394
395	double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
396
397	return (nBinsBelow + bFract) * oOverBins;
398
399	}