TEvtGen/EvtGenModels/EvtSSSCPpng.cxx

   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: EvtSSSCPpng.cc
  12 //
  13 // Description: Routine to decay B -> 2 scalars taking into account penguin
  14 //              contributions (assuming single quark dominance for penguins)
  15 //
  16 // Modification history:
  17 //
  18 //    RYD/NK    December 3, 1997          Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include <stdlib.h>
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtRandom.hh"
  26 #include "EvtGenBase/EvtGenKine.hh"
  27 #include "EvtGenBase/EvtCPUtil.hh"
  28 #include "EvtGenBase/EvtPDL.hh"
  29 #include "EvtGenBase/EvtReport.hh"
  30 #include "EvtGenModels/EvtSSSCPpng.hh"
  31 #include "EvtGenBase/EvtId.hh"
  32 #include <string>
  33 #include "EvtGenBase/EvtConst.hh"
  34
  35 EvtSSSCPpng::~EvtSSSCPpng() {}
  36
  37 std::string EvtSSSCPpng::getName(){
  38
  39   return "SSS_CP_PNG";
  40
  41 }
  42
  43
  44 EvtDecayBase* EvtSSSCPpng::clone(){
  45
  46   return new EvtSSSCPpng;
  47
  48 }
  49
  50 void EvtSSSCPpng::init(){
  51
  52   // check that there are 7 arguments
  53   checkNArg(7);
  54   checkNDaug(2);
  55
  56   checkSpinParent(EvtSpinType::SCALAR);
  57
  58   checkSpinDaughter(0,EvtSpinType::SCALAR);
  59   checkSpinDaughter(1,EvtSpinType::SCALAR);
  60
  61 }
  62
  63 void EvtSSSCPpng::initProbMax(){
  64
  65    setProbMax(getArg(5)*getArg(5)*(1 + getArg(6)*getArg(6)));
  66
  67 }
  68
  69 void EvtSSSCPpng::decay( EvtParticle *p ){
  70
  71   //added by Lange Jan4,2000
  72   static EvtId B0=EvtPDL::getId("B0");
  73   static EvtId B0B=EvtPDL::getId("anti-B0");
  74
  75   double t;
  76   EvtId other_b;
  77
  78   p->initializePhaseSpace(getNDaug(),getDaugs());
  79
  80   EvtComplex amp;
  81
  82   EvtComplex A,Abar;
  83   //EvtComplex ACC, AbarCC;
  84
  85 // assume single (top) quark dominance for the penguin.
  86
  87 // old: a0=alpha, a1=dm, a2=1, a3=1, a4=0, a5=1, a6=0
  88 // new: a0=beta, a1=gamma, a2=delta, a3=dm, a4=1, a5=1=A_{T}, a6=A_{P}/A_{T}
  89
  90 // e.g., for B -> pi pi
  91 // A_{T} = |V_{ub} V_{ud}| T
  92 // A_{P} = |V_{tb} V_{td}| P
  93 // P and T are purely hadronic matrix elements
  94 // P/T = 0.055, A_{P}/A_{T} = 0.2 (see Marrocchesi and Paver, hep-ph/9702353)
  95
  96 // A = A_{T}( exp(i(beta+gamma)) + (A_{P}/A_{T}) exp(i(delta))
  97 // A_bar = same, except for the sign of the weak phases
  98 // here, delta = delta_{p}-delta_{t} (rel. strong phase)
  99
 100    A=getArg(5)*(EvtComplex(cos(-getArg(0)-getArg(1)),sin(-getArg(0)-getArg(1)))+getArg(6)*EvtComplex(cos(getArg(2)),sin(getArg(2))));
 101
 102    Abar=getArg(5)*(EvtComplex(cos(getArg(0)+getArg(1)),sin(getArg(0)+getArg(1)))+getArg(6)*EvtComplex(cos(getArg(2)),sin(getArg(2))));
 103
 104 // get fraction of B0 tags with these amplitudes
 105
 106   //double xd = 0.65;
 107   double ratio = 1/(1 + 0.65*0.65);
 108
 109   EvtComplex rf, rbarf;
 110
 111   rf = EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*Abar/A;
 112   rbarf = EvtComplex(1.0)/rf;
 113
 114   double A2 = real(A)*real(A) + imag(A)*imag(A);
 115   double Abar2 = real(Abar)*real(Abar) + imag(Abar)*imag(Abar);
 116
 117   double rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
 118   double rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
 119
 120   //fraction of B0 _tags_
 121   double fract =(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio))/(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio) + A2*(1+ rf2 + (1 - rf2)*ratio));
 122
 123   EvtCPUtil::OtherB(p,t,other_b,fract);
 124
 125 //this method works just as well -- NK
 126 //randomly generate the tag (B0 or B0B)
 127
 128 //  double tag = EvtRandom::Flat(0.0,1.0);
 129 //  if (tag < 0.5) {
 130 //
 131 //   EvtCPUtil::OtherB(p,t,other_b,1.0);
 132 //   other_b = B0;
 133 //  }
 134 //  else {
 135 //
 136 //   EvtCPUtil::OtherB(p,t,other_b,0.0);
 137 //   other_b = B0B;
 138 //  }
 139
 140 //mixing angle = -beta
 141
 142    if (other_b==B0B){
 143      amp=A*cos(getArg(3)*t/(2*EvtConst::c))+
 144        EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
 145        getArg(4)*EvtComplex(0.0,1.0)*Abar*sin(getArg(3)*t/(2*EvtConst::c));
 146    }
 147    if (other_b==B0){
 148      amp=A*EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
 149        EvtComplex(0.0,1.0)*sin(getArg(3)*t/(2*EvtConst::c))+
 150        getArg(4)*Abar*cos(getArg(3)*t/(2*EvtConst::c));
 151    }
 152
 153    vertex(amp);
 154
 155   return ;
 156 }
 157
 158