TEvtGen/EvtGenModels/EvtVubdGamma.cpp

   1 //-----------------------------------------------------------------------
   2 // File and Version Information:
   3 //
   4 // Copyright Information: See EvtGen/COPYRIGHT
   5 //
   6 //
   7 // Description:
   8 //    3                                         2                2
   9 //   d Gamma                 /   _      _ _2  mb           _2  mb
  10 //  ---------- = 12 Gamma   | (1+x-z)(z-x-p ) -- W  + (1-z+p ) -- W
  11 //         _ 2           0  \                  2  1             2  2
  12 //  dx dz dp                                   2
  13 //                                _   _  _2  mb                 2   \.
  14 //                             + [x(z-x)-p ] -- (W + 2mb W  + mb W ) |
  15 //                                            4   3       4       5 /
  16 //
  17 //   with
  18 //        2 E           2
  19 //           l    _2   p        2 v.p    _
  20 //   x = ------ , p = --- , z = ------ , x = 1-x
  21 //         mb           2         mb
  22 //                    mb
  23 //
  24 //   the triple differential decay rate according to
  25 //   hep-ph/9905351 v2
  26 //
  27 // Environment:
  28 //      Software developed for the BaBar Detector at the SLAC B-Factory.
  29 //
  30 // Author List:
  31 //      Sven Menke
  32 //
  33 //-----------------------------------------------------------------------
  34 //-----------------------
  35 // This Class's Header --
  36 //-----------------------
  37 #include "EvtGenBase/EvtPatches.hh"
  38 #include "EvtGenBase/EvtConst.hh"
  39 #include "EvtGenModels/EvtVubdGamma.hh"
  40 #include "EvtGenBase/EvtDiLog.hh"
  41
  42 //---------------
  43 // C Headers --
  44 //---------------
  45 #include <math.h>
  46
  47 //----------------
  48 // Constructors --
  49 //----------------
  50
  51 EvtVubdGamma::EvtVubdGamma(const double &alphas)
  52 {
  53   _alphas = alphas;
  54
  55   // the range for the delta distribution in p2 is from _epsilon1 to
  56   // _epsilon2. It was checked with the single differential formulae
  57   // in the paper that these values are small enough to imitate p2 = 0
  58   // for the regular terms.
  59   // The ()* distributions, however need further treatment. In order to
  60   // generate the correct spectrum in z a threshold need to be computed
  61   // from the desired value of the coupling alphas. The idea is that
  62   // for z=1 p2=0 is not allowed and therefore the part of dGamma proportional
  63   // to delta(p2) should go to 0 for z->1.
  64   // Using equation (3.1) and (3.2) it is possible to find the correct value
  65   // for log(_epsilon3) from this requirement.
  66
  67   _epsilon1 = 1e-10;
  68   _epsilon2 = 1e-5;
  69   if ( alphas > 0 ) {
  70     double lne3 = 9./16.-2*EvtConst::pi*EvtConst::pi/3.+6*EvtConst::pi/4/alphas;    if ( lne3 > 0 )
  71       lne3 = -7./4. - sqrt(lne3);
  72     else
  73       lne3 = -7./4.;
  74     _epsilon3 = exp(lne3);
  75   }
  76   else
  77     _epsilon3 = 1;
  78 }
  79
  80 //--------------
  81 // Destructor --
  82 //--------------
  83
  84 EvtVubdGamma::~EvtVubdGamma( )
  85 {
  86 }
  87
  88 //-----------
  89 // Methods --
  90 //-----------
  91
  92 double EvtVubdGamma::getdGdxdzdp(const double &x, const double &z, const double &p2)
  93 {
  94   // check phase space
  95
  96   double xb = (1-x);
  97
  98   if ( x < 0 || x > 1 || z < xb || z > (1+xb) )
  99     return 0;
 100
 101   double p2min = (0>z-1.?0:z-1.);
 102   double p2max = (1.-x)*(z-1.+x);
 103
 104   if (p2 < p2min || p2 > p2max)
 105     return 0;
 106
 107   //  // check the phase space
 108   //  return 1.;
 109
 110   double dG;
 111
 112   if ( p2 >_epsilon1 && p2< _epsilon2) {
 113
 114     double W1      = getW1delta(x,z);
 115     double W4plus5 = getW4plus5delta(x,z);
 116
 117     dG = 12. * delta(p2,p2min,p2max) * ((1.+xb-z) * (z-xb) * W1
 118                                         + xb*(z-xb) * (W4plus5));
 119   }
 120   else {
 121
 122     double W1 = getW1nodelta(x,z,p2);
 123     double W2 = getW2nodelta(x,z,p2);
 124     double W3 = getW3nodelta(x,z,p2);
 125     double W4 = getW4nodelta(x,z,p2);
 126     double W5 = getW5nodelta(x,z,p2);
 127
 128     dG = 12. * ((1.+xb-z) * (z-xb-p2) * W1
 129                 + (1.-z+p2) * W2
 130                 + (xb*(z-xb)-p2) * (W3+W4+W5));
 131   }
 132   return dG;
 133 }
 134
 135 double EvtVubdGamma::delta(const double &x, const double &xmin, const double &xmax)
 136 {
 137   if ( xmin > 0 || xmax < 0 ) return 0.;
 138   if ( _epsilon1 < x && x < _epsilon2 ) return 1./(_epsilon2-_epsilon1);
 139   return 0.0;
 140 }
 141
 142 double EvtVubdGamma::getW1delta(const double &, const double &z)
 143 {
 144   double mz = 1.-z;
 145
 146   double lz;
 147   if (z == 1) lz = -1.;
 148   else        lz = log(z)/(1.-z);
 149
 150   // ddilog_(&z) is actually the dilog of (1-z) in maple,
 151   // also in Neuberts paper the limit dilog(1) = pi^2/6 is used
 152   // this corresponds to maple's dilog(0), so
 153   // I take ddilog_(&mz) where mz=1-z in order to satisfy Neubert's definition
 154   // and to compare with Maple the argument in maple should be (1-mz) ...
 155
 156   double dl = 4.*EvtDiLog::DiLog(mz) + 4.*pow(EvtConst::pi,2)/3.;
 157
 158   double w = -(8.*pow(log(z),2) - 10.*log(z) + 2.*lz + dl + 5.)
 159     + (8.*log(z)-7.)*log(_epsilon3) - 2.*pow(log(_epsilon3),2);
 160
 161   return (1. + w*_alphas/3./EvtConst::pi);
 162 }
 163
 164 double EvtVubdGamma::getW1nodelta(const double &, const double &z, const double &p2)
 165 {
 166
 167   double z2 = z*z;
 168   double t2 = 1.-4.*p2/z2;
 169   double t = sqrt(t2);
 170
 171   double w = 0;
 172   if ( p2 > _epsilon2 )
 173     w += 4./p2*(log((1.+t)/(1.-t))/t + log(p2/z2))
 174       + 1. - (8.-z)*(2.-z)/z2/t2
 175       + ((2.-z)/2./z+(8.-z)*(2.-z)/2./z2/t2)*log((1.+t)/(1.-t))/t;
 176   if ( p2 > _epsilon3 )
 177     w += (8.*log(z)-7.)/p2 - 4.*log(p2)/p2;
 178
 179   return w*_alphas/3./EvtConst::pi;
 180 }
 181
 182 double EvtVubdGamma::getW2nodelta(const double &, const double &z, const double &p2)
 183 {
 184
 185   double z2 = z*z;
 186   double t2 = 1.-4.*p2/z2;
 187   double t = sqrt(t2);
 188   double w11 = (32.-8.*z+z2)/4./z/t2;
 189
 190   double w = 0;
 191   if ( p2 > _epsilon2 )
 192     w -=  (z*t2/8. + (4.-z)/4. + w11/2.)*log((1.+t)/(1.-t))/t;
 193   if ( p2 > _epsilon2 )
 194     w += (8.-z)/4. + w11;
 195
 196   return (w*_alphas/3./EvtConst::pi);
 197 }
 198
 199 double EvtVubdGamma::getW3nodelta(const double &, const double &z, const double &p2)
 200 {
 201   double z2 = z*z;
 202   double t2 = 1.-4.*p2/z2;
 203   double t4 = t2*t2;
 204   double t = sqrt(t2);
 205
 206   double w = 0;
 207
 208   if ( p2 > _epsilon2 )
 209     w += (z*t2/16. + 5.*(4.-z)/16. - (64.+56.*z-7.*z2)/16./z/t2
 210           + 3.*(12.-z)/16./t4) * log((1.+t)/(1.-t))/t;
 211   if ( p2 > _epsilon2 )
 212     w += -(8.-3.*z)/8. + (32.+22.*z-3.*z2)/4./z/t2 - 3.*(12.-z)/8./t4;
 213
 214   return (w*_alphas/3./EvtConst::pi);
 215 }
 216
 217 double EvtVubdGamma::getW4nodelta(const double &, const double &z, const double &p2)
 218 {
 219   double z2 = z*z;
 220   double t2 = 1.-4.*p2/z2;
 221   double t4 = t2*t2;
 222   double t = sqrt(t2);
 223
 224   double w = 0;
 225
 226   if ( p2 > _epsilon2 )
 227     w -= ((8.-3.*z)/4./z - (22.-3.*z)/2./z/t2 + 3.*(12.-z)/4./z/t4)
 228       * log((1.+t)/(1.-t))/t;
 229   if ( p2 > _epsilon2 )
 230     w += -1. - (32.-5.*z)/2./z/t2 + 3.*(12.-z)/2./z/t4 ;
 231
 232   return w*_alphas/3./EvtConst::pi;
 233 }
 234
 235 double EvtVubdGamma::getW4plus5delta(const double &, const double &z)
 236 {
 237
 238   double w = 0;
 239
 240   if ( z == 1 )
 241     w = -2;
 242   else
 243     w = 2.*log(z)/(1.-z);
 244
 245   return (w*_alphas/3./EvtConst::pi);
 246 }
 247
 248 double EvtVubdGamma::getW5nodelta(const double &, const double &z, const double &p2)
 249 {
 250   double z2 = z*z;
 251   double t2 = 1.-4.*p2/z2;
 252   double t4 = t2*t2;
 253   double t = sqrt(t2);
 254
 255   double w = 0;
 256   if ( p2 > _epsilon2 )
 257     w += (1./4./z - (2.-z)/2./z2/t2 + 3.*(12.-z)/4./z2/t4)
 258       * log((1.+t)/(1.-t))/t;
 259   if ( p2 > _epsilon2 )
 260     w += -(8.+z)/2./z2/t2 - 3.*(12.-z)/2./z2/t4;
 261
 262   return (w*_alphas/3./EvtConst::pi);
 263
 264 }
 265
 266
 267