TEvtGen/EvtGenModels/EvtBsquark.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: EvtTauScalarnu.cc
  12 //
  13 // Description: The leptonic decay of the tau meson.
  14 //              E.g., tau- -> e- nueb nut
  15 //
  16 // Modification history:
  17 //
  18 //    RYD       January 17, 1997       Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include "EvtGenBase/EvtPatches.hh"
  24 #include <iostream>
  25 #include <string>
  26 #include "EvtGenBase/EvtParticle.hh"
  27 #include "EvtGenBase/EvtDiracParticle.hh"
  28 #include "EvtGenBase/EvtPDL.hh"
  29 #include "EvtGenBase/EvtIdSet.hh"
  30 #include "EvtGenBase/EvtGenKine.hh"
  31 #include "EvtGenModels/EvtBsquark.hh"
  32 #include "EvtGenBase/EvtDiracSpinor.hh"
  33 #include "EvtGenBase/EvtGammaMatrix.hh"
  34 #include "EvtGenBase/EvtReport.hh"
  35
  36 EvtBsquark::~EvtBsquark() {}
  37
  38 std::string EvtBsquark::getName(){
  39
  40   return "BSQUARK";
  41
  42 }
  43
  44
  45 EvtDecayBase* EvtBsquark::clone(){
  46
  47   return new EvtBsquark;
  48
  49 }
  50
  51 void EvtBsquark::init(){
  52
  53   // check that there are 5 arguments
  54   checkNArg(5);
  55 }
  56
  57 void EvtBsquark::initProbMax(){
  58
  59   //For now do not set a maximum.
  60
  61   //SetProbMax(0.000000000005);
  62
  63 }
  64
  65 void EvtBsquark::decay(EvtParticle *p){
  66
  67   static EvtId cquark=EvtPDL::getId("c");
  68   static EvtId anticquark=EvtPDL::getId("anti-c");
  69
  70   static EvtIdSet leptons("e-","mu-","tau-");
  71
  72   p->initializePhaseSpace(getNDaug(),getDaugs());
  73
  74
  75   int charge=1;
  76
  77   EvtParticle* lepton;
  78   lepton = p->getDaug(1);
  79   if (leptons.contains(lepton->getId())){
  80     charge=-1;
  81   }
  82
  83
  84   EvtDiracParticle charmquark;
  85
  86   //this is a very crude approximation...
  87   if (charge==-1){
  88     charmquark.init(cquark,p->getDaug(0)->getP4());
  89   }
  90   else{
  91     charmquark.init(anticquark,p->getDaug(0)->getP4());
  92   }
  93
  94   EvtVector4R p4c = p->getDaug(0)->getP4();
  95
  96   EvtVector4R p4sn = p->getDaug(2)->getP4();
  97
  98   EvtVector4R p4b(p->mass(),0.0,0.0,0.0);
  99
 100   EvtComplex M[2][2];
 101
 102   int il,ic;
 103
 104   //project out the right handed current
 105   EvtGammaMatrix PR=0.5*(EvtGammaMatrix::id()+EvtGammaMatrix::g5());
 106
 107   double tanbeta=getArg(1);
 108   double cosbeta=cos(atan(tanbeta));
 109   double sinbeta=sin(atan(tanbeta));
 110
 111   double mb=4.9;
 112   double mc=1.3;
 113   double mw=80.4;
 114
 115   double Mass=getArg(2);
 116   double mu=getArg(3);
 117   double mchargino=getArg(4);
 118
 119
 120   double tan2phim=2*sqrt(2.0)*mw*(mu*cosbeta+Mass*sinbeta)/
 121     (Mass*Mass-mu*mu+2*mw*mw*cos(2*atan(tanbeta)));
 122
 123   double phim=0.5*atan(tan2phim);
 124
 125   EvtComplex U11=cos(phim);
 126   EvtComplex U12=sin(phim);
 127   EvtComplex U21=-sin(phim);
 128   EvtComplex U22=cos(phim);
 129
 130   double tan2phip=2*sqrt(2.0)*mw*(mu*cosbeta+Mass*sinbeta)/
 131     (Mass*Mass-mu*mu-2*mw*mw*cos(2*atan(tanbeta)));
 132
 133   double phip=0.5*atan(tan2phip);
 134
 135   EvtComplex V11=cos(phip);
 136   EvtComplex V12=sin(phip);
 137   EvtComplex V21=-sin(phip);
 138   EvtComplex V22=cos(phip);
 139
 140
 141   double theta=getArg(0);
 142   double ctheta=cos(theta);
 143   double stheta=sin(theta);
 144
 145   double vcsb=0.08;
 146   double mchi1=mchargino;
 147   double mchi2=mchargino;
 148
 149   //overall scale factor
 150   double g=1.0;
 151
 152   EvtComplex a1=mchi1*(U11*ctheta-mb*U12*stheta/(sqrt(2.0)*mw*cosbeta));
 153   EvtComplex a2=mchi2*(U21*ctheta-mb*U22*stheta/(sqrt(2.0)*mw*cosbeta));
 154
 155   EvtComplex b1=mc*conj(V12)*ctheta/(sqrt(2.0)*mw*sinbeta);
 156   EvtComplex b2=mc*conj(V22)*ctheta/(sqrt(2.0)*mw*sinbeta);
 157
 158   EvtComplex f1=-(g*g*V11*vcsb)/((p4b-p4c).mass2()-mchi1*mchi1);
 159   EvtComplex f2=-(g*g*V21*vcsb)/((p4b-p4c).mass2()-mchi1*mchi2);
 160
 161   //report(INFO,"EvtGen") <<g<<" "<<V11<<" "<<FL<<" "<<vcsb<<" "<<mchi1<<endl;
 162   //report(INFO,"EvtGen") << "f1:"<<f1<<" "<<(p4b-p4c).mass2()<<endl;
 163   //report(INFO,"EvtGen") << "f2:"<<f2<<" "<<(p4b-p4c).mass2()<<endl;
 164
 165   //report(INFO,"EvtGen") << "p4sn:"<<p4sn<<endl;
 166
 167   EvtGammaMatrix pslash=p4sn.get(0)*EvtGammaMatrix::g0()
 168                        -p4sn.get(1)*EvtGammaMatrix::g1()
 169                        -p4sn.get(2)*EvtGammaMatrix::g2()
 170                        -p4sn.get(3)*EvtGammaMatrix::g3();
 171
 172   //report(INFO,"EvtGen") << "pslash:"<<pslash<<endl;
 173
 174
 175
 176   for(il=0;il<2;il++){
 177     for(ic=0;ic<2;ic++){
 178
 179       EvtComplex a=0.0;
 180       EvtComplex b=0.0;
 181
 182       if (charge==-1){
 183         a=charmquark.spParent(ic)*(PR*lepton->spParent(il));
 184         b=charmquark.spParent(ic)*((pslash*PR)*lepton->spParent(il));
 185       }
 186       else{
 187         a=lepton->spParent(il)*(PR*charmquark.spParent(ic));
 188         b=lepton->spParent(il)*((pslash*PR)*charmquark.spParent(ic));
 189       }
 190
 191       //report(INFO,"EvtGen") <<"pslash*PR:"<<pslash*PR<<endl;
 192       //report(INFO,"EvtGen") <<"sp charm:"<<charmquark.spParent(ic)<<endl;
 193       //report(INFO,"EvtGen") <<"sp lepton:"<<lepton->spParent(il)<<endl;
 194
 195       M[ic][il]=f1*(a1*a+b1*b)+f2*(a2*a+b2*b);
 196
 197       //report(INFO,"EvtGen") << "Contr1:"<<a1<<" "<<a<<" "<<b1<<" "<<b<<endl;
 198       //report(INFO,"EvtGen") << "Contr2:"<<a2<<" "<<a<<" "<<b2<<" "<<b<<endl;
 199
 200       //report(INFO,"EvtGen") <<"case1:"<<f1<<" "<<a1<<" "<<b1<<" "<<a<<" "<<b<<endl;
 201       //report(INFO,"EvtGen") <<"case2:"<<f2<<" "<<a2<<" "<<b2<<" "<<a<<" "<<b<<endl;
 202
 203     }
 204   }
 205
 206   double prob=real(M[0][0]*conj(M[0][0])+
 207                    M[1][0]*conj(M[1][0])+
 208                    M[0][1]*conj(M[0][1])+
 209                    M[1][1]*conj(M[1][1]));
 210
 211   //report(INFO,"EvtGen") <<"prob:"<<prob<<endl;
 212
 213   setProb(prob);
 214
 215   return;
 216
 217 }
 218