Line data Source code
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(Severity::Info,"EvtGen")<<coeffs[4]<<endl;
44 0 : 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 0 : 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 0 : std::vector<double> coeffs1(3);
58 0 : std::vector<double> coeffs2(3);
59 :
60 0 : coeffs1[0]=0.2;
61 0 : coeffs1[1]=lambdabar;
62 0 : coeffs1[2]=0.0;
63 :
64 0 : coeffs2[0]=0.2;
65 0 : coeffs2[1]=lambdabar;
66 0 : coeffs2[2]=-lam1/3.;
67 :
68 0 : EvtItgTwoCoeffFcn *lhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnA, -mb, lambdabar, coeffs1, gammaCoeffs);
69 0 : EvtItgTwoCoeffFcn *rhFunc = new EvtItgTwoCoeffFcn(&FermiGaussRootFcnB, -mb, lambdabar, coeffs2, gammaCoeffs);
70 :
71 0 : EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder();
72 :
73 0 : 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 0 : delete rootFinder; rootFinder=0;
76 0 : delete lhFunc; lhFunc=0;
77 0 : delete rhFunc; rhFunc=0;
78 : return root;
79 :
80 0 : }
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 0 : double cp = Gamma((2.0 + coeffs1[0])/2., coeffs2)/Gamma((1.0 + coeffs1[0])/2., coeffs2);
87 :
88 0 : 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 0 : double cp = Gamma((2.0 + coeffs1[0])/2., coeffs2)/Gamma((1.0 + coeffs1[0])/2., coeffs2);
96 0 : 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 0 : tmp = x + 5.5;
110 0 : tmp = tmp - (x+0.5)*log(tmp);
111 : ser=1.000000000190015;
112 :
113 0 : for (j=0;j<6;j++) {
114 0 : y = y +1.0;
115 0 : ser = ser + coeffs[j]/y;
116 : }
117 :
118 0 : 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 0 : if (x<0.0) report(Severity::Info,"EvtGen") <<"x is negative !"<<endl;
127 :
128 : double y, ans;
129 :
130 0 : if (x <= 2.0) {
131 0 : y=x*x/4.0;
132 0 : 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 0 : }
134 : else {
135 0 : y=2.0/x;
136 0 : 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 0 : 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 0 : ax=fabs(x);
150 0 : if ( ax < 3.75) {
151 0 : y=x/3.75;
152 0 : y*=y;
153 0 : 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 0 : }
155 : else {
156 0 : y=3.75/ax;
157 0 : ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1 -y*0.420059e-2));
158 0 : ans=0.398914228+y*(-0.3988024e-1+y*(-0.362018e-2+y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
159 0 : ans*=(exp(ax)/sqrt(ax));
160 : }
161 0 : return x < 0.0 ? -ans:ans;
162 : }
163 :
164 : double EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot(double lambdabar, double lam1) {
165 :
166 0 : EvtItgFunction *lhFunc = new EvtItgFunction(&FermiRomanRootFcnA, -1.e-6, 1.e6);
167 :
168 0 : EvtBtoXsgammaRootFinder *rootFinder = new EvtBtoXsgammaRootFinder();
169 0 : double rhSide = 1.0 - (lam1/(3.0*lambdabar*lambdabar));
170 :
171 0 : double rho = rootFinder->GetRootSingleFunc(lhFunc, rhSide, 0.1, 0.4, 1.0e-6);
172 : //rho=0.250353;
173 0 : report(Severity::Info,"EvtGen")<<"rho/2 "<<rho/2.<<" bessel "<<BesselK1(rho/2.)<<endl;
174 0 : double pF = lambdabar*sqrt(EvtConst::pi)/(rho*exp(rho/2.)*BesselK1(rho/2.));
175 0 : report(Severity::Info,"EvtGen")<<"rho "<<rho<<" pf "<<pF<<endl;
176 :
177 0 : delete lhFunc; lhFunc=0;
178 0 : delete rootFinder; rootFinder=0;
179 0 : return rho;
180 :
181 0 : }
182 :
183 : double EvtBtoXsgammaFermiUtil::FermiRomanRootFcnA(double y) {
184 :
185 0 : 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 0 : 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 0 : double pF = coeffs[4]*sqrt(EvtConst::pi)/(coeffs[3]*exp(coeffs[3]/2.)*BesselK1(coeffs[3]/2.));
193 : // report(Severity::Info,"EvtGen")<<" pf "<<y<<" "<<pF<<" "<<coeffs[1]<<" "<<coeffs[2]<<" "<<coeffs[3]<<" "<<coeffs[4]<<" "<<coeffs[5]<<endl;
194 : //double pF=0.382533;
195 :
196 : //report(Severity::Info,"EvtGen")<<(coeffs[1]-coeffs[2])*(1./(sqrt(EvtConst::pi)*pF))<<endl;
197 : //report(Severity::Info,"EvtGen")<<(1.-y/(coeffs[1]-coeffs[2]))<<endl;
198 : //report(Severity::Info,"EvtGen")<<(coeffs[1]-coeffs[2])<<endl;
199 : //report(Severity::Info,"EvtGen")<<(coeffs[1]-coeffs[2])*(1.-y/(coeffs[1]-coeffs[2]))<<endl;
200 :
201 : //report(Severity::Info,"EvtGen")<<" "<<pF*coeffs[3]/((coeffs[1]-coeffs[2])*(1.-y/(coeffs[1]-coeffs[2])))<<endl;
202 : // report(Severity::Info,"EvtGen")<<" "<<((coeffs[1]-coeffs[2])/pF)*(1. -y/(coeffs[1]-coeffs[2]))<<endl;
203 :
204 : //report(Severity::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(Severity::Info,"EvtGen")<<"leaving"<<endl;
207 0 : 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 : }
|
