//
// RYD March 24, 1998 Module created
//
+// COWAN June 10, 2009 Added methods for getting dGamma(s)
+// and dm(s) using B(s)0H and B(s)0L.
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <assert.h>
using std::endl;
+EvtCPUtil::EvtCPUtil(int mixingType) {
+ _enableFlip = false;
+ _mixingType = mixingType;
+}
+
+EvtCPUtil::~EvtCPUtil() {
+}
+
+EvtCPUtil* EvtCPUtil::getInstance() {
+
+ static EvtCPUtil* theCPUtil = 0;
+
+ if (theCPUtil == 0) {
+ theCPUtil = new EvtCPUtil(1);
+ }
+
+ return theCPUtil;
+
+}
//added two functions for finding the fraction of B0 tags for decays into
//both CP eigenstates and non-CP eigenstates -- NK, Jan. 27th, 1998
void EvtCPUtil::fractB0CP(EvtComplex Af, EvtComplex Abarf,
- double deltam, double beta, double &fract) {
-//this function returns the number of B0 tags for decays into CP-eigenstates
-//(the "probB0" in the new EvtOtherB)
+ double /*deltam*/, double beta, double &fract) {
+
+ //This function returns the number of B0 tags for decays into CP-eigenstates
+ //(the "probB0" in the new EvtOtherB)
//double gamma_B = EvtPDL::getWidth(B0);
//double xd = deltam/gamma_B;
return;
}
+
void EvtCPUtil::fractB0nonCP(EvtComplex Af, EvtComplex Abarf,
EvtComplex Afbar, EvtComplex Abarfbar,
double deltam, double beta,
void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb, double probB0){
+ if (_mixingType == EvtCPUtil::Coherent) {
+
+ OtherCoherentB(p, t, otherb, probB0);
+
+ } else if (_mixingType == EvtCPUtil::Incoherent) {
+
+ OtherIncoherentB(p, t, otherb, probB0);
+
+ }
+
+}
+
+void EvtCPUtil::OtherCoherentB( EvtParticle *p,double &t, EvtId &otherb, double probB0){
+
//Can not call this recursively!!!
static int entryCount=0;
entryCount++;
//kludge!! Lange Mar21, 2003
// if the other B is an alias... don't change the flavor..
if ( other->getId().isAlias() ) {
- OtherB(p,t,otherb);
+ OtherB(p,t,otherb);
entryCount--;
return;
return ;
}
+// ========================================================================
+bool EvtCPUtil::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 EvtCPUtil::isB0Mixed ( EvtParticle * p )
+{
+ if ( ! ( p->getParent() ) ) return false ;
+ static EvtId B0 =EvtPDL::getId("B0");
+ static EvtId B0B=EvtPDL::getId("anti-B0");
-void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb){
+ 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 EvtCPUtil::OtherIncoherentB( EvtParticle * p ,
+ double & t ,
+ EvtId & otherb ,
+ double probB )
+{
+ //std::cout<<"New routine running"<<endl;
+ //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 ;
+}
+//============================================================================
+void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb){
static EvtId BSB=EvtPDL::getId("anti-B_s0");
static EvtId BS0=EvtPDL::getId("B_s0");
return;
}
-
-
-
p->setLifetime();
-
// now get the time between the decay of this B and the other B!
EvtParticle *parent=p->getParent();
return ;
}
-
+// No CP violation is assumed
void EvtCPUtil::incoherentMix(const EvtId id, double &t, int &mix){
int stdHepNum=EvtPDL::getStdHep(id);
double ctauL=EvtPDL::getctau(lId);
double ctauH=EvtPDL::getctau(hId);
- double ctau=0.5*(ctauL+ctauH);
- double y=(ctauH-ctauL)/ctau;
- //need to figure out how to get these parameters into the code...
+ // Bug Fixed: Corrected the average as gamma is the relevent parameter
+ double ctau=2.0*(ctauL*ctauH)/(ctauL+ctauH);
+ //double ctau=0.5*(ctauL+ctauH);
+
+ // Bug Fixed: ctau definition changed above
+ //double y=(ctauH-ctauL)/(2*ctau);
+ double y=(ctauH-ctauL)/(ctauH+ctauL);
+
+ //deltam and qoverp defined in DECAY.DEC
std::string qoverpParmName=std::string("qoverp_incohMix_")+partName;
std::string mdParmName=std::string("dm_incohMix_")+partName;
fac=qoverp*qoverp;
}
- double mixprob=(x*x+y*y)/(x*x+y*y+fac*(2+x*x-y*y));
+ double mixprob=(x*x+y*y)/(x*x+y*y+fac*(2.0+x*x-y*y));
int mixsign;
double prob;
+ // Find the longest of the two lifetimes
+ double ctaulong = ctauL<=ctauH?ctauH:ctauL;
+
+ // Bug fixed: Ensure cosine argument is dimensionless so /ctau
do{
- t=-log(EvtRandom::Flat())*ctauL;
- prob=1.0+exp(2.0*y*t/ctau)+mixsign*2.0*exp(y*t/ctau)*cos(x*t);
- }while(prob<4*EvtRandom::Flat());
+ t=-log(EvtRandom::Flat())*ctaulong;
+ prob=1.0+exp(-2.0*fabs(y)*t/ctau)+mixsign*2.0*exp(-fabs(y)*t/ctau)*cos(x*t/ctau);
+ }while(prob<4.0*EvtRandom::Flat());
mix=0;
}
+double EvtCPUtil::getDeltaGamma(const EvtId id){
+ int stdHepNum = EvtPDL::getStdHep(id);
+ stdHepNum = abs(stdHepNum);
+ EvtId partId = EvtPDL::evtIdFromStdHep(stdHepNum);
+ std::string partName = EvtPDL::name(partId);
+ std::string hname = partName + std::string("H");
+ std::string lname = partName + std::string("L");
+
+ EvtId lId = EvtPDL::getId(lname);
+ EvtId hId = EvtPDL::getId(hname);
+ double ctauL = EvtPDL::getctau(lId);
+ double ctauH = EvtPDL::getctau(hId);
+
+ double dGamma = (1/ctauL - 1/ctauH)*EvtConst::c;
+ return dGamma;
+}
+double EvtCPUtil::getDeltaM(const EvtId id){
+ int stdHepNum = EvtPDL::getStdHep(id);
+ stdHepNum = abs(stdHepNum);
+ EvtId partId = EvtPDL::evtIdFromStdHep(stdHepNum);
+
+ std::string partName = EvtPDL::name(partId);
+ std::string parmName = std::string("dm_incohMix_") + partName;
+ int ierr;
+ double dM = atof(EvtSymTable::get(parmName,ierr).c_str());
+ return dM;
+}
+bool EvtCPUtil::flipIsEnabled() { return _enableFlip ; }
+void EvtCPUtil::enableFlip() { _enableFlip = true ; }
+void EvtCPUtil::disableFlip() { _enableFlip = false ; }
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff