TEvtGen/EvtGenModels/EvtSVSCPLH.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) 2001      Caltech, UCSB
  10 //
  11 // Module: EvtSVSCPLH.cc
  12 //
  13 // Description: The decay of a scalar to a scalar and a vector particle are
  14 //              performed with CP violation and different widths for
  15 //              the cp even and odd states. E.g. B->J/psi K_S.
  16 //
  17 // Modification history:
  18 //
  19 //    Ryd       March 29, 2001       Module created
  20 //
  21 //------------------------------------------------------------------------
  22 //
  23 #include "EvtGenBase/EvtPatches.hh"
  24 #include <stdlib.h>
  25 #include "EvtGenBase/EvtParticle.hh"
  26 #include "EvtGenBase/EvtGenKine.hh"
  27 #include "EvtGenBase/EvtCPUtil.hh"
  28 #include "EvtGenBase/EvtPDL.hh"
  29 #include "EvtGenBase/EvtReport.hh"
  30 #include "EvtGenModels/EvtSVSCPLH.hh"
  31 #include "EvtGenBase/EvtId.hh"
  32 #include <string>
  33 #include "EvtGenBase/EvtConst.hh"
  34 #include "EvtGenBase/EvtVector4C.hh"
  35 using std::endl;
  36
  37 EvtSVSCPLH::~EvtSVSCPLH() {}
  38
  39 std::string EvtSVSCPLH::getName(){
  40
  41   return "SVS_CPLH";
  42
  43 }
  44
  45
  46 EvtDecayBase* EvtSVSCPLH::clone(){
  47
  48   return new EvtSVSCPLH;
  49
  50 }
  51
  52 void EvtSVSCPLH::init(){
  53
  54   // check that there are 8 arguments
  55   checkNArg(8);
  56   checkNDaug(2);
  57
  58   checkSpinParent(EvtSpinType::SCALAR);
  59
  60   checkSpinDaughter(0,EvtSpinType::VECTOR);
  61   checkSpinDaughter(1,EvtSpinType::SCALAR);
  62
  63   static double ctau=EvtPDL::getctau(EvtPDL::getId("B0"));
  64
  65   // hbar/s
  66   _dm=getArg(0);
  67   _dgamma=EvtConst::c*getArg(1)/ctau;
  68
  69   _qop=getArg(2)*exp(EvtComplex(0.0,getArg(3)));
  70
  71   _poq=1.0/_qop;
  72
  73   _Af=getArg(4)*exp(EvtComplex(0.0,getArg(5)));
  74   _Abarf=getArg(6)*exp(EvtComplex(0.0,getArg(7)));
  75
  76   if (verbose()){
  77     report(INFO,"EvtGen")<<":EvtSVSCPLH:dm="<<_dm<<endl;
  78     report(INFO,"EvtGen")<<":EvtSVSCPLH:dGamma="<<_dgamma<<endl;
  79     report(INFO,"EvtGen")<<":EvtSVSCPLH:q/p="<<_qop<<endl;
  80     report(INFO,"EvtGen")<<":EvtSVSCPLH:Af="<<_Af<<endl;
  81     report(INFO,"EvtGen")<<":EvtSVSCPLH:Abarf="<<_Abarf<<endl;
  82   }
  83
  84
  85 }
  86
  87 void EvtSVSCPLH::initProbMax(){
  88
  89   //This is probably not quite right, but it should do as a start...
  90   //Anders
  91
  92   setProbMax(4.0*(getArg(4)*getArg(4)+getArg(6)*getArg(6)));
  93
  94 }
  95
  96 void EvtSVSCPLH::decay( EvtParticle *p){
  97
  98   p->initializePhaseSpace(getNDaug(),getDaugs());
  99
 100   static EvtId B0=EvtPDL::getId("B0");
 101   static EvtId B0B=EvtPDL::getId("anti-B0");
 102
 103   double t;
 104   EvtId other_b;
 105
 106   EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.5);
 107
 108   //convert time from mm to seconds
 109   t/=EvtConst::c;
 110
 111   //sign convention is dm=Mheavy-Mlight
 112   //                   dGamma=Gammalight-Gammaheavy
 113   //such that in the standard model both of these are positive.
 114   EvtComplex gp=0.5*(exp(EvtComplex(0.25*t*_dgamma,-0.5*t*_dm))+exp(EvtComplex(-0.25*t*_dgamma,0.5*t*_dm)));
 115   EvtComplex gm=0.5*(exp(EvtComplex(0.25*t*_dgamma,-0.5*t*_dm))-exp(EvtComplex(-0.25*t*_dgamma,0.5*t*_dm)));
 116
 117   EvtComplex amp;
 118
 119   if (other_b==B0B){
 120     amp=gp*_Af+_qop*gm*_Abarf;
 121   }
 122   else if (other_b==B0){
 123     amp=gp*_Abarf+_poq*gm*_Af;
 124   }
 125   else{
 126     report(ERROR,"EvtGen") << "other_b was not B0 or B0B!"<<endl;
 127     ::abort();
 128   }
 129
 130   EvtVector4R p4_parent=p->getP4Restframe();;
 131
 132   double norm=p->getDaug(0)->mass()/(p->getDaug(0)->getP4().d3mag()*p4_parent.mass());
 133
 134   EvtParticle* v=p->getDaug(0);
 135
 136   vertex(0,amp*norm*(p4_parent*(v->epsParent(0))));
 137   vertex(1,amp*norm*(p4_parent*(v->epsParent(1))));
 138   vertex(2,amp*norm*(p4_parent*(v->epsParent(2))));
 139
 140
 141   return ;
 142 }
 143
 144
 145