00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include "EvtGenBase/EvtPatches.hh"
00038 #include "EvtGenModels/EvtVubdGamma.hh"
00039
00040
00041
00042
00043 #include <math.h>
00044
00045 extern "C" {
00046 double ddilog_(const double *);
00047 }
00048
00049
00050
00051
00052 EvtVubdGamma::EvtVubdGamma(const double &alphas)
00053 {
00054 _alphas = alphas;
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068 _epsilon1 = 1e-10;
00069 _epsilon2 = 1e-5;
00070 if ( alphas > 0 ) {
00071 double lne3 = 9./16.-2*M_PI*M_PI/3.+6*M_PI/4/alphas;
00072 if ( lne3 > 0 )
00073 lne3 = -7./4. - sqrt(lne3);
00074 else
00075 lne3 = -7./4.;
00076 _epsilon3 = exp(lne3);
00077 }
00078 else
00079 _epsilon3 = 1;
00080 }
00081
00082
00083
00084
00085
00086 EvtVubdGamma::~EvtVubdGamma( )
00087 {
00088 }
00089
00090
00091
00092
00093
00094 double EvtVubdGamma::getdGdxdzdp(const double &x, const double &z, const double &p2)
00095 {
00096
00097
00098 double xb = (1-x);
00099
00100 if ( x < 0 || x > 1 || z < xb || z > (1+xb) )
00101 return 0;
00102
00103
00104
00105
00106 double p2min = (0>z-1.?0:z-1.);
00107 double p2max = (1.-x)*(z-1.+x);
00108
00109 if (p2 < p2min || p2 > p2max)
00110 return 0;
00111
00112
00113
00114
00115 double dG;
00116
00117 if ( p2 >_epsilon1 && p2< _epsilon2) {
00118
00119 double W1 = getW1delta(x,z);
00120 double W4plus5 = getW4plus5delta(x,z);
00121
00122 dG = 12. * delta(p2,p2min,p2max) * ((1.+xb-z) * (z-xb) * W1
00123 + xb*(z-xb) * (W4plus5));
00124 }
00125 else {
00126
00127 double W1 = getW1nodelta(x,z,p2);
00128 double W2 = getW2nodelta(x,z,p2);
00129 double W3 = getW3nodelta(x,z,p2);
00130 double W4 = getW4nodelta(x,z,p2);
00131 double W5 = getW5nodelta(x,z,p2);
00132
00133 dG = 12. * ((1.+xb-z) * (z-xb-p2) * W1
00134 + (1.-z+p2) * W2
00135 + (xb*(z-xb)-p2) * (W3+W4+W5));
00136 }
00137 return dG;
00138 }
00139
00140 double EvtVubdGamma::delta(const double &x, const double &xmin, const double &xmax)
00141 {
00142 if ( xmin > 0 || xmax < 0 ) return 0.;
00143 if ( _epsilon1 < x && x < _epsilon2 ) return 1./(_epsilon2-_epsilon1);
00144 return 0.0;
00145 }
00146
00147 double EvtVubdGamma::getW1delta(const double &x, const double &z)
00148 {
00149 double mz = 1.-z;
00150
00151
00152
00153 double lz;
00154 if (z == 1) lz = -1.;
00155 else lz = log(z)/(1.-z);
00156
00157
00158
00159
00160
00161
00162
00163 double dl = 4.*ddilog_(&mz) + 4.*pow(M_PI,2)/3.;
00164
00165 double w = -(8.*pow(log(z),2) - 10.*log(z) + 2.*lz + dl + 5.)
00166 + (8.*log(z)-7.)*log(_epsilon3) - 2.*pow(log(_epsilon3),2);
00167
00168 return (1. + w*_alphas/3./M_PI);
00169 }
00170
00171 double EvtVubdGamma::getW1nodelta(const double &x, const double &z, const double &p2)
00172 {
00173
00174 double z2 = z*z;
00175 double t2 = 1.-4.*p2/z2;
00176 double t = sqrt(t2);
00177
00178 double w = 0;
00179 if ( p2 > _epsilon2 )
00180 w += 4./p2*(log((1.+t)/(1.-t))/t + log(p2/z2))
00181 + 1. - (8.-z)*(2.-z)/z2/t2
00182 + ((2.-z)/2./z+(8.-z)*(2.-z)/2./z2/t2)*log((1.+t)/(1.-t))/t;
00183 if ( p2 > _epsilon3 )
00184 w += (8.*log(z)-7.)/p2 - 4.*log(p2)/p2;
00185
00186 return w*_alphas/3./M_PI;
00187 }
00188
00189 double EvtVubdGamma::getW2nodelta(const double &x, const double &z, const double &p2)
00190 {
00191
00192 double z2 = z*z;
00193 double t2 = 1.-4.*p2/z2;
00194 double t = sqrt(t2);
00195 double w11 = (32.-8.*z+z2)/4./z/t2;
00196
00197 double w = 0;
00198 if ( p2 > _epsilon2 )
00199 w -= (z*t2/8. + (4.-z)/4. + w11/2.)*log((1.+t)/(1.-t))/t;
00200 if ( p2 > _epsilon2 )
00201 w += (8.-z)/4. + w11;
00202
00203 return (w*_alphas/3./M_PI);
00204 }
00205
00206 double EvtVubdGamma::getW3nodelta(const double &x, const double &z, const double &p2)
00207 {
00208 double z2 = z*z;
00209 double t2 = 1.-4.*p2/z2;
00210 double t4 = t2*t2;
00211 double t = sqrt(t2);
00212
00213 double w = 0;
00214
00215 if ( p2 > _epsilon2 )
00216 w += (z*t2/16. + 5.*(4.-z)/16. - (64.+56.*z-7.*z2)/16./z/t2
00217 + 3.*(12.-z)/16./t4) * log((1.+t)/(1.-t))/t;
00218 if ( p2 > _epsilon2 )
00219 w += -(8.-3.*z)/8. + (32.+22.*z-3.*z2)/4./z/t2 - 3.*(12.-z)/8./t4;
00220
00221 return (w*_alphas/3./M_PI);
00222 }
00223
00224 double EvtVubdGamma::getW4nodelta(const double &x, const double &z, const double &p2)
00225 {
00226 double z2 = z*z;
00227 double t2 = 1.-4.*p2/z2;
00228 double t4 = t2*t2;
00229 double t = sqrt(t2);
00230
00231 double w = 0;
00232
00233 if ( p2 > _epsilon2 )
00234 w -= ((8.-3.*z)/4./z - (22.-3.*z)/2./z/t2 + 3.*(12.-z)/4./z/t4)
00235 * log((1.+t)/(1.-t))/t;
00236 if ( p2 > _epsilon2 )
00237 w += -1. - (32.-5.*z)/2./z/t2 + 3.*(12.-z)/2./z/t4 ;
00238
00239 return w*_alphas/3./M_PI;
00240 }
00241
00242 double EvtVubdGamma::getW4plus5delta(const double &x, const double &z)
00243 {
00244
00245 double w = 0;
00246
00247 if ( z == 1 )
00248 w = -2;
00249 else
00250 w = 2.*log(z)/(1.-z);
00251
00252 return (w*_alphas/3./M_PI);
00253 }
00254
00255 double EvtVubdGamma::getW5nodelta(const double &x, const double &z, const double &p2)
00256 {
00257 double z2 = z*z;
00258 double t2 = 1.-4.*p2/z2;
00259 double t4 = t2*t2;
00260 double t = sqrt(t2);
00261
00262 double w = 0;
00263 if ( p2 > _epsilon2 )
00264 w += (1./4./z - (2.-z)/2./z2/t2 + 3.*(12.-z)/4./z2/t4)
00265 * log((1.+t)/(1.-t))/t;
00266 if ( p2 > _epsilon2 )
00267 w += -(8.+z)/2./z2/t2 - 3.*(12.-z)/2./z2/t4;
00268
00269
00270 return (w*_alphas/3./M_PI);
00271 }
00272
00273
00274