TEvtGen/EvtGenModels/EvtBtoXsll.cpp

   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 // Module: EvtBtoXsll.cc
   9 //
  10 // Description: Routine to generate non-resonant B -> Xs l+ l- decays.
  11 // It generates a dilepton mass spectrum according to Kruger and Sehgal
  12 // and then generates the two lepton momenta accoring to Ali et al.
  13 // The resultant X_s particles may be decayed by JETSET.
  14 //
  15 // Modification history:
  16 //
  17 //    Stephane Willocq     Jan  17, 2001    Module created
  18 //    Stephane Willocq     Jul  15, 2003    Input model parameters
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23
  24 #include <stdlib.h>
  25 #include "EvtGenBase/EvtRandom.hh"
  26 #include "EvtGenBase/EvtParticle.hh"
  27 #include "EvtGenBase/EvtGenKine.hh"
  28 #include "EvtGenBase/EvtPDL.hh"
  29 #include "EvtGenBase/EvtReport.hh"
  30 #include "EvtGenModels/EvtbTosllAmp.hh"
  31 #include "EvtGenModels/EvtBtoXsll.hh"
  32 #include "EvtGenModels/EvtBtoXsllUtil.hh"
  33 #include "EvtGenBase/EvtConst.hh"
  34 #include "EvtGenBase/EvtId.hh"
  35 using std::endl;
  36
  37 EvtBtoXsll::~EvtBtoXsll() {
  38   delete _calcprob;
  39 }
  40
  41 std::string EvtBtoXsll::getName(){
  42
  43   return "BTOXSLL";
  44
  45 }
  46
  47 EvtDecayBase* EvtBtoXsll::clone(){
  48
  49   return new EvtBtoXsll;
  50
  51 }
  52
  53
  54 void EvtBtoXsll::init(){
  55
  56   // check that there are no arguments
  57
  58   checkNArg(0,4,5);
  59
  60   checkNDaug(3);
  61
  62   // Check that the two leptons are the same type
  63
  64   EvtId lepton1type = getDaug(1);
  65   EvtId lepton2type = getDaug(2);
  66
  67   int etyp   = 0;
  68   int mutyp  = 0;
  69   int tautyp = 0;
  70   if ( lepton1type == EvtPDL::getId("e+") ||
  71        lepton1type == EvtPDL::getId("e-") ) { etyp++;}
  72   if ( lepton2type == EvtPDL::getId("e+") ||
  73        lepton2type == EvtPDL::getId("e-") ) { etyp++;}
  74   if ( lepton1type == EvtPDL::getId("mu+") ||
  75        lepton1type == EvtPDL::getId("mu-") ) { mutyp++;}
  76   if ( lepton2type == EvtPDL::getId("mu+") ||
  77        lepton2type == EvtPDL::getId("mu-") ) { mutyp++;}
  78   if ( lepton1type == EvtPDL::getId("tau+") ||
  79        lepton1type == EvtPDL::getId("tau-") ) { tautyp++;}
  80   if ( lepton2type == EvtPDL::getId("tau+") ||
  81        lepton2type == EvtPDL::getId("tau-") ) { tautyp++;}
  82
  83   if ( etyp != 2 && mutyp != 2 && tautyp != 2 ) {
  84
  85      report(ERROR,"EvtGen") << "Expect two leptons of the same type in EvtBtoXsll.cc\n";
  86      ::abort();
  87   }
  88
  89   // Check that the second and third entries are leptons with positive
  90   // and negative charge, respectively
  91
  92   int lpos = 0;
  93   int lneg = 0;
  94   if ( lepton1type == EvtPDL::getId("e+")  ||
  95        lepton1type == EvtPDL::getId("mu+") ||
  96        lepton1type == EvtPDL::getId("tau+") ) { lpos++;}
  97   if ( lepton2type == EvtPDL::getId("e-")  ||
  98        lepton2type == EvtPDL::getId("mu-") ||
  99        lepton2type == EvtPDL::getId("tau-") ) { lneg++;}
 100
 101   if ( lpos != 1 || lneg != 1 ) {
 102
 103      report(ERROR,"EvtGen") << "Expect 2nd and 3rd particles to be positive and negative leptons in EvtBtoXsll.cc\n";
 104      ::abort();
 105   }
 106
 107
 108   _mb=4.8;
 109   _ms=0.2;
 110   _mq=0.;
 111   _pf=0.41;
 112   _mxmin=1.1;
 113   if ( getNArg()==4)
 114   {
 115     // b-quark mass
 116     _mb = getArg(0);
 117     // s-quark mass
 118     _ms = getArg(1);
 119     // spectator quark mass
 120     _mq = getArg(2);
 121     // Fermi motion parameter
 122     _pf = getArg(3);
 123   }
 124   if ( getNArg()==5)
 125   {
 126     _mxmin = getArg(4);
 127   }
 128
 129   _calcprob = new EvtBtoXsllUtil;
 130
 131   double ml = EvtPDL::getMeanMass(getDaug(1));
 132
 133   // determine the maximum probability density from dGdsProb
 134
 135   int i, j;
 136   int nsteps  = 100;
 137   double s    = 0.0;
 138   double smin = 4.0 * ml * ml;
 139   double smax = (_mb - _ms)*(_mb - _ms);
 140   double probMax  = -10000.0;
 141   double sProbMax = -10.0;
 142   double uProbMax = -10.0;
 143
 144   for (i=0;i<nsteps;i++)
 145   {
 146     s = smin + (i+0.002)*(smax - smin)/(double)nsteps;
 147     double prob = _calcprob->dGdsProb(_mb, _ms, ml, s);
 148     if (prob > probMax)
 149     {
 150       sProbMax = s;
 151       probMax  = prob;
 152     }
 153   }
 154
 155   _dGdsProbMax = probMax;
 156
 157   if ( verbose() ) {
 158     report(INFO,"EvtGen") << "dGdsProbMax = " << probMax << " for s = "  << sProbMax << endl;
 159   }
 160
 161   // determine the maximum probability density from dGdsdupProb
 162
 163   probMax  = -10000.0;
 164   sProbMax = -10.0;
 165
 166   for (i=0;i<nsteps;i++)
 167   {
 168     s = smin + (i+0.002)*(smax - smin)/(double)nsteps;
 169     double umax = sqrt((s - (_mb + _ms)*(_mb + _ms)) * (s - (_mb - _ms)*(_mb - _ms)));
 170     for (j=0;j<nsteps;j++)
 171     {
 172       double u = -umax + (j+0.002)*(2.0*umax)/(double)nsteps;
 173       double prob = _calcprob->dGdsdupProb(_mb, _ms, ml, s, u);
 174       if (prob > probMax)
 175       {
 176         sProbMax = s;
 177         uProbMax = u;
 178         probMax  = prob;
 179       }
 180     }
 181   }
 182
 183   _dGdsdupProbMax = 2.0*probMax;
 184
 185   if ( verbose() ) {
 186    report(INFO,"EvtGen") << "dGdsdupProbMax = " << probMax << " for s = "  << sProbMax
 187                           << " and u = " << uProbMax << endl;
 188   }
 189
 190 }
 191
 192 void EvtBtoXsll::initProbMax(){
 193
 194   noProbMax();
 195
 196 }
 197
 198 void EvtBtoXsll::decay( EvtParticle *p ){
 199
 200   p->makeDaughters(getNDaug(),getDaugs());
 201
 202   EvtParticle* xhadron = p->getDaug(0);
 203   EvtParticle* leptonp = p->getDaug(1);
 204   EvtParticle* leptonn = p->getDaug(2);
 205
 206   double mass[3];
 207
 208   findMasses( p, getNDaug(), getDaugs(), mass );
 209
 210   double mB = p->mass();
 211   double ml = mass[1];
 212   double pb;
 213
 214   int im = 0;
 215   static int nmsg = 0;
 216   double xhadronMass = -999.0;
 217
 218   EvtVector4R p4xhadron;
 219   EvtVector4R p4leptonp;
 220   EvtVector4R p4leptonn;
 221
 222   // require the hadronic system has mass greater than that of a Kaon pion pair
 223
 224   //  while (xhadronMass < 0.6333)
 225   // the above minimum value of K+pi mass appears to be too close
 226   // to threshold as far as JETSET is concerned
 227   // (JETSET gets caught in an infinite loop)
 228   // so we choose a lightly larger value for the threshold
 229   while (xhadronMass < _mxmin)
 230   {
 231     im++;
 232
 233     // Apply Fermi motion and determine effective b-quark mass
 234
 235     // Old BaBar MC parameters
 236     //    double pf = 0.25;
 237     //    double ms = 0.2;
 238     //    double mq = 0.3;
 239
 240     double mb = 0.0;
 241
 242     double xbox, ybox;
 243
 244     while (mb <= 0.0)
 245     {
 246       pb = _calcprob->FermiMomentum(_pf);
 247
 248       // effective b-quark mass
 249       mb = mB*mB + _mq*_mq - 2.0*mB*sqrt(pb*pb + _mq*_mq);
 250       if ( mb>0. && sqrt(mb)-_ms  < 2.0*ml ) mb= -10.;
 251     }
 252     mb = sqrt(mb);
 253
 254     //    cout << "b-quark momentum = " << pb << " mass = " <<  mb << endl;
 255
 256     // generate a dilepton invariant mass
 257
 258     double s    = 0.0;
 259     double smin = 4.0 * ml * ml;
 260     double smax = (mb - _ms)*(mb - _ms);
 261
 262     while (s == 0.0)
 263       {
 264       xbox = EvtRandom::Flat(smin, smax);
 265       ybox = EvtRandom::Flat(_dGdsProbMax);
 266       double prob= _calcprob->dGdsProb(mb, _ms, ml, xbox);
 267       if ( !(prob>=0.0) && !(prob<=0.0)) {
 268         //      report(INFO,"EvtGen") << "nan from dGdsProb " << prob << " " << mb << " " << _ms << " " << ml << " " << xbox << std::endl;
 269       }
 270       if ( ybox < prob ) s=xbox;
 271       }
 272
 273     //    cout << "dGdsProb(s) = " << _calcprob->dGdsProb(mb, _ms, ml, s)
 274     //         << " for s = " << s << endl;
 275
 276     // two-body decay of b quark at rest into s quark and dilepton pair:
 277     // b -> s (ll)
 278
 279     EvtVector4R p4sdilep[2];
 280
 281     double msdilep[2];
 282     msdilep[0] = _ms;
 283     msdilep[1] = sqrt(s);
 284
 285     EvtGenKine::PhaseSpace(2, msdilep, p4sdilep, mb);
 286
 287     // generate dilepton decay with the expected asymmetry: (ll) -> l+ l-
 288
 289     EvtVector4R p4ll[2];
 290
 291     double mll[2];
 292     mll[0] = ml;
 293     mll[1] = ml;
 294
 295     double tmp = 0.0;
 296
 297     while (tmp == 0.0)
 298     {
 299       // (ll) -> l+ l- decay in dilepton rest frame
 300
 301       EvtGenKine::PhaseSpace(2, mll, p4ll, msdilep[1]);
 302
 303       // boost to b-quark rest frame
 304
 305       p4ll[0] = boostTo(p4ll[0], p4sdilep[1]);
 306       p4ll[1] = boostTo(p4ll[1], p4sdilep[1]);
 307
 308       // compute kinematical variable u
 309
 310       EvtVector4R p4slp = p4sdilep[0] + p4ll[0];
 311       EvtVector4R p4sln = p4sdilep[0] + p4ll[1];
 312
 313       double u = p4slp.mass2() - p4sln.mass2();
 314
 315       ybox = EvtRandom::Flat(_dGdsdupProbMax);
 316
 317       double prob = _calcprob->dGdsdupProb(mb, _ms, ml, s, u);
 318       if ( !(prob>=0.0) && !(prob<=0.0)) {
 319         report(INFO,"EvtGen") << "nan from dGdsProb " << prob << " " << mb << " " << _ms << " " << ml << " " << s << " " << u << std::endl;
 320       }
 321       if (prob > _dGdsdupProbMax && nmsg < 20)
 322       {
 323         report(INFO,"EvtGen") << "d2gdsdup GT d2gdsdup_max:" << prob
 324              << " " << _dGdsdupProbMax
 325              << " for s = " << s << " u = " << u << " mb = " << mb << endl;
 326                  nmsg++;
 327       }
 328       if (ybox < prob)
 329         {
 330           tmp = 1.0;
 331           //        cout << "dGdsdupProb(s) = " << prob
 332           //             << " for u = " << u << endl;
 333         }
 334     }
 335
 336
 337     // assign 4-momenta to valence quarks inside B meson in B rest frame
 338
 339     double phi   = EvtRandom::Flat( EvtConst::twoPi );
 340     double costh = EvtRandom::Flat( -1.0, 1.0 );
 341     double sinth = sqrt(1.0 - costh*costh);
 342
 343     // b-quark four-momentum in B meson rest frame
 344
 345     EvtVector4R p4b(sqrt(mb*mb + pb*pb),
 346                     pb*sinth*sin(phi),
 347                     pb*sinth*cos(phi),
 348                     pb*costh);
 349
 350     // B meson in its rest frame
 351     //
 352     //    EvtVector4R p4B(mB, 0.0, 0.0, 0.0);
 353     //
 354     // boost B meson to b-quark rest frame
 355     //
 356     //    p4B = boostTo(p4B, p4b);
 357     //
 358     //    cout << " B meson mass in b-quark rest frame = " << p4B.mass() << endl;
 359
 360     // boost s, l+ and l- to B meson rest frame
 361
 362     //    EvtVector4R p4s = boostTo(p4sdilep[0], p4B);
 363     //    p4leptonp       = boostTo(p4ll[0],     p4B);
 364     //    p4leptonn       = boostTo(p4ll[1],     p4B);
 365
 366     EvtVector4R p4s = boostTo(p4sdilep[0], p4b);
 367     p4leptonp       = boostTo(p4ll[0],     p4b);
 368     p4leptonn       = boostTo(p4ll[1],     p4b);
 369
 370     // spectator quark in B meson rest frame
 371
 372     EvtVector4R p4q( sqrt(pb*pb + _mq*_mq), -p4b.get(1), -p4b.get(2), -p4b.get(3) );
 373
 374     // hadron system in B meson rest frame
 375
 376     p4xhadron = p4s + p4q;
 377     xhadronMass = p4xhadron.mass();
 378
 379     //    cout << "Xs mass = " << xhadronMass << " trial " << im << endl;
 380   }
 381
 382   // initialize the decay products
 383
 384   xhadron->init(getDaug(0), p4xhadron);
 385
 386   // For B-bar mesons (i.e. containing a b quark) we have the normal
 387   // order of leptons
 388   if ( p->getId() == EvtPDL::getId("anti-B0") ||
 389        p->getId() == EvtPDL::getId("B-") )
 390   {
 391     leptonp->init(getDaug(1), p4leptonp);
 392     leptonn->init(getDaug(2), p4leptonn);
 393   }
 394   // For B mesons (i.e. containing a b-bar quark) we need to flip the
 395   // role of the positive and negative leptons in order to produce the
 396   // correct forward-backward asymmetry between the two leptons
 397   else
 398   {
 399     leptonp->init(getDaug(1), p4leptonn);
 400     leptonn->init(getDaug(2), p4leptonp);
 401   }
 402
 403   return ;
 404 }