TEvtGen/EvtGenModels/EvtbTosllVectorAmp.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) 2000      Caltech, UCSB
  10 //
  11 // Module: EvtbTosllVectorAmp.cc
  12 //
  13 // Description: Routine to implement bTosll decays to vector
  14 //              mesons.
  15 //
  16 // Modification history:
  17 //
  18 //    Ryd       January 5,2000       Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include "EvtGenBase/EvtPatches.hh"
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtGenKine.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenBase/EvtVector4C.hh"
  29 #include "EvtGenBase/EvtTensor4C.hh"
  30 #include "EvtGenBase/EvtDiracSpinor.hh"
  31 #include "EvtGenModels/EvtbTosllVectorAmp.hh"
  32 #include "EvtGenBase/EvtId.hh"
  33 #include "EvtGenBase/EvtIdSet.hh"
  34 #include "EvtGenBase/EvtAmp.hh"
  35 #include "EvtGenModels/EvtbTosllAmp.hh"
  36 #include "EvtGenModels/EvtbTosllFF.hh"
  37
  38 void EvtbTosllVectorAmp::CalcAmp( EvtParticle *parent,
  39                                   EvtAmp& amp,
  40                                   EvtbTosllFF *formFactors ) {
  41
  42   //Add the lepton and neutrino 4 momenta to find q2
  43
  44   EvtVector4R q = parent->getDaug(1)->getP4()
  45                     + parent->getDaug(2)->getP4();
  46   double q2 = (q.mass2());
  47
  48
  49   double a1,a2,a0,v,t1,t2,t3;
  50   double mesonmass = parent->getDaug(0)->mass();
  51   double parentmass = parent->mass();
  52
  53   formFactors->getVectorFF(parent->getId(),
  54                            parent->getDaug(0)->getId(),
  55                            q2,
  56                            mesonmass,
  57                            a1,a2,a0,v,t1,t2,t3);
  58
  59
  60   EvtId daught = parent->getDaug(0)->getId();
  61   bool btod = false;
  62   bool nnlo = true;
  63  if
  64   (
  65    daught == EvtPDL::getId(std::string("rho+")) ||
  66    daught == EvtPDL::getId(std::string("rho-")) ||
  67    daught == EvtPDL::getId(std::string("rho0")) ||
  68    daught == EvtPDL::getId(std::string("omega"))
  69   )
  70    btod = true;
  71
  72   EvtVector4R p4b;
  73   p4b.set(parent->mass(),0.0,0.0,0.0);
  74   EvtVector4R p4meson = parent->getDaug(0)->getP4();
  75
  76   EvtVector4C l11,l12;
  77   EvtVector4C l21,l22;
  78
  79   EvtVector4C a11,a12;
  80   EvtVector4C a21,a22;
  81
  82   EvtId parentID = parent->getId();
  83
  84       //EvtId l_num = parent->getDaug(1)->getId();
  85
  86   EvtVector4R pbhat=p4b/parentmass;
  87   EvtVector4R qhat=q/parentmass;
  88   EvtVector4R pkstarhat=p4meson/parentmass;
  89   EvtVector4R phat=pbhat+pkstarhat;
  90
  91   EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(q2,nnlo);
  92   EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(q2,nnlo,btod);
  93   EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(q2,nnlo);
  94   EvtComplex uniti(0.0,1.0);
  95
  96   double mhatb=4.4/(parentmass);
  97   double mhatkstar=mesonmass/(parentmass);
  98   double shat=q2/(parentmass*parentmass);
  99
 100
 101   EvtComplex a;
 102   a=c9eff*v*2/(1+mhatkstar)+4*mhatb*c7eff*t1/shat;
 103   EvtComplex b;
 104   b=(1+mhatkstar)*(c9eff*a1+2*mhatb*(1-mhatkstar)*c7eff*t2/shat);
 105   EvtComplex c;
 106   c=((1-mhatkstar)*c9eff*a2+
 107             2*mhatb*c7eff*(t3+(1-mhatkstar*mhatkstar)*t2/shat))/
 108     (1-mhatkstar*mhatkstar);
 109   EvtComplex d;
 110   d=(c9eff*((1+mhatkstar)*a1-(1-mhatkstar)*a2-2*mhatkstar*a0)
 111             -2*mhatb*c7eff*t3)/shat;
 112   EvtComplex e;
 113   e=2*c10eff*v/(1+mhatkstar);
 114   EvtComplex f;
 115   f=(1+mhatkstar)*c10eff*a1;
 116   EvtComplex g;
 117   g=c10eff*a2/(1+mhatkstar);
 118   EvtComplex h;
 119   h=c10eff*((1+mhatkstar)*a1-(1-mhatkstar)*a2-2*mhatkstar*a0)/shat;
 120
 121   EvtTensor4C T1,T2;
 122
 123   static EvtIdSet bmesons("B-","anti-B0","anti-B_s0");
 124   static EvtIdSet bbarmesons("B+","B0","B_s0");
 125
 126   EvtParticle* lepPlus=0;
 127   EvtParticle* lepMinus=0;
 128
 129   int charge1 = EvtPDL::chg3(parent->getDaug(1)->getId());
 130   int charge2 = EvtPDL::chg3(parent->getDaug(2)->getId());
 131
 132   lepPlus = (charge1 > charge2) ? parent->getDaug(1) : parent->getDaug(2);
 133   lepMinus = (charge1 < charge2) ? parent->getDaug(1) : parent->getDaug(2);
 134
 135   if (bmesons.contains(parentID)) {
 136
 137     T1=a*dual(directProd(pbhat,pkstarhat))
 138        -b*uniti*EvtTensor4C::g()
 139        +c*uniti*directProd(pbhat,phat)
 140        +d*uniti*directProd(pbhat,qhat);
 141
 142     T2=e*dual(directProd(pbhat,pkstarhat))
 143        -f*uniti*EvtTensor4C::g()
 144        +g*uniti*directProd(pbhat,phat)
 145        +h*uniti*directProd(pbhat,qhat);
 146
 147     l11=EvtLeptonVCurrent(lepPlus->spParent(0),
 148                           lepMinus->spParent(0));
 149     l21=EvtLeptonVCurrent(lepPlus->spParent(1),
 150                           lepMinus->spParent(0));
 151     l12=EvtLeptonVCurrent(lepPlus->spParent(0),
 152                           lepMinus->spParent(1));
 153     l22=EvtLeptonVCurrent(lepPlus->spParent(1),
 154                           lepMinus->spParent(1));
 155
 156     a11=EvtLeptonACurrent(lepPlus->spParent(0),
 157                           lepMinus->spParent(0));
 158     a21=EvtLeptonACurrent(lepPlus->spParent(1),
 159                           lepMinus->spParent(0));
 160     a12=EvtLeptonACurrent(lepPlus->spParent(0),
 161                           lepMinus->spParent(1));
 162     a22=EvtLeptonACurrent(lepPlus->spParent(1),
 163                           lepMinus->spParent(1));
 164
 165   } else {
 166
 167     if (bbarmesons.contains(parentID)) {
 168
 169       T1=-a*dual(directProd(pbhat,pkstarhat))
 170          -b*uniti*EvtTensor4C::g()
 171          +c*uniti*directProd(pbhat,phat)
 172          +d*uniti*directProd(pbhat,qhat);
 173
 174       T2=-e*dual(directProd(pbhat,pkstarhat))
 175          -f*uniti*EvtTensor4C::g()
 176          +g*uniti*directProd(pbhat,phat)
 177          +h*uniti*directProd(pbhat,qhat);
 178
 179       l11=EvtLeptonVCurrent(lepPlus->spParent(1),
 180                             lepMinus->spParent(1));
 181       l21=EvtLeptonVCurrent(lepPlus->spParent(0),
 182                             lepMinus->spParent(1));
 183       l12=EvtLeptonVCurrent(lepPlus->spParent(1),
 184                             lepMinus->spParent(0));
 185       l22=EvtLeptonVCurrent(lepPlus->spParent(0),
 186                             lepMinus->spParent(0));
 187
 188       a11=EvtLeptonACurrent(lepPlus->spParent(1),
 189                             lepMinus->spParent(1));
 190       a21=EvtLeptonACurrent(lepPlus->spParent(0),
 191                             lepMinus->spParent(1));
 192       a12=EvtLeptonACurrent(lepPlus->spParent(1),
 193                             lepMinus->spParent(0));
 194       a22=EvtLeptonACurrent(lepPlus->spParent(0),
 195                             lepMinus->spParent(0));
 196
 197     }
 198     else{
 199       report(ERROR,"EvtGen") << "Wrong lepton number\n";
 200     }
 201   }
 202
 203
 204   int i;
 205
 206   for(i=0;i<3;i++){
 207     EvtVector4C eps=parent->getDaug(0)->epsParent(i).conj();
 208
 209     EvtVector4C E1=T1.cont1(eps);
 210     EvtVector4C E2=T2.cont1(eps);
 211
 212     amp.vertex(i,0,0,l11*E1+a11*E2);
 213     amp.vertex(i,0,1,l12*E1+a12*E2);
 214     amp.vertex(i,1,0,l21*E1+a21*E2);
 215     amp.vertex(i,1,1,l22*E1+a22*E2);
 216   }
 217 }
 218
 219
 220
 221
 222
 223
 224