1 //--------------------------------------------------------------------------
4 // Copyright Information: See EvtGen/COPYRIGHT
7 // This software is part of the EvtGen package developed jointly
8 // for the BaBar and CLEO collaborations. If you use all or part
9 // of it, please give an appropriate acknowledgement.
11 // Module: EvtBtoXsgammaFermiUtil.cc
14 // Class to hold various fermi functions and their helper functions. The
15 // fermi functions are used in EvtBtoXsgammaKagan.
17 // Modification history:
19 // Jane Tinslay March 21, 2001 Module created
20 //------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
23 //-----------------------
24 // This Class's Header --
25 //-----------------------
26 #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
27 #include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
28 #include "EvtGenModels/EvtBtoXsgammaRootFinder.hh"
29 #include "EvtGenModels/EvtItgFunction.hh"
30 #include "EvtGenBase/EvtConst.hh"
31 #include "EvtGenBase/EvtReport.hh"
40 double EvtBtoXsgammaFermiUtil::FermiExpFunc(double y, const std::vector<double> &coeffs) {
42 //coeffs: 1 = lambdabar, 2 = a, 3 = lam1, 4 = norm
43 // report(INFO,"EvtGen")<<coeffs[4]<<endl;
44 return (pow(1. - (y/coeffs[1]),coeffs[2])*exp((-3.*pow(coeffs[1],2.)/coeffs[3])*y/coeffs[1]))/coeffs[4];
48 double EvtBtoXsgammaFermiUtil::FermiGaussFunc(double y, const std::vector<double> &coeffs) {
50 //coeffs: 1 = lambdabar, 2 = a, 3 = c, 4 = norm
51 return (pow(1. - (y/coeffs[1]),coeffs[2])*exp(-pow(coeffs[3],2.)*pow(1. - (y/coeffs[1]),2.)))/coeffs[4];
55 double EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot(double lambdabar, double lam1, double mb, std::vector<double> &gammaCoeffs) {
57 std::vector<double> coeffs1(3);
58 std::vector<double> coeffs2(3);
68 EvtItgTwoCoeffFcn *lhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnA, -mb, lambdabar, coeffs1, gammaCoeffs);
69 EvtItgTwoCoeffFcn *rhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnB, -mb, lambdabar, coeffs2, gammaCoeffs);
71 EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder();
73 double root = rootFinder->GetGaussIntegFcnRoot(lhFunc, rhFunc, 1.0e-4, 1.0e-4, 40, 40, -mb, lambdabar, 0.2, 0.4, 1.0e-6);
75 delete rootFinder; rootFinder=0;
76 delete lhFunc; lhFunc=0;
77 delete rhFunc; rhFunc=0;
82 double EvtBtoXsgammaFermiUtil::FermiGaussRootFcnA(double y, const std::vector<double> &coeffs1, const std::vector<double> &coeffs2) {
85 //coeffs1: 0=ap, 1=lambdabar, coeffs2=gamma function coeffs
86 double cp = Gamma((2.0 + coeffs1[0])/2., coeffs2)/Gamma((1.0 + coeffs1[0])/2., coeffs2);
88 return (y*y)*pow((1. - (y/coeffs1[1])),coeffs1[0])*exp(-pow(cp,2)*pow((1.-(y/coeffs1[1])),2.));
92 double EvtBtoXsgammaFermiUtil::FermiGaussRootFcnB(double y, const std::vector<double> &coeffs1, const std::vector<double> &coeffs2) {
94 //coeffs1: 0=ap, 1=lambdabar, coeffs2=gamma function coeffs
95 double cp = Gamma((2.0 + coeffs1[0])/2., coeffs2)/Gamma((1.0 + coeffs1[0])/2., coeffs2);
96 return pow((1. - (y/coeffs1[1])),coeffs1[0])*exp(-pow(cp,2)*pow((1.-(y/coeffs1[1])),2.));
100 double EvtBtoXsgammaFermiUtil::Gamma(double z, const std::vector<double> &coeffs) {
102 //Lifted from Numerical Recipies in C
103 double x, y, tmp, ser;
110 tmp = tmp - (x+0.5)*log(tmp);
111 ser=1.000000000190015;
115 ser = ser + coeffs[j]/y;
118 return exp(-tmp+log(2.5066282746310005*ser/x));
122 double EvtBtoXsgammaFermiUtil::BesselK1(double x) {
124 //Lifted from Numerical Recipies in C : Returns the modified Bessel
125 //function K_1(x) for positive real x
126 if (x<0.0) report(INFO,"EvtGen") <<"x is negative !"<<endl;
132 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)))))));
136 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)))))));
142 double EvtBtoXsgammaFermiUtil::BesselI1(double x) {
143 //Lifted from Numerical Recipies in C : Returns the modified Bessel
144 //function I_1(x) for any real x
153 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))))));
157 ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1 -y*0.420059e-2));
158 ans=0.398914228+y*(-0.3988024e-1+y*(-0.362018e-2+y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
159 ans*=(exp(ax)/sqrt(ax));
161 return x < 0.0 ? -ans:ans;
164 double EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot(double lambdabar, double lam1) {
166 EvtItgFunction *lhFunc = new EvtItgFunction(&FermiRomanRootFcnA, -1.e-6, 1.e6);
168 EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder();
169 double rhSide = 1.0 - (lam1/(3.0*lambdabar*lambdabar));
171 double rho = rootFinder->GetRootSingleFunc(lhFunc, rhSide, 0.1, 0.4, 1.0e-6);
173 report(INFO,"EvtGen")<<"rho/2 "<<rho/2.<<" bessel "<<BesselK1(rho/2.)<<endl;
174 double pF = lambdabar*sqrt(EvtConst::pi)/(rho*exp(rho/2.)*BesselK1(rho/2.));
175 report(INFO,"EvtGen")<<"rho "<<rho<<" pf "<<pF<<endl;
177 delete lhFunc; lhFunc=0;
178 delete rootFinder; rootFinder=0;
183 double EvtBtoXsgammaFermiUtil::FermiRomanRootFcnA(double y) {
185 return EvtConst::pi*(2. + y)*pow(y,-2.)*exp(-y)*pow(BesselK1(y/2.),-2.);
188 double EvtBtoXsgammaFermiUtil::FermiRomanFunc(double y, const std::vector<double> &coeffs) {
189 if (y == (coeffs[1]-coeffs[2])) y=0.99999999*(coeffs[1]-coeffs[2]);
191 //coeffs: 1 = mB, 2=mb, 3=rho, 4=lambdabar, 5=norm
192 double pF = coeffs[4]*sqrt(EvtConst::pi)/(coeffs[3]*exp(coeffs[3]/2.)*BesselK1(coeffs[3]/2.));
193 // report(INFO,"EvtGen")<<" pf "<<y<<" "<<pF<<" "<<coeffs[1]<<" "<<coeffs[2]<<" "<<coeffs[3]<<" "<<coeffs[4]<<" "<<coeffs[5]<<endl;
194 //double pF=0.382533;
196 //report(INFO,"EvtGen")<<(coeffs[1]-coeffs[2])*(1./(sqrt(EvtConst::pi)*pF))<<endl;
197 //report(INFO,"EvtGen")<<(1.-y/(coeffs[1]-coeffs[2]))<<endl;
198 //report(INFO,"EvtGen")<<(coeffs[1]-coeffs[2])<<endl;
199 //report(INFO,"EvtGen")<<(coeffs[1]-coeffs[2])*(1.-y/(coeffs[1]-coeffs[2]))<<endl;
201 //report(INFO,"EvtGen")<<" "<<pF*coeffs[3]/((coeffs[1]-coeffs[2])*(1.-y/(coeffs[1]-coeffs[2])))<<endl;
202 // report(INFO,"EvtGen")<<" "<<((coeffs[1]-coeffs[2])/pF)*(1. -y/(coeffs[1]-coeffs[2]))<<endl;
204 //report(INFO,"EvtGen")<<"result "<<(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];
206 //report(INFO,"EvtGen")<<"leaving"<<endl;
207 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];