00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022
00023
00024
00025
00026 #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
00027 #include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
00028 #include "EvtGenModels/EvtBtoXsgammaRootFinder.hh"
00029 #include "EvtGenModels/EvtItgFunction.hh"
00030 #include "EvtGenBase/EvtConst.hh"
00031 #include "EvtGenBase/EvtReport.hh"
00032
00033
00034
00035
00036 #include <iostream>
00037 #include <math.h>
00038 using std::endl;
00039
00040 double EvtBtoXsgammaFermiUtil::FermiExpFunc(double y, const std::vector<double> &coeffs) {
00041
00042
00043
00044 return (pow(1. - (y/coeffs[1]),coeffs[2])*exp((-3.*pow(coeffs[1],2.)/coeffs[3])*y/coeffs[1]))/coeffs[4];
00045
00046 }
00047
00048 double EvtBtoXsgammaFermiUtil::FermiGaussFunc(double y, const std::vector<double> &coeffs) {
00049
00050
00051 return (pow(1. - (y/coeffs[1]),coeffs[2])*exp(-pow(coeffs[3],2.)*pow(1. - (y/coeffs[1]),2.)))/coeffs[4];
00052
00053 }
00054
00055 double EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot(double lambdabar, double lam1, double mb, std::vector<double> &gammaCoeffs) {
00056
00057 std::vector<double> coeffs1(3);
00058 std::vector<double> coeffs2(3);
00059
00060 coeffs1[0]=0.2;
00061 coeffs1[1]=lambdabar;
00062 coeffs1[2]=0.0;
00063
00064 coeffs2[0]=0.2;
00065 coeffs2[1]=lambdabar;
00066 coeffs2[2]=-lam1/3.;
00067
00068 EvtItgTwoCoeffFcn *lhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnA, -mb, lambdabar, coeffs1, gammaCoeffs);
00069 EvtItgTwoCoeffFcn *rhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnB, -mb, lambdabar, coeffs2, gammaCoeffs);
00070
00071 EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder();
00072
00073 double root = rootFinder->GetGaussIntegFcnRoot(lhFunc, rhFunc, 1.0e-4, 1.0e-4, 40, 40, -mb, lambdabar, 0.2, 0.4, 1.0e-6);
00074
00075 delete rootFinder; rootFinder=0;
00076 delete lhFunc; lhFunc=0;
00077 delete rhFunc; rhFunc=0;
00078 return root;
00079
00080 }
00081
00082 double EvtBtoXsgammaFermiUtil::FermiGaussRootFcnA(double y, const std::vector<double> &coeffs1, const std::vector<double> &coeffs2) {
00083
00084
00085
00086 double cp = Gamma((2.0 + coeffs1[0])/2., coeffs2)/Gamma((1.0 + coeffs1[0])/2., coeffs2);
00087
00088 return (y*y)*pow((1. - (y/coeffs1[1])),coeffs1[0])*exp(-pow(cp,2)*pow((1.-(y/coeffs1[1])),2.));
00089
00090 }
00091
00092 double EvtBtoXsgammaFermiUtil::FermiGaussRootFcnB(double y, const std::vector<double> &coeffs1, const std::vector<double> &coeffs2) {
00093
00094
00095 double cp = Gamma((2.0 + coeffs1[0])/2., coeffs2)/Gamma((1.0 + coeffs1[0])/2., coeffs2);
00096 return pow((1. - (y/coeffs1[1])),coeffs1[0])*exp(-pow(cp,2)*pow((1.-(y/coeffs1[1])),2.));
00097
00098 }
00099
00100 double EvtBtoXsgammaFermiUtil::Gamma(double z, const std::vector<double> &coeffs) {
00101
00102
00103 double x, y, tmp, ser;
00104
00105 int j;
00106 y = z;
00107 x = z;
00108
00109 tmp = x + 5.5;
00110 tmp = tmp - (x+0.5)*log(tmp);
00111 ser=1.000000000190015;
00112
00113 for (j=0;j<6;j++) {
00114 y = y +1.0;
00115 ser = ser + coeffs[j]/y;
00116 }
00117
00118 return exp(-tmp+log(2.5066282746310005*ser/x));
00119
00120 }
00121
00122 double EvtBtoXsgammaFermiUtil::BesselK1(double x) {
00123
00124
00125
00126 if (x<0.0) report(INFO,"EvtGen") <<"x is negative !"<<endl;
00127
00128 double y, ans;
00129
00130 if (x <= 2.0) {
00131 y=x*x/4.0;
00132 ans = (log(x/2.0)*BesselI1(x))+(1.0/x)*(1.0+y*(0.15443144+y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1+y*(-0.110404e-2+y*(-0.4686e-4)))))));
00133 }
00134 else {
00135 y=2.0/x;
00136 ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619+y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2+y*(0.325614e-2+y*(-0.68245e-3)))))));
00137 }
00138 return ans;
00139
00140 }
00141
00142 double EvtBtoXsgammaFermiUtil::BesselI1(double x) {
00143
00144
00145
00146 double ax, ans;
00147 double y;
00148
00149 if ((ax=fabs(x)) < 3.75) {
00150 y=x/3.75;
00151 y*=y;
00152 ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934+y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
00153 }
00154 else {
00155 y=3.75/ax;
00156 ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1 -y*0.420059e-2));
00157 ans=0.398914228+y*(-0.3988024e-1+y*(-0.362018e-2+y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
00158 ans*=(exp(ax)/sqrt(ax));
00159 }
00160 return x < 0.0 ? -ans:ans;
00161 }
00162
00163 double EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot(double lambdabar, double lam1) {
00164
00165 EvtItgFunction *lhFunc = new EvtItgFunction(&FermiRomanRootFcnA, -1.e-6, 1.e6);
00166
00167 EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder();
00168 double rhSide = 1.0 - (lam1/(3.0*lambdabar*lambdabar));
00169
00170 double rho = rootFinder->GetRootSingleFunc(lhFunc, rhSide, 0.1, 0.4, 1.0e-6);
00171
00172 report(INFO,"EvtGen")<<"rho/2 "<<rho/2.<<" bessel "<<BesselK1(rho/2.)<<endl;
00173 double pF = lambdabar*sqrt(EvtConst::pi)/(rho*exp(rho/2.)*BesselK1(rho/2.));
00174 report(INFO,"EvtGen")<<"rho "<<rho<<" pf "<<pF<<endl;
00175
00176 delete lhFunc; lhFunc=0;
00177 delete rootFinder; rootFinder=0;
00178 return rho;
00179
00180 }
00181
00182 double EvtBtoXsgammaFermiUtil::FermiRomanRootFcnA(double y) {
00183
00184 return EvtConst::pi*(2. + y)*pow(y,-2.)*exp(-y)*pow(BesselK1(y/2.),-2.);
00185
00186 }
00187 double EvtBtoXsgammaFermiUtil::FermiRomanFunc(double y, const std::vector<double> &coeffs) {
00188 if (y == (coeffs[1]-coeffs[2])) y=0.99999999*(coeffs[1]-coeffs[2]);
00189
00190
00191 double pF = coeffs[4]*sqrt(EvtConst::pi)/(coeffs[3]*exp(coeffs[3]/2.)*BesselK1(coeffs[3]/2.));
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206 return (coeffs[1]-coeffs[2])*(1./(sqrt(EvtConst::pi)*pF))*exp(-(1./4.)*pow(pF*(coeffs[3]/((coeffs[1]-coeffs[2])*(1.-y/(coeffs[1]-coeffs[2])))) - ((coeffs[1]-coeffs[2])/pF)*(1. -y/(coeffs[1]-coeffs[2])),2.))/coeffs[5];
00207
00208
00209 }