TEvtGen/EvtGenModels/EvtBTo3piCP.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) 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   // perform common blocks initialization before
  81   // first use
  82   double alpha=getArg(1);
  83   int iset;
  84
  85   iset=10000;
  86
  87   double p4piplus[4],p4piminus[4],p4gamm1[4],p4gamm2[4];
  88
  89   double realA,imgA,realbarA,imgbarA;
  90
  91 #ifdef WIN32
  92   EVT3PIONS(&alpha,&iset,p4piplus,p4piminus,p4gamm1,p4gamm2,
  93              &realA,&imgA,&realbarA,&imgbarA);
  94 #else
  95   evt3pions_(&alpha,&iset,p4piplus,p4piminus,p4gamm1,p4gamm2,
  96              &realA,&imgA,&realbarA,&imgbarA);
  97 #endif
  98
  99   setProbMax(1.5);
 100
 101 }
 102
 103 void EvtBTo3piCP::decay( EvtParticle *p){
 104
 105   //added by Lange Jan4,2000
 106   static EvtId B0=EvtPDL::getId("B0");
 107   static EvtId B0B=EvtPDL::getId("anti-B0");
 108
 109   double t;
 110   EvtId other_b;
 111
 112   EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.5);
 113
 114   EvtParticle *pip,*pim,*pi0;
 115
 116   p->makeDaughters(getNDaug(),getDaugs());
 117
 118   //  p->init_daug(SCALAR,&pip,SCALAR,&pim,SCALAR,&pi0);
 119   pip=p->getDaug(0);
 120   pim=p->getDaug(1);
 121   pi0=p->getDaug(2);
 122
 123   EvtVector4R p4[3];
 124
 125   double dm=getArg(0);
 126   double alpha=getArg(1);
 127   int iset;
 128
 129   iset=0;
 130
 131   double p4piplus[4],p4piminus[4],p4gamm1[4],p4gamm2[4];
 132
 133   double realA,imgA,realbarA,imgbarA;
 134
 135 #ifdef WIN32
 136   EVT3PIONS(&alpha,&iset,p4piplus,p4piminus,p4gamm1,p4gamm2,
 137             &realA,&imgA,&realbarA,&imgbarA);
 138 #else
 139   evt3pions_(&alpha,&iset,p4piplus,p4piminus,p4gamm1,p4gamm2,
 140              &realA,&imgA,&realbarA,&imgbarA);
 141 #endif
 142
 143   p4[0].set(p4piplus[3],p4piplus[0],p4piplus[1],p4piplus[2]);
 144   p4[1].set(p4piminus[3],p4piminus[0],p4piminus[1],p4piminus[2]);
 145   p4[2].set(p4gamm1[3]+p4gamm2[3],p4gamm1[0]+p4gamm2[0],
 146             p4gamm1[1]+p4gamm2[1],p4gamm1[2]+p4gamm2[2]);
 147
 148   if (pip->getId()==EvtPDL::getId("pi+")) {
 149     pip->init( getDaug(0), p4[0] );
 150     pim->init( getDaug(1), p4[1] );
 151   }
 152   else {
 153     pip->init( getDaug(0), p4[1] );
 154     pim->init( getDaug(1), p4[0] );
 155   }
 156
 157    pi0->init( getDaug(2), p4[2] );
 158
 159    EvtComplex amp;
 160
 161    EvtComplex A(realA,imgA);
 162    EvtComplex Abar(realbarA,imgbarA);
 163
 164    if (other_b==B0B){
 165      amp=A*cos(dm*t/(2*EvtConst::c))+
 166        EvtComplex(0.,1.)*Abar*sin(dm*t/(2*EvtConst::c));
 167    }
 168    if (other_b==B0){
 169      amp=Abar*cos(dm*t/(2*EvtConst::c))+
 170        EvtComplex(0.,1.)*A*sin(dm*t/(2*EvtConst::c));
 171    }
 172
 173    vertex(amp);
 174
 175   return ;
 176 }
 177