]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //-------------------------------------------------------------------------- |
2 | // | |
3 | // | |
4 | // Copyright Information: See EvtGen/COPYRIGHT | |
5 | // | |
6 | // Environment: | |
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. | |
10 | // | |
11 | // Module: EvtBtoXsgammaFermiUtil.cc | |
12 | // | |
13 | // Description: | |
14 | // Class to hold various fermi functions and their helper functions. The | |
15 | // fermi functions are used in EvtBtoXsgammaKagan. | |
16 | // | |
17 | // Modification history: | |
18 | // | |
19 | // Jane Tinslay March 21, 2001 Module created | |
20 | //------------------------------------------------------------------------ | |
21 | #include "EvtGenBase/EvtPatches.hh" | |
22 | ||
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" | |
32 | ||
33 | //--------------- | |
34 | // C++ Headers -- | |
35 | //--------------- | |
36 | #include <iostream> | |
37 | #include <math.h> | |
38 | using std::endl; | |
39 | ||
40 | double EvtBtoXsgammaFermiUtil::FermiExpFunc(double y, const std::vector<double> &coeffs) { | |
41 | ||
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]; | |
45 | ||
46 | } | |
47 | ||
48 | double EvtBtoXsgammaFermiUtil::FermiGaussFunc(double y, const std::vector<double> &coeffs) { | |
49 | ||
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]; | |
52 | ||
53 | } | |
54 | ||
55 | double EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot(double lambdabar, double lam1, double mb, std::vector<double> &gammaCoeffs) { | |
56 | ||
57 | std::vector<double> coeffs1(3); | |
58 | std::vector<double> coeffs2(3); | |
59 | ||
60 | coeffs1[0]=0.2; | |
61 | coeffs1[1]=lambdabar; | |
62 | coeffs1[2]=0.0; | |
63 | ||
64 | coeffs2[0]=0.2; | |
65 | coeffs2[1]=lambdabar; | |
66 | coeffs2[2]=-lam1/3.; | |
67 | ||
68 | EvtItgTwoCoeffFcn *lhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnA, -mb, lambdabar, coeffs1, gammaCoeffs); | |
69 | EvtItgTwoCoeffFcn *rhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnB, -mb, lambdabar, coeffs2, gammaCoeffs); | |
70 | ||
71 | EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder(); | |
72 | ||
73 | double root = rootFinder->GetGaussIntegFcnRoot(lhFunc, rhFunc, 1.0e-4, 1.0e-4, 40, 40, -mb, lambdabar, 0.2, 0.4, 1.0e-6); | |
74 | ||
75 | delete rootFinder; rootFinder=0; | |
76 | delete lhFunc; lhFunc=0; | |
77 | delete rhFunc; rhFunc=0; | |
78 | return root; | |
79 | ||
80 | } | |
81 | ||
82 | double EvtBtoXsgammaFermiUtil::FermiGaussRootFcnA(double y, const std::vector<double> &coeffs1, const std::vector<double> &coeffs2) { | |
83 | ||
84 | ||
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); | |
87 | ||
88 | return (y*y)*pow((1. - (y/coeffs1[1])),coeffs1[0])*exp(-pow(cp,2)*pow((1.-(y/coeffs1[1])),2.)); | |
89 | ||
90 | } | |
91 | ||
92 | double EvtBtoXsgammaFermiUtil::FermiGaussRootFcnB(double y, const std::vector<double> &coeffs1, const std::vector<double> &coeffs2) { | |
93 | ||
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.)); | |
97 | ||
98 | } | |
99 | ||
100 | double EvtBtoXsgammaFermiUtil::Gamma(double z, const std::vector<double> &coeffs) { | |
101 | ||
102 | //Lifted from Numerical Recipies in C | |
103 | double x, y, tmp, ser; | |
104 | ||
105 | int j; | |
106 | y = z; | |
107 | x = z; | |
108 | ||
109 | tmp = x + 5.5; | |
110 | tmp = tmp - (x+0.5)*log(tmp); | |
111 | ser=1.000000000190015; | |
112 | ||
113 | for (j=0;j<6;j++) { | |
114 | y = y +1.0; | |
115 | ser = ser + coeffs[j]/y; | |
116 | } | |
117 | ||
118 | return exp(-tmp+log(2.5066282746310005*ser/x)); | |
119 | ||
120 | } | |
121 | ||
122 | double EvtBtoXsgammaFermiUtil::BesselK1(double x) { | |
123 | ||
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; | |
127 | ||
128 | double y, ans; | |
129 | ||
130 | if (x <= 2.0) { | |
131 | y=x*x/4.0; | |
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))))))); | |
133 | } | |
134 | else { | |
135 | y=2.0/x; | |
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))))))); | |
137 | } | |
138 | return ans; | |
139 | ||
140 | } | |
141 | ||
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 | |
145 | ||
146 | double ax, ans; | |
147 | double y; | |
148 | ||
149 | ax=fabs(x); | |
150 | if ( ax < 3.75) { | |
151 | y=x/3.75; | |
152 | y*=y; | |
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)))))); | |
154 | } | |
155 | else { | |
156 | y=3.75/ax; | |
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)); | |
160 | } | |
161 | return x < 0.0 ? -ans:ans; | |
162 | } | |
163 | ||
164 | double EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot(double lambdabar, double lam1) { | |
165 | ||
166 | EvtItgFunction *lhFunc = new EvtItgFunction(&FermiRomanRootFcnA, -1.e-6, 1.e6); | |
167 | ||
168 | EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder(); | |
169 | double rhSide = 1.0 - (lam1/(3.0*lambdabar*lambdabar)); | |
170 | ||
171 | double rho = rootFinder->GetRootSingleFunc(lhFunc, rhSide, 0.1, 0.4, 1.0e-6); | |
172 | //rho=0.250353; | |
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; | |
176 | ||
177 | delete lhFunc; lhFunc=0; | |
178 | delete rootFinder; rootFinder=0; | |
179 | return rho; | |
180 | ||
181 | } | |
182 | ||
183 | double EvtBtoXsgammaFermiUtil::FermiRomanRootFcnA(double y) { | |
184 | ||
185 | return EvtConst::pi*(2. + y)*pow(y,-2.)*exp(-y)*pow(BesselK1(y/2.),-2.); | |
186 | ||
187 | } | |
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]); | |
190 | ||
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; | |
195 | ||
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; | |
200 | ||
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; | |
203 | ||
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]; | |
205 | ||
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]; | |
208 | ||
209 | ||
210 | } |