TEvtGen/EvtGenModels/EvtBHadronic.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: EvtBHadronic.cc
  12 //
  13 // Description: Model for hadronic B decays. Uses naive factorization.
  14 //
  15 // Modification history:
  16 //
  17 //    RYD     June 14, 1997         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include "EvtGenBase/EvtParticle.hh"
  24 #include "EvtGenModels/EvtISGW2FF.hh"
  25 #include "EvtGenBase/EvtGenKine.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenModels/EvtISGW2FF.hh"
  28 #include "EvtGenBase/EvtTensor4C.hh"
  29 #include "EvtGenBase/EvtVector4C.hh"
  30 #include "EvtGenBase/EvtVector4R.hh"
  31 #include "EvtGenModels/EvtBHadronic.hh"
  32 #include "EvtGenBase/EvtReport.hh"
  33 #include <string>
  34 using std::endl;
  35
  36 EvtBHadronic::~EvtBHadronic() {}
  37
  38 std::string EvtBHadronic::getName(){
  39
  40   return "BHADRONIC";
  41
  42 }
  43
  44
  45 EvtDecayBase* EvtBHadronic::clone(){
  46
  47   return new EvtBHadronic;
  48
  49 }
  50
  51
  52 void EvtBHadronic::init(){
  53
  54   // check that there are 2 argument
  55   checkNArg(2);
  56
  57
  58 }
  59
  60 void EvtBHadronic::decay( EvtParticle *p){
  61
  62   //added by Lange Jan4,2000
  63   static EvtId B0=EvtPDL::getId("B0");
  64   static EvtId D0=EvtPDL::getId("D0");
  65   static EvtId DST0=EvtPDL::getId("D*0");
  66   static EvtId D3P10=EvtPDL::getId("D'_10");
  67   static EvtId D3P20=EvtPDL::getId("D_2*0");
  68   static EvtId D3P00=EvtPDL::getId("D_0*0");
  69   static EvtId D1P10=EvtPDL::getId("D_10");
  70
  71   static EvtISGW2FF ffmodel;
  72
  73   p->initializePhaseSpace(getNDaug(),getDaugs());
  74
  75   EvtVector4R p4[MAX_DAUG];
  76   double m;
  77
  78   m = p->mass();
  79
  80
  81   int i,j;
  82
  83   for ( i=0; i<getNDaug(); i++) {
  84     p4[i]=p->getDaug(i)->getP4();
  85   }
  86
  87   int bcurrent,wcurrent;
  88   int nbcurrent=0;
  89   int nwcurrent=0;
  90
  91   bcurrent=(int)getArg(0);
  92   wcurrent=(int)getArg(1);
  93
  94   EvtVector4C jb[5],jw[5];
  95   EvtTensor4C g,tds;
  96
  97   EvtVector4R p4b;
  98   p4b.set(p->mass(),0.0,0.0,0.0);
  99
 100   EvtVector4R q;
 101   double q2;
 102
 103   EvtComplex ep_meson_bb[5];
 104   double f,gf,ap,am;
 105   double fp,fm;
 106   double hf,kf,bp,bm;
 107   EvtVector4R pp,pm;
 108   EvtVector4C ep_meson_b[5];
 109
 110   switch (bcurrent) {
 111
 112   // D0
 113   case 1:
 114     q=p4b-p4[0];
 115     q2=q*q;
 116     nbcurrent=1;
 117     ffmodel.getscalarff(B0,D0,EvtPDL::getMeanMass(D0),
 118                         EvtPDL::getMeanMass(getDaugs()[1]),&fp,&fm);
 119     jb[0]=EvtVector4C(fp*(p4b+p4[0])-fm*q);
 120     break;
 121     // D*
 122   case 2:
 123     q=p4b-p4[0];
 124     q2=q*q;
 125     nbcurrent=3;
 126     ffmodel.getvectorff(B0,DST0,EvtPDL::getMeanMass(DST0),q2,&f,&gf,&ap,&am);
 127
 128     g.setdiag(1.0,-1.0,-1.0,-1.0);
 129     tds = -f*g
 130       -ap*(directProd(p4b,p4b)+directProd(p4b,p4[0]))
 131       -gf*EvtComplex(0.0,1.0)*dual(directProd(p4[0]+p4b,p4b-p4[0]))
 132       -am*((directProd(p4b,p4b)-directProd(p4b,p4[0])));
 133     jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
 134     jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
 135     jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
 136     break;
 137     // D1
 138   case 3:
 139     q=p4b-p4[0];
 140     q2=q*q;
 141     nbcurrent=3;
 142     ffmodel.getvectorff(B0,D3P10,EvtPDL::getMeanMass(D3P10),q2,&f,&gf,&ap,&am);
 143
 144     g.setdiag(1.0,-1.0,-1.0,-1.0);
 145     tds = -f*g
 146       -ap*(directProd(p4b,p4b)+directProd(p4b,p4[0]))
 147       -gf*EvtComplex(0.0,1.0)*dual(directProd(p4[0]+p4b,p4b-p4[0]))
 148       -am*((directProd(p4b,p4b)-directProd(p4b,p4[0])));
 149     jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
 150     jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
 151     jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
 152
 153     break;
 154     // D2*
 155   case 4:
 156     q=p4b-p4[0];
 157     q2=q*q;
 158     nbcurrent=5;
 159     ffmodel.gettensorff(B0,D3P20,EvtPDL::getMeanMass(D3P20),q2,&hf,&kf,&bp,&bm);
 160     g.setdiag(1.0,-1.0,-1.0,-1.0);
 161
 162     ep_meson_b[0] = ((p->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
 163     ep_meson_b[1] = ((p->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
 164     ep_meson_b[2] = ((p->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
 165     ep_meson_b[3] = ((p->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
 166     ep_meson_b[4] = ((p->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();
 167
 168     pp=p4b+p4[0];
 169     pm=p4b-p4[0];
 170
 171     ep_meson_bb[0]=ep_meson_b[0]*(p4b);
 172     ep_meson_bb[1]=ep_meson_b[1]*(p4b);
 173     ep_meson_bb[2]=ep_meson_b[2]*(p4b);
 174     ep_meson_bb[3]=ep_meson_b[3]*(p4b);
 175     ep_meson_bb[4]=ep_meson_b[4]*(p4b);
 176
 177     jb[0]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[0])
 178        -kf*ep_meson_b[0]
 179        -bp*ep_meson_bb[0]*pp-bm*ep_meson_bb[0]*pm;
 180
 181     jb[1]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[1])
 182        -kf*ep_meson_b[1]
 183        -bp*ep_meson_bb[1]*pp-bm*ep_meson_bb[1]*pm;
 184
 185     jb[2]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[2])
 186        -kf*ep_meson_b[2]
 187        -bp*ep_meson_bb[2]*pp-bm*ep_meson_bb[2]*pm;
 188
 189     jb[3]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[3])
 190        -kf*ep_meson_b[3]
 191        -bp*ep_meson_bb[3]*pp-bm*ep_meson_bb[3]*pm;
 192
 193     jb[4]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[4])
 194        -kf*ep_meson_b[4]
 195        -bp*ep_meson_bb[4]*pp-bm*ep_meson_bb[4]*pm;
 196     break;
 197     // D_0*
 198   case 5:
 199     q=p4b-p4[0];
 200     q2=q*q;
 201     double f,gf,ap,am;
 202     nbcurrent=3;
 203     ffmodel.getvectorff(B0,D1P10,EvtPDL::getMeanMass(D1P10),q2,&f,&gf,&ap,&am);
 204     g.setdiag(1.0,-1.0,-1.0,-1.0);
 205     tds = -f*g
 206       -ap*(directProd(p4b,p4b)+directProd(p4b,p4[0]))
 207       +gf*EvtComplex(0.0,1.0)*dual(directProd(p4[0]+p4b,p4b-p4[0]))
 208       -am*((directProd(p4b,p4b)-directProd(p4b,p4[0])));
 209     jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
 210     jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
 211     jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
 212     break;
 213     // D_1
 214   case 6:
 215     q=p4b-p4[0];
 216     q2=q*q;
 217     nbcurrent=1;
 218     ffmodel.getscalarff(B0,D3P00,EvtPDL::getMeanMass(D3P00),q2,&fp,&fm);
 219     jb[0]=fp*(p4b+p4[0])+fm*q;
 220     break;
 221   default:
 222     report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown hadronic current."<<endl;
 223
 224   }
 225
 226   double norm;
 227
 228   switch (wcurrent) {
 229
 230
 231   case 1: case 3: case 4:
 232     nwcurrent=1;
 233     jw[0]=p4[getNDaug()-1];
 234     break;
 235
 236   case 2: case 5: case 6:
 237     nwcurrent=3;
 238     norm=1.0/sqrt(p4[1].get(0)*p4[1].get(0)/p4[1].mass2()-1.0);
 239     jw[0]=norm*p->getDaug(getNDaug()-1)->epsParent(0);
 240     jw[1]=norm*p->getDaug(getNDaug()-1)->epsParent(1);
 241     jw[2]=norm*p->getDaug(getNDaug()-1)->epsParent(2);
 242     break;
 243
 244
 245   default:
 246     report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown W current."<<endl;
 247
 248   }
 249
 250   if (nbcurrent==1&&nwcurrent==1){
 251     vertex(jb[0]*jw[0]);
 252     return;
 253   }
 254
 255   if (nbcurrent==1){
 256     for(j=0;j<nwcurrent;j++){
 257       vertex(j,jb[0]*jw[j]);
 258     }
 259     return;
 260   }
 261
 262   if (nwcurrent==1){
 263     for(j=0;j<nbcurrent;j++){
 264       vertex(j,jb[j]*jw[0]);
 265     }
 266     return;
 267   }
 268
 269   for(j=0;j<nbcurrent;j++){
 270     for(i=0;i<nwcurrent;i++){
 271       vertex(j,i,jb[j]*jw[i]);
 272     }
 273   }
 274   return;
 275
 276 }
 277
 278
 279
 280
 281
 282
 283