--- /dev/null
+// $Id: EvtIncoherentMixing.cpp,v 1.13 2009-11-27 09:09:41 mwhitehe Exp $
+// Include files
+
+
+// local
+#include "EvtGenBase/EvtIncoherentMixing.hh"
+#include <stdlib.h>
+#include "EvtGenBase/EvtPDL.hh"
+#include "EvtGenBase/EvtId.hh"
+#include "EvtGenBase/EvtRandom.hh"
+
+//-----------------------------------------------------------------------------
+// Implementation file for class : EvtIncoherentMixing
+//
+// 2003-10-09 : Patrick Robbe
+//-----------------------------------------------------------------------------
+
+
+bool EvtIncoherentMixing::_doB0Mixing = false ;
+bool EvtIncoherentMixing::_doBsMixing = false ;
+bool EvtIncoherentMixing::_enableFlip = false ;
+double EvtIncoherentMixing::_dGammad = 0. ;
+double EvtIncoherentMixing::_deltamd = 0.502e12 ;
+// dGamma_s corresponds to DeltaGamma / Gamma = 10 %
+double EvtIncoherentMixing::_dGammas = 6.852e10 ;
+double EvtIncoherentMixing::_deltams = 20.e12 ;
+
+//=============================================================================
+// Standard constructor, initializes variables
+//=============================================================================
+EvtIncoherentMixing::EvtIncoherentMixing( ) {
+ _doB0Mixing = false ;
+ _doBsMixing = false ;
+ _dGammad = 0. ;
+ // dGammas corresponds to DeltaGamma / Gamma = 10 %
+ _dGammas = 6.852e10 ;
+ _deltamd = 0.502e12 ;
+ _deltams = 20.e12 ;
+ _enableFlip = false ;
+}
+//=============================================================================
+EvtIncoherentMixing::~EvtIncoherentMixing( )
+{
+}
+// ============================================================================
+void EvtIncoherentMixing::incoherentB0Mix( const EvtId id, double &t ,
+ int &mix )
+{
+ static EvtId B0 = EvtPDL::getId( "B0" ) ;
+ static EvtId B0B = EvtPDL::getId( "anti-B0" ) ;
+
+ if ( ( B0 != id ) && ( B0B != id ) ) {
+ report(ERROR,"EvtGen") << "Bad configuration in incoherentB0Mix"
+ << std::endl ;
+ ::abort() ;
+ }
+
+ double x = getdeltamd() * EvtPDL::getctau( B0 ) / EvtConst::c ;
+
+ double y = getdGammad() * ( EvtPDL::getctau( B0 ) / EvtConst::c ) / 2. ;
+
+ double fac = 1. ; // No CP violation
+
+ double mixprob = ( x*x + y*y ) / ( x*x + y*y + ( 1./fac ) *
+ ( 2. + x*x - y*y ) ) ;
+
+ int mixsign ;
+
+ // decide if state is mixed
+ mixsign = ( mixprob > EvtRandom::Flat( 0. , 1. ) ) ? -1 : 1 ;
+
+ double prob ;
+
+ do {
+ t = -log( EvtRandom::Flat() ) * EvtPDL::getctau( B0 ) / ( 1. - y ) ;
+ prob = ( 1. + exp( -2. * y * t / EvtPDL::getctau( B0 ) ) +
+ mixsign * 2. * exp( -y * t / EvtPDL::getctau( B0 ) ) *
+ cos( getdeltamd() * t / EvtConst::c ) ) / 2. ;
+ } while ( prob < 2. * EvtRandom::Flat() ) ;
+
+ mix = 0 ;
+ if ( mixsign == -1 ) mix = 1 ;
+
+ return ;
+}
+// ============================================================================
+void EvtIncoherentMixing::incoherentBsMix( const EvtId id, double &t ,
+ int &mix )
+{
+ static EvtId BS = EvtPDL::getId( "B_s0" ) ;
+ static EvtId BSB = EvtPDL::getId( "anti-B_s0" ) ;
+
+ if ( ( BS != id ) && ( BSB != id ) ) {
+ report(ERROR,"EvtGen") << "Bad configuration in incoherentBsMix"
+ << std::endl ;
+ ::abort() ;
+ }
+
+ double x = getdeltams() * EvtPDL::getctau( BS ) / EvtConst::c ;
+
+ double y = getdGammas() * ( EvtPDL::getctau( BS ) / EvtConst::c ) / 2. ;
+
+ double fac = 1. ; // No CP violation
+
+ double mixprob = ( x*x + y*y ) / ( x*x + y*y + ( 1./fac ) *
+ ( 2. + x*x - y*y ) ) ;
+
+ int mixsign ;
+
+ // decide if state is mixed
+ mixsign = ( mixprob > EvtRandom::Flat( 0. , 1. ) ) ? -1 : 1 ;
+
+ double prob ;
+
+ do {
+ t = -log( EvtRandom::Flat() ) * EvtPDL::getctau( BS ) / ( 1. - y ) ;
+ prob = ( 1. + exp( -2. * y * t / EvtPDL::getctau( BS ) ) +
+ mixsign * 2. * exp( -y * t / EvtPDL::getctau( BS ) ) *
+ cos( getdeltams() * t / EvtConst::c ) ) / 2. ;
+ } while ( prob < 2. * EvtRandom::Flat() ) ;
+
+ mix = 0 ;
+ if ( mixsign == -1 ) mix = 1 ;
+
+ return ;
+}
+
+// ========================================================================
+bool EvtIncoherentMixing::isBsMixed ( EvtParticle * p )
+{
+ if ( ! ( p->getParent() ) ) return false ;
+
+ static EvtId BS0=EvtPDL::getId("B_s0");
+ static EvtId BSB=EvtPDL::getId("anti-B_s0");
+
+ if ( ( p->getId() != BS0 ) && ( p->getId() != BSB ) ) return false ;
+
+ if ( ( p->getParent()->getId() == BS0 ) ||
+ ( p->getParent()->getId() == BSB ) ) return true ;
+
+ return false ;
+}
+
+// ========================================================================
+bool EvtIncoherentMixing::isB0Mixed ( EvtParticle * p )
+{
+ if ( ! ( p->getParent() ) ) return false ;
+
+ static EvtId B0 =EvtPDL::getId("B0");
+ static EvtId B0B=EvtPDL::getId("anti-B0");
+
+ if ( ( p->getId() != B0 ) && ( p->getId() != B0B ) ) return false ;
+
+ if ( ( p->getParent()->getId() == B0 ) ||
+ ( p->getParent()->getId() == B0B ) ) return true ;
+
+ return false ;
+}
+//============================================================================
+// Return the tag of the event (ie the anti-flavour of the produced
+// B meson). Flip the flavour of the event with probB probability
+//============================================================================
+void EvtIncoherentMixing::OtherB( EvtParticle * p ,
+ double & t ,
+ EvtId & otherb ,
+ double probB )
+{
+ //if(p->getId() == B0 || p->getId() == B0B)
+ //added by liming Zhang
+ enableFlip();
+ if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
+ p->getParent()->setLifetime() ;
+ t = p->getParent()->getLifetime() ;
+ }
+ else {
+ p->setLifetime() ;
+ t = p->getLifetime() ;
+ }
+
+ if ( flipIsEnabled() ) {
+ //std::cout << " liming << flipIsEnabled " << std::endl;
+ // Flip the flavour of the particle with probability probB
+ bool isFlipped = ( EvtRandom::Flat( 0. , 1. ) < probB ) ;
+
+ if ( isFlipped ) {
+ if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
+ p->getParent()
+ ->setId( EvtPDL::chargeConj( p->getParent()->getId() ) ) ;
+ p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
+ }
+ else {
+ p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
+ }
+ }
+ }
+
+ if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
+ // if B has mixed, tag flavour is charge conjugate of parent of B-meson
+ otherb = EvtPDL::chargeConj( p->getParent()->getId() ) ;
+ }
+ else {
+ // else it is opposite flavour than this B hadron
+ otherb = EvtPDL::chargeConj( p->getId() ) ;
+ }
+
+ return ;
+}
+//============================================================================
+// Return the tag of the event (ie the anti-flavour of the produced
+// B meson). No flip
+//============================================================================
+void EvtIncoherentMixing::OtherB( EvtParticle * p ,
+ double & t ,
+ EvtId & otherb )
+{
+ if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
+ p->getParent()->setLifetime() ;
+ t = p->getParent()->getLifetime() ;
+ }
+ else {
+ p->setLifetime() ;
+ t = p->getLifetime() ;
+ }
+
+ if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
+ // if B has mixed, tag flavour is charge conjugate of parent of B-meson
+ otherb = EvtPDL::chargeConj( p->getParent()->getId() ) ;
+ }
+ else {
+ // else it is opposite flavour than this B hadron
+ otherb = EvtPDL::chargeConj( p->getId() ) ;
+ }
+
+ return ;
+}
+
+
+// activate or desactivate the Bs mixing
+void EvtIncoherentMixing::setB0Mixing() { _doB0Mixing = true ; }
+void EvtIncoherentMixing::unsetB0Mixing() { _doB0Mixing = false ; }
+
+// activate or desactivate the B0 mixing
+void EvtIncoherentMixing::setBsMixing() { _doBsMixing = true ; }
+void EvtIncoherentMixing::unsetBsMixing() { _doBsMixing = false ; }
+
+// is mixing activated ?
+bool EvtIncoherentMixing::doB0Mixing() { return _doB0Mixing ; }
+bool EvtIncoherentMixing::doBsMixing() { return _doBsMixing ; }
+
+// set values for the mixing
+void EvtIncoherentMixing::setdGammad( double value ) { _dGammad = value ; }
+void EvtIncoherentMixing::setdeltamd( double value ) { _deltamd = value ; }
+void EvtIncoherentMixing::setdGammas( double value ) { _dGammas = value ; }
+void EvtIncoherentMixing::setdeltams( double value ) { _deltams = value ; }
+
+// get parameters for mixing
+double EvtIncoherentMixing::getdGammad() { return _dGammad ; }
+double EvtIncoherentMixing::getdeltamd() { return _deltamd ; }
+double EvtIncoherentMixing::getdGammas() { return _dGammas ; }
+double EvtIncoherentMixing::getdeltams() { return _deltams ; }
+
+bool EvtIncoherentMixing::flipIsEnabled() { return _enableFlip ; }
+void EvtIncoherentMixing::enableFlip() { _enableFlip = true ; }
+void EvtIncoherentMixing::disableFlip() { _enableFlip = false ; }