TEvtGen/EvtGenModels/EvtGoityRoberts.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: EvtGoityRoberts.cc
  12 //
  13 // Description: Routine to decay vector-> scalar scalar
  14 //
  15 // Modification history:
  16 //
  17 //    RYD     November 24, 1996        Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include "EvtGenBase/EvtParticle.hh"
  24 #include "EvtGenBase/EvtGenKine.hh"
  25 #include "EvtGenBase/EvtPDL.hh"
  26 #include "EvtGenBase/EvtReport.hh"
  27 #include "EvtGenModels/EvtGoityRoberts.hh"
  28 #include "EvtGenBase/EvtTensor4C.hh"
  29 #include "EvtGenBase/EvtDiracSpinor.hh"
  30 #include <string>
  31 #include "EvtGenBase/EvtVector4C.hh"
  32
  33 EvtGoityRoberts::~EvtGoityRoberts() {}
  34
  35 std::string EvtGoityRoberts::getName(){
  36
  37   return "GOITY_ROBERTS";
  38
  39 }
  40
  41
  42 EvtDecayBase* EvtGoityRoberts::clone(){
  43
  44   return new EvtGoityRoberts;
  45
  46 }
  47
  48 void EvtGoityRoberts::init(){
  49
  50   // check that there are 0 arguments
  51   checkNArg(0);
  52   checkNDaug(4);
  53
  54   checkSpinParent(EvtSpinType::SCALAR);
  55   checkSpinDaughter(1,EvtSpinType::SCALAR);
  56   checkSpinDaughter(2,EvtSpinType::DIRAC);
  57   checkSpinDaughter(3,EvtSpinType::NEUTRINO);
  58
  59 }
  60
  61
  62 void EvtGoityRoberts::initProbMax() {
  63
  64    setProbMax( 3000.0);
  65 }
  66
  67 void EvtGoityRoberts::decay( EvtParticle *p){
  68
  69   //added by Lange Jan4,2000
  70   static EvtId DST0=EvtPDL::getId("D*0");
  71   static EvtId DSTB=EvtPDL::getId("anti-D*0");
  72   static EvtId DSTP=EvtPDL::getId("D*+");
  73   static EvtId DSTM=EvtPDL::getId("D*-");
  74   static EvtId D0=EvtPDL::getId("D0");
  75   static EvtId D0B=EvtPDL::getId("anti-D0");
  76   static EvtId DP=EvtPDL::getId("D+");
  77   static EvtId DM=EvtPDL::getId("D-");
  78
  79
  80
  81   EvtId meson=getDaug(0);
  82
  83   if (meson==DST0||meson==DSTP||meson==DSTM||meson==DSTB) {
  84     DecayBDstarpilnuGR(p,getDaug(0),getDaug(2),getDaug(3));
  85   }
  86   else{
  87     if (meson==D0||meson==DP||meson==DM||meson==D0B) {
  88       DecayBDpilnuGR(p,getDaug(0),getDaug(2),getDaug(3));
  89     }
  90     else{
  91       report(ERROR,"EvtGen") << "Wrong daugther in EvtGoityRoberts!\n";
  92     }
  93   }
  94   return ;
  95 }
  96
  97 void EvtGoityRoberts::DecayBDstarpilnuGR(EvtParticle *pb,EvtId ndstar,
  98                                          EvtId nlep, EvtId /*nnu*/)
  99 {
 100
 101   pb->initializePhaseSpace(getNDaug(),getDaugs());
 102
 103   //added by Lange Jan4,2000
 104   static EvtId EM=EvtPDL::getId("e-");
 105   static EvtId EP=EvtPDL::getId("e+");
 106   static EvtId MUM=EvtPDL::getId("mu-");
 107   static EvtId MUP=EvtPDL::getId("mu+");
 108
 109   EvtParticle *dstar, *pion, *lepton, *neutrino;
 110
 111   // pb->makeDaughters(getNDaug(),getDaugs());
 112   dstar=pb->getDaug(0);
 113   pion=pb->getDaug(1);
 114   lepton=pb->getDaug(2);
 115   neutrino=pb->getDaug(3);
 116
 117   EvtVector4C l1, l2, et0, et1, et2;
 118
 119   EvtVector4R v,vp,p4_pi;
 120   double w;
 121
 122   v.set(1.0,0.0,0.0,0.0);       //4-velocity of B meson
 123   vp=(1.0/dstar->getP4().mass())*dstar->getP4();  //4-velocity of D*
 124   p4_pi=pion->getP4();
 125
 126   w=v*vp;                       //four velocity transfere.
 127
 128   EvtTensor4C omega;
 129
 130   double mb=EvtPDL::getMeanMass(pb->getId());     //B mass
 131   double md=EvtPDL::getMeanMass(ndstar);   //D* mass
 132
 133   EvtComplex dmb(0.0460,-0.5*0.00001);   // B*-B mass splitting ?
 134   EvtComplex dmd(0.1421,-0.5*0.00006);
 135                              // The last two sets of numbers should
 136                              // be correctly calculated from the
 137                              // dstar and pion charges.
 138   double g = 0.5;         // EvtAmplitude proportional to these coupling constants
 139   double alpha3 =  0.690; // See table I in G&R's paper
 140   double alpha1 = -1.430;
 141   double alpha2 = -0.140;
 142   double f0 = 0.093;      // The pion decay constants set to 93 MeV
 143
 144   EvtComplex dmt3(0.563,-0.5*0.191); // Mass splitting = dmt - iGamma/2
 145   EvtComplex dmt1(0.392,-0.5*1.040);
 146   EvtComplex dmt2(0.709,-0.5*0.405);
 147
 148   double betas=0.285;      // magic number for meson wave function ground state
 149   double betap=0.280;      // magic number for meson wave function state "1"
 150   double betad=0.260;      // magic number for meson wave function state "2"
 151   double betasp=betas*betas+betap*betap;
 152   double betasd=betas*betas+betad*betad;
 153
 154   double lambdabar=0.750;  //M(0-,1-) - mQ From Goity&Roberts's code
 155
 156 // Isgur&Wise fct
 157   double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
 158   double xi1= -1.0*sqrt(2.0/3.0)*(
 159        lambdabar*lambdabar*(w*w-1.0)/(4*betas*betas))*
 160        exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
 161   double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
 162                pow((2*betas*betap/(betasp)),2.5)*
 163                exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
 164   double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
 165                pow((2*betas*betad/(betasd)),3.5)*
 166                exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
 167
 168   //report(INFO,"EvtGen") <<"rho's:"<<rho1<<rho2<<endl;
 169
 170   EvtComplex h1nr,h2nr,h3nr,f1nr,f2nr;
 171   EvtComplex f3nr,f4nr,f5nr,f6nr,knr,g1nr,g2nr,g3nr,g4nr,g5nr;
 172   EvtComplex h1r,h2r,h3r,f1r,f2r,f3r,f4r,f5r,f6r,kr,g1r,g2r,g3r,g4r,g5r;
 173   EvtComplex h1,h2,h3,f1,f2,f3,f4,f5,f6,k,g1,g2,g3,g4,g5;
 174
 175   // Non-resonance part
 176   h1nr = -g*xi*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmb));
 177   h2nr = -g*xi/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmb));
 178   h3nr = -(g*xi/(f0*md))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
 179                          -EvtComplex((1.0+w)/(p4_pi*vp),0.0));
 180
 181   f1nr = -(g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb) -
 182          1.0/(EvtComplex(p4_pi*vp,0.0)+dmd));
 183   f2nr = f1nr*mb/md;
 184   f3nr = EvtComplex(0.0);
 185   f4nr = EvtComplex(0.0);
 186   f5nr = (g*xi/(2*f0*mb*md))*(EvtComplex(1.0,0.0)
 187                  +(p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmb));
 188   f6nr = (g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
 189                  -EvtComplex(1.0/(p4_pi*vp),0.0));
 190
 191   knr = (g*xi/(2*f0))*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmb) +
 192                  EvtComplex((p4_pi*(v-w*vp))/(p4_pi*vp),0.0));
 193
 194   g1nr = EvtComplex(0.0);
 195   g2nr = EvtComplex(0.0);
 196   g3nr = EvtComplex(0.0);
 197   g4nr = (g*xi)/(f0*md*EvtComplex(p4_pi*vp));
 198   g5nr = EvtComplex(0.0);
 199
 200   // Resonance part (D** removed by hand - alainb)
 201   h1r = -alpha1*rho1*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
 202         alpha2*rho2*(p4_pi*(v+2.0*w*v-vp))
 203         /(3*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
 204         alpha3*xi1*(p4_pi*v)/(f0*mb*md*EvtComplex(p4_pi*v,0.0)+dmt3);
 205   h2r = -alpha2*(1+w)*rho2/(3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
 206         alpha3*xi1/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
 207   h3r = alpha2*rho2*(1+w)/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
 208         alpha3*xi1/(f0*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
 209
 210   f1r = -alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
 211         alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
 212   f2r = f1r*mb/md;
 213   f3r = EvtComplex(0.0);
 214   f4r = EvtComplex(0.0);
 215   f5r = alpha1*rho1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
 216         alpha2*rho2*(p4_pi*(vp-v/3.0-2.0/3.0*w*v))/
 217         (2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
 218         alpha3*xi1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
 219   f6r = alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
 220         alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
 221
 222   kr = -alpha1*rho1*(w-1.0)*(p4_pi*v)/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt1)) -
 223        alpha2*rho2*(w-1.0)*(p4_pi*(vp-w*v))
 224        /(3*f0*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
 225        alpha3*xi1*(p4_pi*(vp-w*v))/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt3));
 226
 227   g1r = EvtComplex(0.0);
 228   g2r = EvtComplex(0.0);
 229   g3r = -g2r;
 230   g4r = 2.0*alpha2*rho2/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2));
 231   g5r = EvtComplex(0.0);
 232
 233   //Sum
 234   h1 = h1nr + h1r;
 235   h2 = h2nr + h2r;
 236   h3 = h3nr + h3r;
 237
 238   f1 = f1nr + f1r;
 239   f2 = f2nr + f2r;
 240   f3 = f3nr + f3r;
 241   f4 = f4nr + f4r;
 242   f5 = f5nr + f5r;
 243   f6 = f6nr + f6r;
 244
 245   k = knr+kr;
 246
 247   g1 = g1nr + g1r;
 248   g2 = g2nr + g2r;
 249   g3 = g3nr + g3r;
 250   g4 = g4nr + g4r;
 251   g5 = g5nr + g5r;
 252
 253   EvtTensor4C g_metric;
 254   g_metric.setdiag(1.0,-1.0,-1.0,-1.0);
 255
 256   if (nlep==EM||nlep==MUM){
 257     omega=EvtComplex(0.0,0.5)*dual(h1*mb*md*EvtGenFunctions::directProd(v,vp)+
 258                              h2*mb*EvtGenFunctions::directProd(v,p4_pi)+
 259                              h3*md*EvtGenFunctions::directProd(vp,p4_pi))+
 260         f1*mb*EvtGenFunctions::directProd(v,p4_pi)+f2*md*EvtGenFunctions::directProd(vp,p4_pi)+
 261                        f3*EvtGenFunctions::directProd(p4_pi,p4_pi)+f4*mb*mb*EvtGenFunctions::directProd(v,v)+
 262         f5*mb*md*EvtGenFunctions::directProd(vp,v)+f6*mb*EvtGenFunctions::directProd(p4_pi,v)+k*g_metric+
 263         EvtComplex(0.0,0.5)*EvtGenFunctions::directProd(dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v),
 264                               (g1*p4_pi+g2*mb*v))+
 265         EvtComplex(0.0,0.5)*EvtGenFunctions::directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
 266                              dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v));
 267
 268    l1=EvtLeptonVACurrent(lepton->spParent(0),neutrino->spParentNeutrino());
 269    l2=EvtLeptonVACurrent(lepton->spParent(1),neutrino->spParentNeutrino());
 270   }
 271   else{
 272     if (nlep==EP||nlep==MUP){
 273       omega=EvtComplex(0.0,-0.5)*dual(h1*mb*md*EvtGenFunctions::directProd(v,vp)+
 274                              h2*mb*EvtGenFunctions::directProd(v,p4_pi)+
 275                                       h3*md*EvtGenFunctions::directProd(vp,p4_pi))+
 276         f1*mb*EvtGenFunctions::directProd(v,p4_pi)+f2*md*EvtGenFunctions::directProd(vp,p4_pi)+
 277                        f3*EvtGenFunctions::directProd(p4_pi,p4_pi)+f4*mb*mb*EvtGenFunctions::directProd(v,v)+
 278         f5*mb*md*EvtGenFunctions::directProd(vp,v)+f6*mb*EvtGenFunctions::directProd(p4_pi,v)+k*g_metric+
 279         EvtComplex(0.0,-0.5)*EvtGenFunctions::directProd(dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v),
 280                               (g1*p4_pi+g2*mb*v))+
 281         EvtComplex(0.0,-0.5)*EvtGenFunctions::directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
 282                              dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v));
 283
 284    l1=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(0));
 285    l2=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(1));
 286     }
 287     else{
 288    report(DEBUG,"EvtGen") << "42387dfs878w wrong lepton number\n";
 289     }
 290  }
 291
 292   et0=omega.cont2( dstar->epsParent(0).conj() );
 293   et1=omega.cont2( dstar->epsParent(1).conj() );
 294   et2=omega.cont2( dstar->epsParent(2).conj() );
 295
 296   vertex(0,0,l1.cont(et0));
 297   vertex(0,1,l2.cont(et0));
 298
 299   vertex(1,0,l1.cont(et1));
 300   vertex(1,1,l2.cont(et1));
 301
 302   vertex(2,0,l1.cont(et2));
 303   vertex(2,1,l2.cont(et2));
 304
 305   return;
 306
 307 }
 308
 309 void EvtGoityRoberts::DecayBDpilnuGR(EvtParticle *pb,EvtId nd,
 310                                      EvtId nlep, EvtId /*nnu*/)
 311
 312 {
 313   //added by Lange Jan4,2000
 314   static EvtId EM=EvtPDL::getId("e-");
 315   static EvtId EP=EvtPDL::getId("e+");
 316   static EvtId MUM=EvtPDL::getId("mu-");
 317   static EvtId MUP=EvtPDL::getId("mu+");
 318
 319   EvtParticle *d, *pion, *lepton, *neutrino;
 320
 321   pb->initializePhaseSpace(getNDaug(),getDaugs());
 322   d=pb->getDaug(0);
 323   pion=pb->getDaug(1);
 324   lepton=pb->getDaug(2);
 325   neutrino=pb->getDaug(3);
 326
 327   EvtVector4C l1, l2, et0, et1, et2;
 328
 329   EvtVector4R v,vp,p4_pi;
 330   double w;
 331
 332   v.set(1.0,0.0,0.0,0.0);       //4-velocity of B meson
 333   vp=(1.0/d->getP4().mass())*d->getP4();  //4-velocity of D
 334   p4_pi=pion->getP4();                  //4-momentum of pion
 335   w=v*vp;                       //four velocity transfer.
 336
 337   double mb=EvtPDL::getMeanMass(pb->getId());     //B mass
 338   double md=EvtPDL::getMeanMass(nd);   //D* mass
 339   EvtComplex dmb(0.0460,-0.5*0.00001);   //B mass splitting ?
 340                       //The last two numbers should be
 341                       //correctly calculated from the
 342                       //dstar and pion particle number.
 343
 344   double g = 0.5;         // Amplitude proportional to these coupling constants
 345   double alpha3 =  0.690; // See table I in G&R's paper
 346   double alpha1 = -1.430;
 347   double alpha2 = -0.140;
 348   double f0=0.093;        // The pion decay constant set to 93 MeV
 349
 350   EvtComplex dmt3(0.563,-0.5*0.191); // Mass splitting = dmt - iGamma/2
 351   EvtComplex dmt1(0.392,-0.5*1.040);
 352   EvtComplex dmt2(0.709,-0.5*0.405);
 353
 354   double betas=0.285;      // magic number for meson wave function ground state
 355   double betap=0.280;      // magic number for meson wave function state "1"
 356   double betad=0.260;      // magic number for meson wave function state "2"
 357   double betasp=betas*betas+betap*betap;
 358   double betasd=betas*betas+betad*betad;
 359
 360   double lambdabar=0.750;  //M(0-,1-) - mQ From Goity&Roberts's code
 361
 362   // Isgur&Wise fct
 363   double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
 364   double xi1= -1.0*sqrt(2.0/3.0)*(lambdabar*lambdabar*(w*w-1.0)/(4*betas*betas))*
 365               exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
 366   double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
 367                pow((2*betas*betap/(betasp)),2.5)*
 368                exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
 369   double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
 370                pow((2*betas*betad/(betasd)),3.5)*
 371                exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
 372
 373   EvtComplex h,a1,a2,a3;
 374   EvtComplex hnr,a1nr,a2nr,a3nr;
 375   EvtComplex hr,a1r,a2r,a3r;
 376
 377 // Non-resonance part (D* and D** removed by hand - alainb)
 378   hnr = g*xi*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb*md);
 379   a1nr= -1.0*g*xi*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0);
 380   a2nr= g*xi*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb);
 381   a3nr=EvtComplex(0.0,0.0);
 382
 383 // Resonance part (D** remove by hand - alainb)
 384   hr = alpha2*rho2*(w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0*mb*md) +
 385        alpha3*xi1*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb*md);
 386   a1r= -1.0*alpha2*rho2*(w*w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0) -
 387        alpha3*xi1*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0);
 388   a2r= alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*mb) +
 389        alpha2*rho2*(0.5*p4_pi*(w*vp-v)+p4_pi*(vp-w*v))/
 390                   (3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
 391        alpha3*xi1*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb);
 392   a3r= -1.0*alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*md) -
 393        alpha2*rho2*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmt2))/(2*f0*md);
 394
 395 // Sum
 396   h=hnr+hr;
 397   a1=a1nr+a1r;
 398   a2=a2nr+a2r;
 399   a3=a3nr+a3r;
 400
 401   EvtVector4C omega;
 402
 403   if ( nlep==EM|| nlep==MUM ) {
 404     omega=EvtComplex(0.0,-1.0)*h*mb*md*dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v)+
 405                  a1*p4_pi+a2*mb*v+a3*md*vp;
 406     l1=EvtLeptonVACurrent(
 407              lepton->spParent(0),neutrino->spParentNeutrino());
 408     l2=EvtLeptonVACurrent(
 409              lepton->spParent(1),neutrino->spParentNeutrino());
 410   }
 411   else{
 412     if ( nlep==EP|| nlep==MUP ) {
 413      omega=EvtComplex(0.0,1.0)*h*mb*md*dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v)+
 414                  a1*p4_pi+a2*mb*v+a3*md*vp;
 415      l1=EvtLeptonVACurrent(
 416               neutrino->spParentNeutrino(),lepton->spParent(0));
 417      l2=EvtLeptonVACurrent(
 418               neutrino->spParentNeutrino(),lepton->spParent(1));
 419     }
 420     else{
 421      report(ERROR,"EvtGen") << "42387dfs878w wrong lepton number\n";
 422     }
 423   }
 424
 425   vertex(0,l1*omega);
 426   vertex(1,l2*omega);
 427
 428 return;
 429
 430 }
 431