EvtBtoXsll.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// 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.
//
// Module: EvtBtoXsll.cc
//
// Description: Routine to generate non-resonant B -> Xs l+ l- decays.
// It generates a dilepton mass spectrum according to Kruger and Sehgal
// and then generates the two lepton momenta accoring to Ali et al.
// The resultant X_s particles may be decayed by JETSET.
//
// Modification history:
//
//    Stephane Willocq     Jan  17, 2001    Module created
//    Stephane Willocq     Jul  15, 2003    Input model parameters
//
//------------------------------------------------------------------------
//
#include "../EvtGenBase/EvtPatches.hh"
 
#include "../EvtGenBase/EvtConst.hh"
#include "../EvtGenBase/EvtGenKine.hh"
#include "../EvtGenBase/EvtId.hh"
#include "../EvtGenBase/EvtPDL.hh"
#include "../EvtGenBase/EvtParticle.hh"
#include "../EvtGenBase/EvtRandom.hh"
#include "../EvtGenBase/EvtReport.hh"
#include "EvtBtoXsll.hh"
#include "EvtBtoXsllUtil.hh"
#include "EvtbTosllAmp.hh"
#include <stdlib.h>
using std::endl;
 
EvtBtoXsll::~EvtBtoXsll() {}
 
void EvtBtoXsll::getName( std::string& model_name ) { model_name = "BTOXSLL"; }
 
EvtDecayBase* EvtBtoXsll::clone() { return new EvtBtoXsll; }
 
void EvtBtoXsll::init() {
 
  // check that there are no arguments
 
  checkNArg( 0, 4, 5 );
 
  checkNDaug( 3 );
 
  // Check that the two leptons are the same type
 
  EvtId lepton1type = getDaug( 1 );
  EvtId lepton2type = getDaug( 2 );
 
  int etyp   = 0;
  int mutyp  = 0;
  int tautyp = 0;
  if ( lepton1type == EvtPDL::getId( "e+" ) || lepton1type == EvtPDL::getId( "e-" ) )
  { etyp++; }
  if ( lepton2type == EvtPDL::getId( "e+" ) || lepton2type == EvtPDL::getId( "e-" ) )
  { etyp++; }
  if ( lepton1type == EvtPDL::getId( "mu+" ) || lepton1type == EvtPDL::getId( "mu-" ) )
  { mutyp++; }
  if ( lepton2type == EvtPDL::getId( "mu+" ) || lepton2type == EvtPDL::getId( "mu-" ) )
  { mutyp++; }
  if ( lepton1type == EvtPDL::getId( "tau+" ) || lepton1type == EvtPDL::getId( "tau-" ) )
  { tautyp++; }
  if ( lepton2type == EvtPDL::getId( "tau+" ) || lepton2type == EvtPDL::getId( "tau-" ) )
  { tautyp++; }
 
  if ( etyp != 2 && mutyp != 2 && tautyp != 2 )
  {
 
    report( ERROR, "EvtGen" ) << "Expect two leptons of the same type in EvtBtoXsll.cc\n";
    ::abort();
  }
 
  // Check that the second and third entries are leptons with positive
  // and negative charge, respectively
 
  int lpos = 0;
  int lneg = 0;
  if ( lepton1type == EvtPDL::getId( "e+" ) || lepton1type == EvtPDL::getId( "mu+" ) ||
       lepton1type == EvtPDL::getId( "tau+" ) )
  { lpos++; }
  if ( lepton2type == EvtPDL::getId( "e-" ) || lepton2type == EvtPDL::getId( "mu-" ) ||
       lepton2type == EvtPDL::getId( "tau-" ) )
  { lneg++; }
 
  if ( lpos != 1 || lneg != 1 )
  {
 
    report( ERROR, "EvtGen" ) << "Expect 2nd and 3rd particles to be positive and negative "
                                 "leptons in EvtBtoXsll.cc\n";
    ::abort();
  }
 
  _mb    = 4.8;
  _ms    = 0.2;
  _mq    = 0.;
  _pf    = 0.41;
  _mxmin = 1.1;
  if ( getNArg() == 4 )
  {
    // b-quark mass
    _mb = getArg( 0 );
    // s-quark mass
    _ms = getArg( 1 );
    // spectator quark mass
    _mq = getArg( 2 );
    // Fermi motion parameter
    _pf = getArg( 3 );
  }
  if ( getNArg() == 5 ) { _mxmin = getArg( 4 ); }
 
  _calcprob = new EvtBtoXsllUtil;
 
  double ml = EvtPDL::getMeanMass( getDaug( 1 ) );
 
  // determine the maximum probability density from dGdsProb
 
  int    i, j;
  int    nsteps   = 100;
  double s        = 0.0;
  double smin     = 4.0 * ml * ml;
  double smax     = ( _mb - _ms ) * ( _mb - _ms );
  double probMax  = -10000.0;
  double sProbMax = -10.0;
  double uProbMax = -10.0;
 
  for ( i = 0; i < nsteps; i++ )
  {
    s           = smin + ( i + 0.0005 ) * ( smax - smin ) / (double)nsteps;
    double prob = _calcprob->dGdsProb( _mb, _ms, ml, s );
    if ( prob > probMax )
    {
      sProbMax = s;
      probMax  = prob;
    }
  }
 
  _dGdsProbMax = probMax;
 
  if ( verbose() )
  {
    report( INFO, "EvtGen" ) << "dGdsProbMax = " << probMax << " for s = " << sProbMax << endl;
  }
 
  // determine the maximum probability density from dGdsdupProb
 
  probMax  = -10000.0;
  sProbMax = -10.0;
 
  for ( i = 0; i < nsteps; i++ )
  {
    s = smin + ( i + 0.0005 ) * ( smax - smin ) / (double)nsteps;
    double umax =
        sqrt( ( s - ( _mb + _ms ) * ( _mb + _ms ) ) * ( s - ( _mb - _ms ) * ( _mb - _ms ) ) );
    for ( j = 0; j < nsteps; j++ )
    {
      double u    = -umax + ( j + 0.0005 ) * ( 2.0 * umax ) / (double)nsteps;
      double prob = _calcprob->dGdsdupProb( _mb, _ms, ml, s, u );
      if ( prob > probMax )
      {
        sProbMax = s;
        uProbMax = u;
        probMax  = prob;
      }
    }
  }
 
  _dGdsdupProbMax = probMax;
 
  if ( verbose() )
  {
    report( INFO, "EvtGen" ) << "dGdsdupProbMax = " << probMax << " for s = " << sProbMax
                             << " and u = " << uProbMax << endl;
  }
}
 
void EvtBtoXsll::initProbMax() { noProbMax(); }
 
void EvtBtoXsll::decay( EvtParticle* p ) {
 
  p->makeDaughters( getNDaug(), getDaugs() );
 
  EvtParticle* xhadron = p->getDaug( 0 );
  EvtParticle* leptonp = p->getDaug( 1 );
  EvtParticle* leptonn = p->getDaug( 2 );
 
  double mass[3];
 
  findMasses( p, getNDaug(), getDaugs(), mass );
 
  double mB = p->mass();
  double ml = mass[1];
  double pb;
 
  int        im          = 0;
  static int nmsg        = 0;
  double     xhadronMass = -999.0;
 
  EvtVector4R p4xhadron;
  EvtVector4R p4leptonp;
  EvtVector4R p4leptonn;
 
  // require the hadronic system has mass greater than that of a Kaon pion pair
 
  //  while (xhadronMass < 0.6333)
  // the above minimum value of K+pi mass appears to be too close
  // to threshold as far as JETSET is concerned
  // (JETSET gets caught in an infinite loop)
  // so we choose a lightly larger value for the threshold
  while ( xhadronMass < _mxmin )
  {
    im++;
 
    // Apply Fermi motion and determine effective b-quark mass
 
    // Old BaBar MC parameters
    //    double pf = 0.25;
    //    double ms = 0.2;
    //    double mq = 0.3;
 
    double mb = 0.0;
 
    double xbox, ybox;
 
    while ( mb <= 0.0 )
    {
      pb = _calcprob->FermiMomentum( _pf );
 
      // effective b-quark mass
      mb = mB * mB + _mq * _mq - 2.0 * mB * sqrt( pb * pb + _mq * _mq );
    }
    mb = sqrt( mb );
 
    //    cout << "b-quark momentum = " << pb << " mass = " <<  mb << endl;
 
    // generate a dilepton invariant mass
 
    double s    = 0.0;
    double smin = 4.0 * ml * ml;
    double smax = ( mb - _ms ) * ( mb - _ms );
 
    while ( s == 0.0 )
    {
      xbox = EvtRandom::Flat( smin, smax );
      ybox = EvtRandom::Flat( _dGdsProbMax );
      if ( ybox < _calcprob->dGdsProb( mb, _ms, ml, xbox ) ) { s = xbox; }
    }
 
    //    cout << "dGdsProb(s) = " << _calcprob->dGdsProb(mb, _ms, ml, s)
    //         << " for s = " << s << endl;
 
    // two-body decay of b quark at rest into s quark and dilepton pair:
    // b -> s (ll)
 
    EvtVector4R p4sdilep[2];
 
    double msdilep[2];
    msdilep[0] = _ms;
    msdilep[1] = sqrt( s );
 
    EvtGenKine::PhaseSpace( 2, msdilep, p4sdilep, mb );
 
    // generate dilepton decay with the expected asymmetry: (ll) -> l+ l-
 
    EvtVector4R p4ll[2];
 
    double mll[2];
    mll[0] = ml;
    mll[1] = ml;
 
    double tmp = 0.0;
 
    while ( tmp == 0.0 )
    {
      // (ll) -> l+ l- decay in dilepton rest frame
 
      EvtGenKine::PhaseSpace( 2, mll, p4ll, msdilep[1] );
 
      // boost to b-quark rest frame
 
      p4ll[0] = boostTo( p4ll[0], p4sdilep[1] );
      p4ll[1] = boostTo( p4ll[1], p4sdilep[1] );
 
      // compute kinematical variable u
 
      EvtVector4R p4slp = p4sdilep[0] + p4ll[0];
      EvtVector4R p4sln = p4sdilep[0] + p4ll[1];
 
      double u = p4slp.mass2() - p4sln.mass2();
 
      ybox = EvtRandom::Flat( _dGdsdupProbMax );
 
      double prob = _calcprob->dGdsdupProb( mb, _ms, ml, s, u );
      if ( prob > _dGdsdupProbMax && nmsg < 20 )
      {
        report( INFO, "EvtGen" )
            << "d2gdsdup GT d2gdsdup_max:" << prob << " " << _dGdsdupProbMax
            << " for s = " << s << " u = " << u << " mb = " << mb << endl;
        nmsg++;
      }
      if ( ybox < prob )
      {
        tmp = 1.0;
        //        cout << "dGdsdupProb(s) = " << prob
        //             << " for u = " << u << endl;
      }
    }
 
    // assign 4-momenta to valence quarks inside B meson in B rest frame
 
    double phi   = EvtRandom::Flat( EvtConst::twoPi );
    double costh = EvtRandom::Flat( -1.0, 1.0 );
    double sinth = sqrt( 1.0 - costh * costh );
 
    // b-quark four-momentum in B meson rest frame
 
    EvtVector4R p4b( sqrt( mb * mb + pb * pb ), pb * sinth * sin( phi ),
                     pb * sinth * cos( phi ), pb * costh );
 
    // B meson in its rest frame
    //
    //    EvtVector4R p4B(mB, 0.0, 0.0, 0.0);
    //
    // boost B meson to b-quark rest frame
    //
    //    p4B = boostTo(p4B, p4b);
    //
    //    cout << " B meson mass in b-quark rest frame = " << p4B.mass() << endl;
 
    // boost s, l+ and l- to B meson rest frame
 
    //    EvtVector4R p4s = boostTo(p4sdilep[0], p4B);
    //    p4leptonp       = boostTo(p4ll[0],     p4B);
    //    p4leptonn       = boostTo(p4ll[1],     p4B);
 
    EvtVector4R p4s = boostTo( p4sdilep[0], p4b );
    p4leptonp       = boostTo( p4ll[0], p4b );
    p4leptonn       = boostTo( p4ll[1], p4b );
 
    // spectator quark in B meson rest frame
 
    EvtVector4R p4q( sqrt( pb * pb + _mq * _mq ), -p4b.get( 1 ), -p4b.get( 2 ),
                     -p4b.get( 3 ) );
 
    // hadron system in B meson rest frame
 
    p4xhadron   = p4s + p4q;
    xhadronMass = p4xhadron.mass();
 
    //    cout << "Xs mass = " << xhadronMass << " trial " << im << endl;
  }
 
  // initialize the decay products
 
  xhadron->init( getDaug( 0 ), p4xhadron );
 
  // For B-bar mesons (i.e. containing a b quark) we have the normal
  // order of leptons
  if ( p->getId() == EvtPDL::getId( "anti-B0" ) || p->getId() == EvtPDL::getId( "B-" ) )
  {
    leptonp->init( getDaug( 1 ), p4leptonp );
    leptonn->init( getDaug( 2 ), p4leptonn );
  }
  // For B mesons (i.e. containing a b-bar quark) we need to flip the
  // role of the positive and negative leptons in order to produce the
  // correct forward-backward asymmetry between the two leptons
  else
  {
    leptonp->init( getDaug( 1 ), p4leptonn );
    leptonn->init( getDaug( 2 ), p4leptonp );
  }
 
  return;
}