1 //---------------------------------------------------------------------------
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.
8 // Copyright Information:
9 // Copyright (C) 1998 Caltech, UCSB
11 // Module: EvtVubHybrid.cc
13 // Description: Routine to decay a particle according to phase space.
15 // Modification history:
17 // Jochen Dingfelder February 1, 2005 Created Module as update of the
18 // original module EvtVub by Sven Menke
19 //---------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
23 #include "EvtGenBase/EvtParticle.hh"
24 #include "EvtGenBase/EvtGenKine.hh"
25 #include "EvtGenBase/EvtPDL.hh"
26 #include "EvtGenBase/EvtReport.hh"
27 #include "EvtGenModels/EvtVubHybrid.hh"
28 #include "EvtGenBase/EvtVector4R.hh"
29 #include "EvtGenModels/EvtPFermi.hh"
30 #include "EvtGenModels/EvtVubdGamma.hh"
31 #include "EvtGenBase/EvtRandom.hh"
41 // _noHybrid will be set TRUE if the DECAY.DEC file has no binning or weights
42 // _storeQplus should alwasy be TRUE: writes out Fermi motion parameter
44 EvtVubHybrid::EvtVubHybrid()
45 : _noHybrid(false), _storeQplus(true),
46 _mb(4.62), _a(2.27), _alphas(0.22), _dGMax(3.),
47 _nbins_mX(0), _nbins_q2(0), _nbins_El(0), _nbins(0),
48 _masscut(0.28), _bins_mX(0), _bins_q2(0), _bins_El(0),
49 _weights(0), _dGamma(0)
52 EvtVubHybrid::~EvtVubHybrid() {
61 std::string EvtVubHybrid::getName(){
67 EvtDecayBase* EvtVubHybrid::clone(){
69 return new EvtVubHybrid;
73 void EvtVubHybrid::init(){
75 // check that there are at least 3 arguments
76 if (getNArg() < EvtVubHybrid::nParameters) {
77 report(ERROR,"EvtVubHybrid") << "EvtVub generator expected "
78 << "at least " << EvtVubHybrid::nParameters
79 << " arguments but found: " << getNArg()
80 << "\nWill terminate execution!"<<endl;
82 } else if (getNArg() == EvtVubHybrid::nParameters) {
83 report(WARNING,"EvtVubHybrid") << "EvtVub: generate B -> Xu l nu events "
84 << "without using the hybrid reweighting."
87 } else if (getNArg() < EvtVubHybrid::nParameters+EvtVubHybrid::nVariables) {
88 report(ERROR,"EvtVubHybrid") << "EvtVub could not read number of bins for "
89 << "all variables used in the reweighting\n"
90 << "Will terminate execution!" << endl;
94 // check that there are 3 daughters
97 // read minimum required parameters from decay.dec
102 // the maximum dGamma*p2 value depends on alpha_s only:
103 const double dGMax0 = 3.;
104 _dGMax = 0.21344+8.905*_alphas;
105 if ( _dGMax < dGMax0 ) _dGMax = dGMax0;
107 // for the Fermi Motion we need a B-Meson mass - but it's not critical
108 // to get an exact value; in order to stay in the phase space for
109 // B+- and B0 use the smaller mass
111 static double mB0 = EvtPDL::getMaxMass(EvtPDL::getId("B0"));
112 static double mBP = EvtPDL::getMaxMass(EvtPDL::getId("B+"));
113 static double mB = (mB0<mBP?mB0:mBP);
115 const double xlow = -_mb;
116 const double xhigh = mB-_mb;
117 const int aSize = 10000;
119 EvtPFermi pFermi(_a,mB,_mb);
120 // pf is the cumulative distribution normalized to 1.
122 for(int i=0;i<aSize;i++){
123 double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
125 _pf[i] = pFermi.getFPFermi(kplus);
127 _pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus);
129 for (size_t index=0; index<_pf.size(); index++) {
130 _pf[index]/=_pf[_pf.size()-1];
133 _dGamma = new EvtVubdGamma(_alphas);
135 if (_noHybrid) return; // Without hybrid weighting, nothing else to do
137 _nbins_mX = abs((int)getArg(3));
138 _nbins_q2 = abs((int)getArg(4));
139 _nbins_El = abs((int)getArg(5));
141 int nextArg = EvtVubHybrid::nParameters + EvtVubHybrid::nVariables;
143 _nbins = _nbins_mX*_nbins_q2*_nbins_El; // Binning of weight table
145 int expectArgs = nextArg + _nbins_mX +_nbins_q2 + _nbins_El + _nbins;
147 if (getNArg() < expectArgs) {
148 report(ERROR,"EvtVubHybrid")
149 << " finds " << getNArg() << " arguments, expected " << expectArgs
150 << ". Something is wrong with the tables of weights or thresholds."
151 << "\nWill terminate execution!" << endl;
155 // read bin boundaries from decay.dec
158 _bins_mX = new double[_nbins_mX];
159 for (i = 0; i < _nbins_mX; i++,nextArg++) {
160 _bins_mX[i] = getArg(nextArg);
162 _masscut = _bins_mX[0];
164 _bins_q2 = new double[_nbins_q2];
165 for (i = 0; i < _nbins_q2; i++,nextArg++) {
166 _bins_q2[i] = getArg(nextArg);
169 _bins_El = new double[_nbins_El];
170 for (i = 0; i < _nbins_El; i++,nextArg++) {
171 _bins_El[i] = getArg(nextArg);
174 // read in weights (and rescale to range 0..1)
175 readWeights(nextArg);
178 void EvtVubHybrid::initProbMax(){
182 void EvtVubHybrid::decay( EvtParticle *p ){
185 // B+ -> u-bar specflav l+ nu
187 EvtParticle *xuhad, *lepton, *neutrino;
189 // R. Faccini 21/02/03
190 // move the reweighting up , before also shooting the fermi distribution
193 double mB,ml,xlow,xhigh,qplus;
199 const double lp2epsilon=-10;
203 p->initializePhaseSpace(getNDaug(),getDaugs());
206 lepton=p->getDaug(1);
207 neutrino=p->getDaug(2);
215 // Fermi motion does not need to be computed inside the
216 // tryit loop as m_b in Gamma0 does not need to be replaced by (m_b+kplus).
217 // The difference however should be of the Order (lambda/m_b)^2 which is
218 // beyond the considered orders in the paper anyway ...
220 // for alpha_S = 0 and a mass cut on X_u not all values of kplus are
221 // possible. The maximum value is mB/2-_mb + sqrt(mB^2/4-_masscut^2)
224 while( kplus >= xhigh || kplus <= xlow
225 || (_alphas == 0 && kplus >= mB/2-_mb
226 + sqrt(mB*mB/4-_masscut*_masscut))) {
227 kplus = findPFermi(); //_pFermi->shoot();
228 kplus = xlow + kplus*(xhigh-xlow);
230 qplus = mB-_mb-kplus;
231 if( (mB-qplus)/2.<=ml) continue;
236 x = EvtRandom::Flat();
237 z = EvtRandom::Flat(0,2);
238 p2=EvtRandom::Flat();
239 p2 = pow(10,lp2epsilon*p2);
242 if ( El > ml && El < mB/2) {
244 Eh = z*(mB-qplus)/2+qplus;
245 if ( Eh > 0 && Eh < mB ) {
247 sh = p2*pow(mB-qplus,2)+2*qplus*(Eh-qplus)+qplus*qplus;
248 if ( sh > _masscut*_masscut
249 && mB*mB + sh - 2*mB*Eh > ml*ml) {
251 double xran = EvtRandom::Flat();
253 double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2;
255 if ( y > 1 ) report(WARNING,"EvtVubHybrid") <<"EvtVubHybrid decay probability > 1 found: " << y << endl;
256 if ( y >= xran ) tryit = 0;
262 // compute all kinematic variables needed for reweighting (J. Dingfelder)
264 q2 = mB*mB + sh - 2*mB*Eh;
266 // Reweighting in bins of mX, q2, El (J. Dingfelder)
268 double xran1 = EvtRandom::Flat();
270 if (!_noHybrid) w = getWeight(mX, q2, El);
271 if ( w >= xran1 ) rew = false;
278 // o.k. we have the three kineamtic variables
279 // now calculate a flat cos Theta_H [-1,1] distribution of the
280 // hadron flight direction w.r.t the B flight direction
281 // because the B is a scalar and should decay isotropic.
282 // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
283 // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
284 // W flight direction.
286 double ctH = EvtRandom::Flat(-1,1);
287 double phH = EvtRandom::Flat(0,2*M_PI);
288 double phL = EvtRandom::Flat(0,2*M_PI);
290 // now compute the four vectors in the B Meson restframe
293 // calculate the hadron 4 vector in the B Meson restframe
295 sttmp = sqrt(1-ctH*ctH);
296 ptmp = sqrt(Eh*Eh-sh);
297 double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
298 p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
299 xuhad->init( getDaug(0), p4);
302 // cludge to store the hidden parameter q+ with the decay;
303 // the lifetime of the Xu is abused for this purpose.
304 // tau = 1 ps corresponds to ctau = 0.3 mm -> in order to
305 // stay well below BaBars sensitivity we take q+/(10000 GeV) which
306 // goes up to 0.0005 in the most extreme cases as ctau in mm.
307 // To extract q+ back from the StdHepTrk its necessary to get
308 // delta_ctau = Xu->anyDaughter->getVertexTime()-Xu->getVertexTime()
309 // where these pseudo calls refere to the StdHep time stored at
310 // the production vertex in the lab for each particle. The boost
311 // has to be reversed and the result is:
313 // q+ = delta_ctau * 10000 GeV/mm * Mass_Xu/Energy_Xu
315 xuhad->setLifetime(qplus/10000.);
318 // calculate the W 4 vector in the B Meson restrframe
321 double pWB[4] = {mB-Eh,-pHB[1],-pHB[2],-pHB[3]};
323 // first go in the W restframe and calculate the lepton and
324 // the neutrino in the W frame
326 double mW2 = mB*mB + sh - 2*mB*Eh;
327 double beta = ptmp/pWB[0];
328 double gamma = pWB[0]/sqrt(mW2);
332 ptmp = (mW2-ml*ml)/2/sqrt(mW2);
333 pLW[0] = sqrt(ml*ml + ptmp*ptmp);
335 double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
336 if ( ctL < -1 ) ctL = -1;
337 if ( ctL > 1 ) ctL = 1;
338 sttmp = sqrt(1-ctL*ctL);
341 double xW[3] = {-pWB[2],pWB[1],0};
343 double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
345 double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
350 double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
351 double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
355 // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX'
356 // + sin(Theta) * sin(Phi) * eY'
357 // + cos(Theta) * eZ')
359 pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
360 + sttmp*sin(phL)*ptmp*yW[j]
364 // calculate the neutrino 4 vector in the W restframe
365 //double pNW[4] = {sqrt(mW2)-pLW[0],-pLW[1],-pLW[2],-pLW[3]};
367 // boost them back in the B Meson restframe
369 double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
371 ptmp = sqrt(El*El-ml*ml);
372 double ctLL = appLB/ptmp;
374 if ( ctLL > 1 ) ctLL = 1;
375 if ( ctLL < -1 ) ctLL = -1;
377 double pLB[4] = {El,0,0,0};
378 double pNB[4] = {pWB[0]-El,0,0,0};
381 pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
382 pNB[j] = pWB[j] - pLB[j];
385 p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
386 lepton->init( getDaug(1), p4);
388 p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
389 neutrino->init( getDaug(2), p4);
395 double EvtVubHybrid::findPFermi() {
397 double ranNum=EvtRandom::Flat();
398 double oOverBins= 1.0/(float(_pf.size()));
399 int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand
400 int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand
403 while (nBinsAbove > nBinsBelow+1) {
404 middle = (nBinsAbove + nBinsBelow+1)>>1;
405 if (ranNum >= _pf[middle]) {
412 double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
413 // binMeasure is always aProbFunc[nBinsBelow],
416 // rand lies right in a bin of measure 0. Simply return the center
417 // of the range of that bin. (Any value between k/N and (k+1)/N is
418 // equally good, in this rare case.)
419 return (nBinsBelow + .5) * oOverBins;
422 double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
424 return (nBinsBelow + bFract) * oOverBins;
429 double EvtVubHybrid::getWeight(double mX, double q2, double El) {
435 for (int i = 0; i < _nbins_mX; i++) {
436 if (mX >= _bins_mX[i]) ibin_mX = i;
438 for (int i = 0; i < _nbins_q2; i++) {
439 if (q2 >= _bins_q2[i]) ibin_q2 = i;
441 for (int i = 0; i < _nbins_El; i++) {
442 if (El >= _bins_El[i]) ibin_El = i;
444 int ibin = ibin_mX + ibin_q2*_nbins_mX + ibin_El*_nbins_mX*_nbins_q2;
446 if ( (ibin_mX < 0) || (ibin_q2 < 0) || (ibin_El < 0) ) {
447 report(ERROR,"EvtVubHybrid") << "Cannot determine hybrid weight "
449 << "-> assign weight = 0" << endl;
453 return _weights[ibin];
457 void EvtVubHybrid::readWeights(int startArg) {
458 _weights = new double[_nbins];
461 for (int i = 0; i < _nbins; i++, startArg++) {
462 _weights[i] = getArg(startArg);
463 if (_weights[i] > maxw) maxw = _weights[i];
467 report(ERROR,"EvtVubHybrid") << "EvtVub generator expected at least one "
468 << " weight > 0, but found none! "
469 << "Will terminate execution!"<<endl;
473 // rescale weights (to be in range 0..1)
474 for (int i = 0; i < _nbins; i++) {