TEvtGen/EvtGenBase/EvtIncoherentMixing.cpp

   1 // $Id: EvtIncoherentMixing.cpp,v 1.13 2009-11-27 09:09:41 mwhitehe Exp $
   2 // Include files
   3
   4
   5 // local
   6 #include "EvtGenBase/EvtIncoherentMixing.hh"
   7 #include <stdlib.h>
   8 #include "EvtGenBase/EvtPDL.hh"
   9 #include "EvtGenBase/EvtId.hh"
  10 #include "EvtGenBase/EvtRandom.hh"
  11
  12 //-----------------------------------------------------------------------------
  13 // Implementation file for class : EvtIncoherentMixing
  14 //
  15 // 2003-10-09 : Patrick Robbe
  16 //-----------------------------------------------------------------------------
  17
  18
  19 bool EvtIncoherentMixing::_doB0Mixing = false ;
  20 bool EvtIncoherentMixing::_doBsMixing = false ;
  21 bool EvtIncoherentMixing::_enableFlip = false ;
  22 double EvtIncoherentMixing::_dGammad = 0. ;
  23 double EvtIncoherentMixing::_deltamd = 0.502e12 ;
  24 // dGamma_s corresponds to DeltaGamma / Gamma = 10 %
  25 double EvtIncoherentMixing::_dGammas = 6.852e10 ;
  26 double EvtIncoherentMixing::_deltams = 20.e12 ;
  27
  28 //=============================================================================
  29 // Standard constructor, initializes variables
  30 //=============================================================================
  31 EvtIncoherentMixing::EvtIncoherentMixing(  ) {
  32   _doB0Mixing = false ;
  33   _doBsMixing = false ;
  34   _dGammad = 0. ;
  35   // dGammas corresponds to DeltaGamma / Gamma = 10 %
  36   _dGammas = 6.852e10 ;
  37   _deltamd = 0.502e12 ;
  38   _deltams = 20.e12 ;
  39   _enableFlip = false ;
  40 }
  41 //=============================================================================
  42 EvtIncoherentMixing::~EvtIncoherentMixing( )
  43 {
  44 }
  45 // ============================================================================
  46 void EvtIncoherentMixing::incoherentB0Mix( const EvtId id, double &t ,
  47                                            int &mix )
  48 {
  49   static EvtId B0  = EvtPDL::getId( "B0" ) ;
  50   static EvtId B0B = EvtPDL::getId( "anti-B0" ) ;
  51
  52   if ( ( B0 != id ) && ( B0B != id ) ) {
  53     report(ERROR,"EvtGen") << "Bad configuration in incoherentB0Mix"
  54                            << std::endl ;
  55     ::abort() ;
  56   }
  57
  58   double x = getdeltamd() * EvtPDL::getctau( B0 ) / EvtConst::c ;
  59
  60   double y = getdGammad() * ( EvtPDL::getctau( B0 ) / EvtConst::c ) / 2. ;
  61
  62   double fac = 1. ; // No CP violation
  63
  64   double mixprob = ( x*x + y*y ) / ( x*x + y*y + ( 1./fac ) *
  65                                      ( 2. + x*x - y*y ) ) ;
  66
  67   int mixsign ;
  68
  69   // decide if state is mixed
  70   mixsign = ( mixprob > EvtRandom::Flat( 0. , 1. ) ) ? -1 : 1 ;
  71
  72   double prob ;
  73
  74   do {
  75     t = -log( EvtRandom::Flat() ) * EvtPDL::getctau( B0 ) / ( 1. - y ) ;
  76     prob = ( 1. + exp( -2. * y * t / EvtPDL::getctau( B0 ) ) +
  77       mixsign * 2. * exp( -y * t / EvtPDL::getctau( B0 ) ) *
  78       cos( getdeltamd() * t / EvtConst::c ) ) / 2. ;
  79   } while ( prob < 2. * EvtRandom::Flat() ) ;
  80
  81   mix = 0 ;
  82   if ( mixsign == -1 ) mix = 1 ;
  83
  84   return ;
  85 }
  86 // ============================================================================
  87 void EvtIncoherentMixing::incoherentBsMix( const EvtId id, double &t ,
  88                                            int &mix )
  89 {
  90   static EvtId BS  = EvtPDL::getId( "B_s0" ) ;
  91   static EvtId BSB = EvtPDL::getId( "anti-B_s0" ) ;
  92
  93   if ( ( BS != id ) && ( BSB != id ) ) {
  94     report(ERROR,"EvtGen") << "Bad configuration in incoherentBsMix"
  95                            << std::endl ;
  96     ::abort() ;
  97   }
  98
  99   double x = getdeltams() * EvtPDL::getctau( BS ) / EvtConst::c ;
 100
 101   double y = getdGammas() * ( EvtPDL::getctau( BS ) / EvtConst::c ) / 2. ;
 102
 103   double fac = 1. ; // No CP violation
 104
 105   double mixprob = ( x*x + y*y ) / ( x*x + y*y + ( 1./fac ) *
 106                                      ( 2. + x*x - y*y ) ) ;
 107
 108   int mixsign ;
 109
 110   // decide if state is mixed
 111   mixsign = ( mixprob > EvtRandom::Flat( 0. , 1. ) ) ? -1 : 1 ;
 112
 113   double prob ;
 114
 115   do {
 116     t = -log( EvtRandom::Flat() ) * EvtPDL::getctau( BS ) / ( 1. - y ) ;
 117     prob = ( 1. + exp( -2. * y * t / EvtPDL::getctau( BS ) ) +
 118       mixsign * 2. * exp( -y * t / EvtPDL::getctau( BS ) ) *
 119       cos( getdeltams() * t / EvtConst::c ) ) / 2. ;
 120   } while ( prob < 2. * EvtRandom::Flat() ) ;
 121
 122   mix = 0 ;
 123   if ( mixsign == -1 ) mix = 1 ;
 124
 125   return ;
 126 }
 127
 128 // ========================================================================
 129 bool EvtIncoherentMixing::isBsMixed ( EvtParticle * p )
 130 {
 131   if ( ! ( p->getParent() ) ) return false ;
 132
 133   static EvtId BS0=EvtPDL::getId("B_s0");
 134   static EvtId BSB=EvtPDL::getId("anti-B_s0");
 135
 136   if ( ( p->getId() != BS0 ) && ( p->getId() != BSB ) ) return false ;
 137
 138   if ( ( p->getParent()->getId() == BS0 ) ||
 139        ( p->getParent()->getId() == BSB ) ) return true ;
 140
 141   return false ;
 142 }
 143
 144 // ========================================================================
 145 bool EvtIncoherentMixing::isB0Mixed ( EvtParticle * p )
 146 {
 147   if ( ! ( p->getParent() ) ) return false ;
 148
 149   static EvtId B0 =EvtPDL::getId("B0");
 150   static EvtId B0B=EvtPDL::getId("anti-B0");
 151
 152   if ( ( p->getId() != B0 ) && ( p->getId() != B0B ) ) return false ;
 153
 154   if ( ( p->getParent()->getId() == B0 ) ||
 155        ( p->getParent()->getId() == B0B ) ) return true ;
 156
 157   return false ;
 158 }
 159 //============================================================================
 160 // Return the tag of the event (ie the anti-flavour of the produced
 161 // B meson). Flip the flavour of the event with probB probability
 162 //============================================================================
 163 void EvtIncoherentMixing::OtherB( EvtParticle * p ,
 164                                   double & t ,
 165                                   EvtId & otherb ,
 166                                   double probB )
 167 {
 168   //if(p->getId() == B0 || p->getId() == B0B)
 169   //added by liming Zhang
 170   enableFlip();
 171   if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
 172     p->getParent()->setLifetime() ;
 173     t = p->getParent()->getLifetime() ;
 174   }
 175   else {
 176     p->setLifetime() ;
 177     t = p->getLifetime() ;
 178   }
 179
 180   if ( flipIsEnabled() ) {
 181     //std::cout << " liming << flipIsEnabled " << std::endl;
 182     // Flip the flavour of the particle with probability probB
 183     bool isFlipped = ( EvtRandom::Flat( 0. , 1. ) < probB ) ;
 184
 185     if ( isFlipped ) {
 186       if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
 187         p->getParent()
 188           ->setId( EvtPDL::chargeConj( p->getParent()->getId() ) ) ;
 189         p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
 190       }
 191       else {
 192         p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
 193       }
 194     }
 195   }
 196
 197   if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
 198     // if B has mixed, tag flavour is charge conjugate of parent of B-meson
 199     otherb = EvtPDL::chargeConj( p->getParent()->getId() ) ;
 200   }
 201   else {
 202     // else it is opposite flavour than this B hadron
 203     otherb = EvtPDL::chargeConj( p->getId() ) ;
 204   }
 205
 206   return ;
 207 }
 208 //============================================================================
 209 // Return the tag of the event (ie the anti-flavour of the produced
 210 // B meson). No flip
 211 //============================================================================
 212 void EvtIncoherentMixing::OtherB( EvtParticle * p ,
 213                                   double & t ,
 214                                   EvtId & otherb )
 215 {
 216   if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
 217     p->getParent()->setLifetime() ;
 218     t = p->getParent()->getLifetime() ;
 219   }
 220   else {
 221     p->setLifetime() ;
 222     t = p->getLifetime() ;
 223   }
 224
 225   if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
 226     // if B has mixed, tag flavour is charge conjugate of parent of B-meson
 227     otherb = EvtPDL::chargeConj( p->getParent()->getId() ) ;
 228   }
 229   else {
 230     // else it is opposite flavour than this B hadron
 231     otherb = EvtPDL::chargeConj( p->getId() ) ;
 232   }
 233
 234   return ;
 235 }
 236
 237
 238 // activate or desactivate the Bs mixing
 239 void EvtIncoherentMixing::setB0Mixing()   { _doB0Mixing = true ; }
 240 void EvtIncoherentMixing::unsetB0Mixing() { _doB0Mixing = false ; }
 241
 242 // activate or desactivate the B0 mixing
 243 void EvtIncoherentMixing::setBsMixing()   { _doBsMixing = true ; }
 244 void EvtIncoherentMixing::unsetBsMixing() { _doBsMixing = false ; }
 245
 246 // is mixing activated ?
 247 bool EvtIncoherentMixing::doB0Mixing()  { return _doB0Mixing ; }
 248 bool EvtIncoherentMixing::doBsMixing()  { return _doBsMixing ; }
 249
 250 // set values for the mixing
 251 void EvtIncoherentMixing::setdGammad( double value )  { _dGammad = value ; }
 252 void EvtIncoherentMixing::setdeltamd( double value )  { _deltamd = value ; }
 253 void EvtIncoherentMixing::setdGammas( double value )  { _dGammas = value ; }
 254 void EvtIncoherentMixing::setdeltams( double value )  { _deltams = value ; }
 255
 256 // get parameters for mixing
 257 double EvtIncoherentMixing::getdGammad() { return _dGammad ; }
 258 double EvtIncoherentMixing::getdeltamd() { return _deltamd ; }
 259 double EvtIncoherentMixing::getdGammas() { return _dGammas ; }
 260 double EvtIncoherentMixing::getdeltams() { return _deltams ; }
 261
 262 bool EvtIncoherentMixing::flipIsEnabled() { return _enableFlip ; }
 263 void EvtIncoherentMixing::enableFlip() { _enableFlip = true ; }
 264 void EvtIncoherentMixing::disableFlip() { _enableFlip = false ; }