TEvtGen/EvtGenBase/EvtCPUtil.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 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtCPUtil.cc
  12 //
  13 // Description: Utilities needed for generation of CP violating
  14 //              decays.
  15 //
  16 // Modification history:
  17 //
  18 //    RYD     March 24, 1998         Module created
  19 //
  20 //    COWAN   June 10, 2009          Added methods for getting dGamma(s)
  21 //                                   and dm(s) using B(s)0H and B(s)0L.
  22 //------------------------------------------------------------------------
  23 //
  24 #include "EvtGenBase/EvtPatches.hh"
  25 #include "EvtGenBase/EvtPatches.hh"
  26 #include "EvtGenBase/EvtParticle.hh"
  27 #include "EvtGenBase/EvtScalarParticle.hh"
  28 #include "EvtGenBase/EvtRandom.hh"
  29 #include "EvtGenBase/EvtCPUtil.hh"
  30 #include "EvtGenBase/EvtPDL.hh"
  31 #include "EvtGenBase/EvtReport.hh"
  32 #include "EvtGenBase/EvtSymTable.hh"
  33 #include "EvtGenBase/EvtConst.hh"
  34 #include <stdio.h>
  35 #include <stdlib.h>
  36
  37 #include <assert.h>
  38 using std::endl;
  39
  40 EvtCPUtil::EvtCPUtil(int mixingType) {
  41   _enableFlip = false;
  42   _mixingType = mixingType;
  43 }
  44
  45 EvtCPUtil::~EvtCPUtil() {
  46 }
  47
  48 EvtCPUtil* EvtCPUtil::getInstance() {
  49
  50   static EvtCPUtil* theCPUtil = 0;
  51
  52   if (theCPUtil == 0) {
  53     theCPUtil = new EvtCPUtil(1);
  54   }
  55
  56   return theCPUtil;
  57
  58 }
  59
  60 //added two functions for finding the fraction of B0 tags for decays into
  61 //both CP eigenstates and non-CP eigenstates -- NK, Jan. 27th, 1998
  62
  63 void EvtCPUtil::fractB0CP(EvtComplex Af, EvtComplex Abarf,
  64                           double /*deltam*/, double beta, double &fract) {
  65
  66   //This function returns the number of B0 tags for decays into CP-eigenstates
  67   //(the "probB0" in the new EvtOtherB)
  68
  69   //double gamma_B = EvtPDL::getWidth(B0);
  70   //double xd = deltam/gamma_B;
  71   //double xd = 0.65;
  72   double ratio = 1/(1 + 0.65*0.65);
  73
  74   EvtComplex rf, rbarf;
  75
  76   rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
  77   rbarf = EvtComplex(1.0)/rf;
  78
  79   double A2 = real(Af)*real(Af) + imag(Af)*imag(Af);
  80   double Abar2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
  81
  82   double rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
  83   double rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
  84
  85   fract = (Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio))/(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio) + A2*(1+ rf2 + (1 - rf2)*ratio));
  86   return;
  87
  88 }
  89
  90 void EvtCPUtil::fractB0nonCP(EvtComplex Af, EvtComplex Abarf,
  91                              EvtComplex Afbar, EvtComplex Abarfbar,
  92                              double deltam, double beta,
  93                              int flip, double &fract) {
  94
  95   //this function returns the number of B0 tags for decays into non-CP eigenstates
  96   //(the "probB0" in the new EvtOtherB)
  97   //this needs more thought...
  98
  99   //double gamma_B = EvtPDL::getWidth(B0);
 100   //report(INFO,"EvtGen") << "gamma " << gamma_B<< endl;
 101   //double xd = deltam/gamma_B;
 102
 103   //why is the width of B0 0 in PDL??
 104
 105   double xd = 0.65;
 106   double gamma_B = deltam/xd;
 107   double IAf, IAfbar, IAbarf, IAbarfbar;
 108   EvtComplex rf, rfbar, rbarf, rbarfbar;
 109   double rf2, rfbar2, rbarf2, rbarfbar2;
 110   double Af2, Afbar2, Abarf2, Abarfbar2;
 111
 112   rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
 113   rfbar = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarfbar/Afbar;
 114   rbarf = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Af/Abarf;
 115   rbarfbar = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Afbar/Abarfbar;
 116
 117
 118   rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
 119   rfbar2 = real(rfbar)*real(rfbar) + imag(rfbar)*imag(rfbar);
 120   rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
 121   rbarfbar2 = real(rbarfbar)*real(rbarfbar) + imag(rbarfbar)*imag(rbarfbar);
 122
 123   Af2 = real(Af)*real(Af) + imag(Af)*imag(Af);
 124   Afbar2 = real(Afbar)*real(Afbar) + imag(Afbar)*imag(Afbar);
 125   Abarf2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
 126   Abarfbar2 = real(Abarfbar)*real(Abarfbar) + imag(Abarfbar)*imag(Abarfbar);
 127
 128
 129   //
 130   //IAf = integral(gamma(B0->f)), etc.
 131   //
 132
 133   IAf = (Af2/(2*gamma_B))*(1+rf2+(1-rf2)/(1+xd*xd));
 134   IAfbar = (Afbar2/(2*gamma_B))*(1+rfbar2+(1-rfbar2)/(1+xd*xd));
 135   IAbarf = (Abarf2/(2*gamma_B))*(1+rbarf2+(1-rbarf2)/(1+xd*xd));
 136   IAbarfbar = (Abarfbar2/(2*gamma_B))*(1+rbarfbar2+(1-rbarfbar2)/(1+xd*xd));
 137
 138   //flip specifies the relative fraction of fbar events
 139
 140   fract = IAbarf/(IAbarf+IAf) + flip*IAbarfbar/(IAfbar+IAbarfbar);
 141
 142
 143   return;
 144 }
 145
 146 void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb, double probB0){
 147
 148   if (_mixingType == EvtCPUtil::Coherent) {
 149
 150     OtherCoherentB(p, t, otherb, probB0);
 151
 152   } else if (_mixingType == EvtCPUtil::Incoherent) {
 153
 154     OtherIncoherentB(p, t, otherb, probB0);
 155
 156   }
 157
 158 }
 159
 160 void EvtCPUtil::OtherCoherentB( EvtParticle *p,double &t, EvtId &otherb, double probB0){
 161
 162   //Can not call this recursively!!!
 163   static int entryCount=0;
 164   entryCount++;
 165
 166   //added by Lange Jan4,2000
 167   static EvtId B0B=EvtPDL::getId("anti-B0");
 168   static EvtId B0=EvtPDL::getId("B0");
 169   static EvtId BSB=EvtPDL::getId("anti-B_s0");
 170   static EvtId BS=EvtPDL::getId("B_s0");
 171
 172   static EvtId UPS4S=EvtPDL::getId("Upsilon(4S)");
 173
 174   int isB0=EvtRandom::Flat(0.0,1.0)<probB0;
 175
 176   int idaug;
 177
 178   p->setLifetime();
 179
 180   // now get the time between the decay of this B and the other B!
 181
 182   EvtParticle *parent=p->getParent();
 183
 184   EvtParticle *other;
 185
 186   bool incoherentmix=false;
 187
 188   if ((parent!=0)&&(parent->getId()==B0||
 189                     parent->getId()==B0B||
 190                     parent->getId()==BS||
 191                     parent->getId()==BSB)) {
 192     incoherentmix=true;
 193   }
 194
 195   if (incoherentmix) parent=parent->getParent();
 196
 197   if (parent==0||parent->getId()!=UPS4S) {
 198     //Need to make this more general, but for now
 199     //assume no parent. If we have parent of B we
 200     //need to charge conj. full decay tree.
 201
 202
 203     if (parent!=0) {
 204       report(INFO,"EvtGen") << "p="<<EvtPDL::name(p->getId())
 205                             << " parent="<<EvtPDL::name(parent->getId())
 206                             << endl;
 207     }
 208     assert(parent==0);
 209     p->setLifetime();
 210     t=p->getLifetime();
 211     bool needToChargeConj=false;
 212     if (p->getId()==B0B&&isB0) needToChargeConj=true;
 213     if (p->getId()==B0&&!isB0) needToChargeConj=true;
 214     if (p->getId()==BSB&&isB0) needToChargeConj=true;
 215     if (p->getId()==BS&&!isB0) needToChargeConj=true;
 216
 217     if (needToChargeConj) {
 218       p->setId( EvtPDL::chargeConj(p->getId()));
 219       if (incoherentmix) {
 220         p->getDaug(0)->setId(EvtPDL::chargeConj(p->getDaug(0)->getId()));
 221       }
 222     }
 223     otherb=EvtPDL::chargeConj(p->getId());
 224
 225     entryCount--;
 226     return;
 227   }
 228   else{
 229     if (parent->getDaug(0)!=p){
 230       other=parent->getDaug(0);
 231       idaug=0;
 232     }
 233     else{
 234       other=parent->getDaug(1);
 235       idaug=1;
 236     }
 237   }
 238
 239   if (parent != 0 ) {
 240
 241     //if (entryCount>1){
 242     //  report(INFO,"EvtGen") << "Double CP decay:"<<entryCount<<endl;
 243     //}
 244
 245     //kludge!! Lange Mar21, 2003
 246     // if the other B is an alias... don't change the flavor..
 247     if ( other->getId().isAlias() ) {
 248       OtherB(p,t,otherb);
 249       entryCount--;
 250       return;
 251
 252     }
 253
 254     if (entryCount==1){
 255
 256       EvtVector4R p_init=other->getP4();
 257       //int decayed=other->getNDaug()>0;
 258       bool decayed = other->isDecayed();
 259
 260       other->deleteTree();
 261
 262       EvtScalarParticle* scalar_part;
 263
 264       scalar_part=new EvtScalarParticle;
 265       if (isB0) {
 266         scalar_part->init(B0,p_init);
 267       }
 268       else{
 269         scalar_part->init(B0B,p_init);
 270       }
 271       other=(EvtParticle *)scalar_part;
 272       //    other->set_type(EvtSpinType::SCALAR);
 273       other->setDiagonalSpinDensity();
 274
 275       parent->insertDaugPtr(idaug,other);
 276
 277       if (decayed){
 278         //report(INFO,"EvtGen") << "In CP Util calling decay \n";
 279         other->decay();
 280       }
 281
 282     }
 283
 284     otherb=other->getId();
 285
 286     other->setLifetime();
 287     t=p->getLifetime()-other->getLifetime();
 288
 289     otherb = other->getId();
 290
 291   }
 292   else {
 293     report(INFO,"EvtGen") << "We have an error here!!!!"<<endl;
 294     otherb = EvtId(-1,-1);
 295   }
 296
 297   entryCount--;
 298   return ;
 299 }
 300
 301 // ========================================================================
 302 bool EvtCPUtil::isBsMixed ( EvtParticle * p )
 303 {
 304   if ( ! ( p->getParent() ) ) return false ;
 305
 306   static EvtId BS0=EvtPDL::getId("B_s0");
 307   static EvtId BSB=EvtPDL::getId("anti-B_s0");
 308
 309   if ( ( p->getId() != BS0 ) && ( p->getId() != BSB ) ) return false ;
 310
 311   if ( ( p->getParent()->getId() == BS0 ) ||
 312        ( p->getParent()->getId() == BSB ) ) return true ;
 313
 314   return false ;
 315 }
 316
 317 // ========================================================================
 318 bool EvtCPUtil::isB0Mixed ( EvtParticle * p )
 319 {
 320   if ( ! ( p->getParent() ) ) return false ;
 321
 322   static EvtId B0 =EvtPDL::getId("B0");
 323   static EvtId B0B=EvtPDL::getId("anti-B0");
 324
 325   if ( ( p->getId() != B0 ) && ( p->getId() != B0B ) ) return false ;
 326
 327   if ( ( p->getParent()->getId() == B0 ) ||
 328        ( p->getParent()->getId() == B0B ) ) return true ;
 329
 330   return false ;
 331 }
 332 //============================================================================
 333 // Return the tag of the event (ie the anti-flavour of the produced
 334 // B meson). Flip the flavour of the event with probB probability
 335 //============================================================================
 336 void EvtCPUtil::OtherIncoherentB( EvtParticle * p ,
 337                                   double & t ,
 338                                   EvtId & otherb ,
 339                                   double probB )
 340 {
 341   //std::cout<<"New routine running"<<endl;
 342   //if(p->getId() == B0 || p->getId() == B0B)
 343   //added by liming Zhang
 344   enableFlip();
 345   if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
 346     p->getParent()->setLifetime() ;
 347     t = p->getParent()->getLifetime() ;
 348   }
 349   else {
 350     p->setLifetime() ;
 351     t = p->getLifetime() ;
 352   }
 353
 354   if ( flipIsEnabled() ) {
 355     //std::cout << " liming << flipIsEnabled " << std::endl;
 356     // Flip the flavour of the particle with probability probB
 357     bool isFlipped = ( EvtRandom::Flat( 0. , 1. ) < probB ) ;
 358
 359     if ( isFlipped ) {
 360       if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
 361         p->getParent()
 362           ->setId( EvtPDL::chargeConj( p->getParent()->getId() ) ) ;
 363         p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
 364       }
 365       else {
 366         p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
 367       }
 368     }
 369   }
 370
 371   if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
 372     // if B has mixed, tag flavour is charge conjugate of parent of B-meson
 373     otherb = EvtPDL::chargeConj( p->getParent()->getId() ) ;
 374   }
 375   else {
 376     // else it is opposite flavour than this B hadron
 377     otherb = EvtPDL::chargeConj( p->getId() ) ;
 378   }
 379
 380   return ;
 381 }
 382 //============================================================================
 383 void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb){
 384
 385   static EvtId BSB=EvtPDL::getId("anti-B_s0");
 386   static EvtId BS0=EvtPDL::getId("B_s0");
 387   static EvtId B0B=EvtPDL::getId("anti-B0");
 388   static EvtId B0=EvtPDL::getId("B0");
 389   static EvtId D0B=EvtPDL::getId("anti-D0");
 390   static EvtId D0=EvtPDL::getId("D0");
 391   static EvtId UPS4=EvtPDL::getId("Upsilon(4S)");
 392
 393   if (p->getId()==BS0||p->getId()==BSB){
 394     static double ctauL=EvtPDL::getctau(EvtPDL::getId("B_s0L"));
 395     static double ctauH=EvtPDL::getctau(EvtPDL::getId("B_s0H"));
 396     static double ctau=ctauL<ctauH?ctauH:ctauL;
 397     t=-log(EvtRandom::Flat())*ctau;
 398     EvtParticle* parent=p->getParent();
 399     if (parent!=0&&(parent->getId()==BS0||parent->getId()==BSB)){
 400       if (parent->getId()==BS0) otherb=BSB;
 401       if (parent->getId()==BSB) otherb=BS0;
 402       parent->setLifetime(t);
 403       return;
 404     }
 405     if (p->getId()==BS0) otherb=BSB;
 406     if (p->getId()==BSB) otherb=BS0;
 407     p->setLifetime(t);
 408     return;
 409   }
 410
 411   if (p->getId()==D0||p->getId()==D0B){
 412     static double ctauL=EvtPDL::getctau(EvtPDL::getId("D0L"));
 413     static double ctauH=EvtPDL::getctau(EvtPDL::getId("D0H"));
 414     static double ctau=ctauL<ctauH?ctauH:ctauL;
 415     t=-log(EvtRandom::Flat())*ctau;
 416     EvtParticle* parent=p->getParent();
 417     if (parent!=0&&(parent->getId()==D0||parent->getId()==D0B)){
 418       if (parent->getId()==D0) otherb=D0B;
 419       if (parent->getId()==D0B) otherb=D0;
 420       parent->setLifetime(t);
 421       return;
 422     }
 423     if (p->getId()==D0) otherb=D0B;
 424     if (p->getId()==D0B) otherb=D0;
 425     p->setLifetime(t);
 426     return;
 427   }
 428
 429   p->setLifetime();
 430
 431   // now get the time between the decay of this B and the other B!
 432
 433   EvtParticle *parent=p->getParent();
 434
 435   if (parent==0||parent->getId()!=UPS4) {
 436     //report(ERROR,"EvtGen") <<
 437     //  "Warning CP violation with B having no parent!"<<endl;
 438     t=p->getLifetime();
 439     if (p->getId()==B0) otherb=B0B;
 440     if (p->getId()==B0B) otherb=B0;
 441     if (p->getId()==BS0) otherb=BSB;
 442     if (p->getId()==BSB) otherb=BS0;
 443     return;
 444   }
 445   else{
 446     if (parent->getDaug(0)!=p){
 447       otherb=parent->getDaug(0)->getId();
 448       parent->getDaug(0)->setLifetime();
 449       t=p->getLifetime()-parent->getDaug(0)->getLifetime();
 450     }
 451     else{
 452      otherb=parent->getDaug(1)->getId();
 453       parent->getDaug(1)->setLifetime();
 454       t=p->getLifetime()-parent->getDaug(1)->getLifetime();
 455    }
 456   }
 457
 458
 459   return ;
 460 }
 461
 462 // No CP violation is assumed
 463 void EvtCPUtil::incoherentMix(const EvtId id, double &t, int &mix){
 464
 465   int stdHepNum=EvtPDL::getStdHep(id);
 466   stdHepNum=abs(stdHepNum);
 467
 468   EvtId partId=EvtPDL::evtIdFromStdHep(stdHepNum);
 469
 470   std::string partName=EvtPDL::name(partId);
 471   std::string hname=partName+std::string("H");
 472   std::string lname=partName+std::string("L");
 473
 474   EvtId lId=EvtPDL::getId(lname);
 475   EvtId hId=EvtPDL::getId(hname);
 476
 477   double ctauL=EvtPDL::getctau(lId);
 478   double ctauH=EvtPDL::getctau(hId);
 479
 480   // Bug Fixed: Corrected the average as gamma is the relevent parameter
 481   double ctau=2.0*(ctauL*ctauH)/(ctauL+ctauH);
 482   //double ctau=0.5*(ctauL+ctauH);
 483
 484   // Bug Fixed: ctau definition changed above
 485   //double y=(ctauH-ctauL)/(2*ctau);
 486   double y=(ctauH-ctauL)/(ctauH+ctauL);
 487
 488   //deltam and qoverp defined in DECAY.DEC
 489
 490   std::string qoverpParmName=std::string("qoverp_incohMix_")+partName;
 491   std::string mdParmName=std::string("dm_incohMix_")+partName;
 492   int ierr;
 493   double qoverp=atof(EvtSymTable::get(qoverpParmName,ierr).c_str());
 494   double x=atof(EvtSymTable::get(mdParmName,ierr).c_str())*ctau/EvtConst::c;
 495   double fac;
 496
 497   if(id==partId){
 498     fac=1.0/(qoverp*qoverp);
 499   }
 500   else{
 501     fac=qoverp*qoverp;
 502   }
 503
 504   double mixprob=(x*x+y*y)/(x*x+y*y+fac*(2.0+x*x-y*y));
 505
 506   int mixsign;
 507
 508   mixsign=(mixprob>EvtRandom::Flat(0.0,1.0))?-1:1;
 509
 510   double prob;
 511
 512   // Find the longest of the two lifetimes
 513   double ctaulong = ctauL<=ctauH?ctauH:ctauL;
 514
 515   // Bug fixed: Ensure cosine argument is dimensionless so /ctau
 516   do{
 517     t=-log(EvtRandom::Flat())*ctaulong;
 518     prob=1.0+exp(-2.0*fabs(y)*t/ctau)+mixsign*2.0*exp(-fabs(y)*t/ctau)*cos(x*t/ctau);
 519   }while(prob<4.0*EvtRandom::Flat());
 520
 521   mix=0;
 522
 523   if (mixsign==-1) mix=1;
 524
 525   return;
 526 }
 527
 528
 529 double EvtCPUtil::getDeltaGamma(const EvtId id){
 530
 531   int stdHepNum = EvtPDL::getStdHep(id);
 532   stdHepNum = abs(stdHepNum);
 533   EvtId partId = EvtPDL::evtIdFromStdHep(stdHepNum);
 534
 535   std::string partName = EvtPDL::name(partId);
 536   std::string hname = partName + std::string("H");
 537   std::string lname = partName + std::string("L");
 538
 539   EvtId lId = EvtPDL::getId(lname);
 540   EvtId hId = EvtPDL::getId(hname);
 541
 542   double ctauL = EvtPDL::getctau(lId);
 543   double ctauH = EvtPDL::getctau(hId);
 544
 545   double dGamma = (1/ctauL - 1/ctauH)*EvtConst::c;
 546   return dGamma;
 547 }
 548
 549 double EvtCPUtil::getDeltaM(const EvtId id){
 550
 551   int stdHepNum = EvtPDL::getStdHep(id);
 552   stdHepNum = abs(stdHepNum);
 553   EvtId partId = EvtPDL::evtIdFromStdHep(stdHepNum);
 554
 555   std::string partName = EvtPDL::name(partId);
 556   std::string parmName = std::string("dm_incohMix_") + partName;
 557
 558   int ierr;
 559   double dM = atof(EvtSymTable::get(parmName,ierr).c_str());
 560   return dM;
 561 }
 562
 563 bool EvtCPUtil::flipIsEnabled() { return _enableFlip ; }
 564 void EvtCPUtil::enableFlip() { _enableFlip = true ; }
 565 void EvtCPUtil::disableFlip() { _enableFlip = false ; }
 566