TEvtGen/EvtGenModels/EvtBTo3piCP.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: EvtBTo3piCP.cc
  12 //
  13 // Description: Routine to decay B->pi+ pi- pi0
  14 //              and has CP violation.
  15 //
  16 // Modification history:
  17 //
  18 //    RYD/VERSILLE     March 2, 1997        Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include <stdlib.h>
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtCPUtil.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenModels/EvtBTo3piCP.hh"
  29 #include "EvtGenBase/EvtId.hh"
  30 #include <string>
  31 #include "EvtGenBase/EvtConst.hh"
  32
  33 #ifdef WIN32
  34 extern "C" {
  35   extern void __stdcall EVT3PIONS(double *,int *,double *,
  36                                   double *,double *,double *,double *,
  37                                   double *,double *,double *);
  38 }
  39 #else
  40 extern "C" {
  41   extern void evt3pions_(double *,int *,double *,
  42                          double *,double *,double *,double *,
  43                          double *,double *,double *);
  44 }
  45 #endif
  46
  47 EvtBTo3piCP::~EvtBTo3piCP() {}
  48
  49
  50 std::string EvtBTo3piCP::getName(){
  51
  52   return "BTO3PI_CP";
  53
  54 }
  55
  56
  57 EvtDecayBase* EvtBTo3piCP::clone(){
  58
  59   return new EvtBTo3piCP;
  60
  61 }
  62
  63 void EvtBTo3piCP::init(){
  64
  65   // check that there are 2 arguments
  66   checkNArg(2);
  67   checkNDaug(3);
  68
  69   checkSpinParent(EvtSpinType::SCALAR);
  70
  71   checkSpinDaughter(0,EvtSpinType::SCALAR);
  72   checkSpinDaughter(1,EvtSpinType::SCALAR);
  73   checkSpinDaughter(2,EvtSpinType::SCALAR);
  74 }
  75
  76
  77
  78 void EvtBTo3piCP::initProbMax(){
  79
  80   setProbMax(1.5);
  81
  82 }
  83
  84 void EvtBTo3piCP::decay( EvtParticle *p){
  85
  86   //added by Lange Jan4,2000
  87   static EvtId B0=EvtPDL::getId("B0");
  88   static EvtId B0B=EvtPDL::getId("anti-B0");
  89
  90   double t;
  91   EvtId other_b;
  92
  93   EvtCPUtil::OtherB(p,t,other_b);
  94
  95   EvtParticle *pip,*pim,*pi0;
  96
  97   p->makeDaughters(getNDaug(),getDaugs());
  98
  99   //  p->init_daug(SCALAR,&pip,SCALAR,&pim,SCALAR,&pi0);
 100   pip=p->getDaug(0);
 101   pim=p->getDaug(1);
 102   pi0=p->getDaug(2);
 103
 104   EvtVector4R p4[3];
 105
 106   double dm=getArg(0);
 107   double alpha=getArg(1);
 108   int iset;
 109
 110   static int first=1;
 111
 112   if (first==1) {
 113     iset=10000;
 114     first=0;
 115   }
 116   else{
 117     iset=0;
 118   }
 119
 120   double p4piplus[4],p4piminus[4],p4gamm1[4],p4gamm2[4];
 121
 122   double realA,imgA,realbarA,imgbarA;
 123
 124 #ifdef WIN32
 125   EVT3PIONS(&alpha,&iset,p4piplus,p4piminus,p4gamm1,p4gamm2,
 126             &realA,&imgA,&realbarA,&imgbarA);
 127 #else
 128   evt3pions_(&alpha,&iset,p4piplus,p4piminus,p4gamm1,p4gamm2,
 129              &realA,&imgA,&realbarA,&imgbarA);
 130 #endif
 131
 132   p4[0].set(p4piplus[3],p4piplus[0],p4piplus[1],p4piplus[2]);
 133   p4[1].set(p4piminus[3],p4piminus[0],p4piminus[1],p4piminus[2]);
 134   p4[2].set(p4gamm1[3]+p4gamm2[3],p4gamm1[0]+p4gamm2[0],
 135             p4gamm1[1]+p4gamm2[1],p4gamm1[2]+p4gamm2[2]);
 136
 137   if (pip->getId()==EvtPDL::getId("pi+")) {
 138     pip->init( getDaug(0), p4[0] );
 139     pim->init( getDaug(1), p4[1] );
 140   }
 141   else {
 142     pip->init( getDaug(0), p4[1] );
 143     pim->init( getDaug(1), p4[0] );
 144   }
 145
 146    pi0->init( getDaug(2), p4[2] );
 147
 148    EvtComplex amp;
 149
 150    EvtComplex A(realA,imgA);
 151    EvtComplex Abar(realbarA,imgbarA);
 152
 153    if (other_b==B0B){
 154      amp=A*cos(dm*t/(2*EvtConst::c))+
 155        EvtComplex(0.,1.)*Abar*sin(dm*t/(2*EvtConst::c));
 156    }
 157    if (other_b==B0){
 158      amp=Abar*cos(dm*t/(2*EvtConst::c))+
 159        EvtComplex(0.,1.)*A*sin(dm*t/(2*EvtConst::c));
 160    }
 161
 162    vertex(amp);
 163
 164   return ;
 165 }
 166