bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenModels/EvtSSSCPpng.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtSSSCPpng.cc
//
// Description: Routine to decay B -> 2 scalars taking into account penguin 
//              contributions (assuming single quark dominance for penguins)
//
// Modification history:
//
//    RYD/NK    December 3, 1997          Module created
//
//------------------------------------------------------------------------
// 
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtCPUtil.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtSSSCPpng.hh"
#include "EvtGenBase/EvtId.hh"
#include <string>
#include "EvtGenBase/EvtConst.hh"
 
EvtSSSCPpng::~EvtSSSCPpng() {}
 
void EvtSSSCPpng::getName(std::string& model_name){
 
  model_name="SSS_CP_PNG";     
 
}
 
 
EvtDecayBase* EvtSSSCPpng::clone(){
 
  return new EvtSSSCPpng;
 
}
 
void EvtSSSCPpng::init(){
 
  // check that there are 7 arguments
  checkNArg(7);
  checkNDaug(2);
 
  checkSpinParent(EvtSpinType::SCALAR);
 
  checkSpinDaughter(0,EvtSpinType::SCALAR);
  checkSpinDaughter(1,EvtSpinType::SCALAR);
 
}
 
void EvtSSSCPpng::initProbMax(){
 
   setProbMax(getArg(5)*getArg(5)*(1 + getArg(6)*getArg(6)));
 
}
 
void EvtSSSCPpng::decay( EvtParticle *p ){
 
  //added by Lange Jan4,2000
  static EvtId B0=EvtPDL::getId("B0");
  static EvtId B0B=EvtPDL::getId("anti-B0");
 
  double t;
  EvtId other_b;
 
  p->initializePhaseSpace(getNDaug(),getDaugs());
 
  EvtComplex amp;
 
  EvtComplex A,Abar;
  //EvtComplex ACC, AbarCC;
 
// assume single (top) quark dominance for the penguin.
 
// old: a0=alpha, a1=dm, a2=1, a3=1, a4=0, a5=1, a6=0
// new: a0=beta, a1=gamma, a2=delta, a3=dm, a4=1, a5=1=A_{T}, a6=A_{P}/A_{T}
 
// e.g., for B -> pi pi  
// A_{T} = |V_{ub} V_{ud}| T
// A_{P} = |V_{tb} V_{td}| P
// P and T are purely hadronic matrix elements       
// P/T = 0.055, A_{P}/A_{T} = 0.2 (see Marrocchesi and Paver, hep-ph/9702353)
 
// A = A_{T}( exp(i(beta+gamma)) + (A_{P}/A_{T}) exp(i(delta))
// A_bar = same, except for the sign of the weak phases 
// here, delta = delta_{p}-delta_{t} (rel. strong phase)
 
   A=getArg(5)*(EvtComplex(cos(-getArg(0)-getArg(1)),sin(-getArg(0)-getArg(1)))+getArg(6)*EvtComplex(cos(getArg(2)),sin(getArg(2))));
 
   Abar=getArg(5)*(EvtComplex(cos(getArg(0)+getArg(1)),sin(getArg(0)+getArg(1)))+getArg(6)*EvtComplex(cos(getArg(2)),sin(getArg(2))));
 
// get fraction of B0 tags with these amplitudes
 
  //double xd = 0.65;
  double ratio = 1/(1 + 0.65*0.65);
  
  EvtComplex rf, rbarf;
 
  rf = EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*Abar/A;
  rbarf = EvtComplex(1.0)/rf;
 
  double A2 = real(A)*real(A) + imag(A)*imag(A);
  double Abar2 = real(Abar)*real(Abar) + imag(Abar)*imag(Abar);
  
  double rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);    
  double rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);    
 
  //fraction of B0 _tags_
  double fract =(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio))/(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio) + A2*(1+ rf2 + (1 - rf2)*ratio)); 
  
  EvtCPUtil::OtherB(p,t,other_b,fract);
 
//this method works just as well -- NK
//randomly generate the tag (B0 or B0B) 
 
//  double tag = EvtRandom::Flat(0.0,1.0);
//  if (tag < 0.5) {
//
//   EvtCPUtil::OtherB(p,t,other_b,1.0);
//   other_b = B0;
//  }
//  else {
//   
//   EvtCPUtil::OtherB(p,t,other_b,0.0);
//   other_b = B0B;
//  }
 
//mixing angle = -beta
 
   if (other_b==B0B){
     amp=A*cos(getArg(3)*t/(2*EvtConst::c))+
       EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
       getArg(4)*EvtComplex(0.0,1.0)*Abar*sin(getArg(3)*t/(2*EvtConst::c));
   }
   if (other_b==B0){
     amp=A*EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
       EvtComplex(0.0,1.0)*sin(getArg(3)*t/(2*EvtConst::c))+       
       getArg(4)*Abar*cos(getArg(3)*t/(2*EvtConst::c));
   }
 
   vertex(amp);
 
  return ;
}