TEvtGen/EvtGenModels/EvtVub.cxx

   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(mB*mB/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 = p2*pow(mB-qplus,2)+2*qplus*(Eh-qplus)+qplus*qplus;
 228           if ( sh > _masses[0]*_masses[0]
 229                && mB*mB + sh - 2*mB*Eh > ml*ml) {
 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,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(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   = mB*mB + sh - 2*mB*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(ml*ml + ptmp*ptmp);
 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] = sttmp*cos(phL)*ptmp*xW[j]
 336       +        sttmp*sin(phL)*ptmp*yW[j]
 337       +          ctL         *ptmp*zW[j];
 338
 339   double apLW = ptmp;
 340
 341   // boost them back in the B Meson restframe
 342   double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
 343
 344   ptmp = sqrt(El*El-ml*ml);
 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] + (ctLL*ptmp - ctL*apLW)/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 }