TEvtGen/EvtGenModels/EvtSSD_DirectCP.cpp

   1 // $Id: EvtSSD_DirectCP.cpp,v 1.2 2009-03-16 16:24:05 robbep Exp $
   2 // Generation of direct CP violation in hadronic environment
   3 // Patrick Robbe, LHCb,  08 Nov 2006
   4 //
   5 #include <stdlib.h>
   6 #include "EvtGenBase/EvtParticle.hh"
   7 #include "EvtGenBase/EvtRandom.hh"
   8 #include "EvtGenBase/EvtGenKine.hh"
   9 #include "EvtGenBase/EvtCPUtil.hh"
  10 #include "EvtGenBase/EvtPDL.hh"
  11 #include "EvtGenBase/EvtReport.hh"
  12 #include "EvtGenBase/EvtVector4C.hh"
  13 #include "EvtGenBase/EvtTensor4C.hh"
  14 #include "EvtGenModels/EvtSSD_DirectCP.hh"
  15 #include <string>
  16 #include "EvtGenBase/EvtConst.hh"
  17
  18 EvtSSD_DirectCP::~EvtSSD_DirectCP() {}
  19
  20 std::string EvtSSD_DirectCP::getName( ){
  21
  22   return "SSD_DirectCP" ;
  23
  24 }
  25
  26
  27 EvtDecayBase* EvtSSD_DirectCP::clone(){
  28
  29   return new EvtSSD_DirectCP;
  30
  31 }
  32
  33 void EvtSSD_DirectCP::init(){
  34
  35   // check that there is 1 argument and 2-body decay
  36
  37   checkNArg(1);
  38   checkNDaug(2);
  39
  40   EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
  41   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
  42
  43   if ( (!(d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR))||
  44        (!((d2type==EvtSpinType::SCALAR)||(d2type==EvtSpinType::VECTOR)||
  45           (d2type==EvtSpinType::TENSOR)))||
  46        (!((d1type==EvtSpinType::SCALAR)||(d1type==EvtSpinType::VECTOR)||
  47           (d1type==EvtSpinType::TENSOR)))
  48        ) {
  49     report(ERROR,"EvtGen") << "EvtSSD_DirectCP generator expected "
  50                            << "one of the daugters to be a scalar, "
  51                            << "the other either scalar, vector, or tensor, "
  52                            << "found:"
  53                            << EvtPDL::name(getDaug(0)).c_str()
  54                            <<" and "
  55                            <<EvtPDL::name(getDaug(1)).c_str()<<std::endl;
  56     report(ERROR,"EvtGen") << "Will terminate execution!"<<std::endl;
  57     ::abort();
  58   }
  59
  60   _acp = getArg( 0 ) ; // A_CP defined as A_CP = (BR(fbar)-BR(f))/(BR(fbar)+BR(f))
  61
  62 }
  63
  64 void EvtSSD_DirectCP::initProbMax() {
  65   double theProbMax = 1. ;
  66
  67   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
  68   EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
  69   if (d1type==EvtSpinType::TENSOR||d2type==EvtSpinType::TENSOR) theProbMax*=10;
  70
  71   setProbMax(theProbMax);
  72 }
  73
  74 void EvtSSD_DirectCP::decay( EvtParticle *p) {
  75
  76   bool flip = false ;
  77   EvtId daugs[2];
  78
  79   // decide it is B or Bbar:
  80   if ( EvtRandom::Flat(0.,1.) < ( ( 1. - _acp ) / 2. ) ) {
  81     // it is a B
  82     if ( EvtPDL::getStdHep( getParentId() ) < 0 ) flip = true ;
  83   } else {
  84     // it is a Bbar
  85     if ( EvtPDL::getStdHep( getParentId() ) > 0 ) flip = true ;
  86   }
  87
  88   if ( flip ) {
  89     if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
  90       p->getParent()
  91         ->setId( EvtPDL::chargeConj( p->getParent()->getId() ) ) ;
  92       p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
  93     }
  94     else {
  95       p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
  96     }
  97   }
  98
  99   if (!flip) {
 100     daugs[0]=getDaug(0);
 101     daugs[1]=getDaug(1);
 102   }
 103   else{
 104     daugs[0]=EvtPDL::chargeConj(getDaug(0));
 105     daugs[1]=EvtPDL::chargeConj(getDaug(1));
 106   }
 107
 108   EvtParticle *d;
 109   p->initializePhaseSpace(2, daugs);
 110
 111   EvtVector4R p4_parent=p->getP4Restframe();
 112   double m_parent=p4_parent.mass();
 113
 114   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
 115
 116   EvtVector4R momv;
 117   EvtVector4R moms;
 118
 119   if (d2type==EvtSpinType::SCALAR){
 120     d2type=EvtPDL::getSpinType(getDaug(0));
 121     d= p->getDaug(0);
 122     momv = d->getP4();
 123     moms = p->getDaug(1)->getP4();
 124   }
 125   else{
 126     d= p->getDaug(1);
 127     momv = d->getP4();
 128     moms = p->getDaug(0)->getP4();
 129   }
 130
 131   if (d2type==EvtSpinType::SCALAR) {
 132     vertex(1.);
 133   }
 134
 135   if (d2type==EvtSpinType::VECTOR) {
 136
 137     double norm=momv.mass()/(momv.d3mag()*p->mass());
 138
 139     vertex(0,norm*p4_parent*(d->epsParent(0)));
 140     vertex(1,norm*p4_parent*(d->epsParent(1)));
 141     vertex(2,norm*p4_parent*(d->epsParent(2)));
 142
 143   }
 144
 145   if (d2type==EvtSpinType::TENSOR) {
 146
 147     double norm=
 148       d->mass()*d->mass()/(m_parent*d->getP4().d3mag()*d->getP4().d3mag());
 149
 150
 151    vertex(0,norm*d->epsTensorParent(0).cont1(p4_parent)*p4_parent);
 152    vertex(1,norm*d->epsTensorParent(1).cont1(p4_parent)*p4_parent);
 153    vertex(2,norm*d->epsTensorParent(2).cont1(p4_parent)*p4_parent);
 154    vertex(3,norm*d->epsTensorParent(3).cont1(p4_parent)*p4_parent);
 155    vertex(4,norm*d->epsTensorParent(4).cont1(p4_parent)*p4_parent);
 156   }
 157 }
 158
 159 bool EvtSSD_DirectCP::isB0Mixed ( EvtParticle * p ) {
 160   if ( ! ( p->getParent() ) ) return false ;
 161
 162   static EvtId B0 =EvtPDL::getId("B0");
 163   static EvtId B0B=EvtPDL::getId("anti-B0");
 164
 165   if ( ( p->getId() != B0 ) && ( p->getId() != B0B ) ) return false ;
 166
 167   if ( ( p->getParent()->getId() == B0 ) ||
 168        ( p->getParent()->getId() == B0B ) ) return true ;
 169
 170   return false ;
 171 }
 172
 173 bool EvtSSD_DirectCP::isBsMixed ( EvtParticle * p ) {
 174   if ( ! ( p->getParent() ) ) return false ;
 175
 176   static EvtId BS0=EvtPDL::getId("B_s0");
 177   static EvtId BSB=EvtPDL::getId("anti-B_s0");
 178
 179   if ( ( p->getId() != BS0 ) && ( p->getId() != BSB ) ) return false ;
 180
 181   if ( ( p->getParent()->getId() == BS0 ) ||
 182        ( p->getParent()->getId() == BSB ) ) return true ;
 183
 184   return false ;
 185 }
 186
 187 std::string EvtSSD_DirectCP::getParamName(int i) {
 188   switch(i) {
 189   case 0:
 190     return "ACP";
 191   default:
 192     return "";
 193   }
 194 }