]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //----------------------------------------------------------------------- |
2 | // File and Version Information: | |
3 | // | |
4 | // Copyright Information: See EvtGen/COPYRIGHT | |
5 | // | |
6 | // | |
7 | // Description: | |
8 | // 3 2 2 | |
9 | // d Gamma / _ _ _2 mb _2 mb | |
10 | // ---------- = 12 Gamma | (1+x-z)(z-x-p ) -- W + (1-z+p ) -- W | |
11 | // _ 2 0 \ 2 1 2 2 | |
12 | // dx dz dp 2 | |
13 | // _ _ _2 mb 2 \. | |
14 | // + [x(z-x)-p ] -- (W + 2mb W + mb W ) | | |
15 | // 4 3 4 5 / | |
16 | // | |
17 | // with | |
18 | // 2 E 2 | |
19 | // l _2 p 2 v.p _ | |
20 | // x = ------ , p = --- , z = ------ , x = 1-x | |
21 | // mb 2 mb | |
22 | // mb | |
23 | // | |
24 | // the triple differential decay rate according to | |
25 | // hep-ph/9905351 v2 | |
26 | // | |
27 | // Environment: | |
28 | // Software developed for the BaBar Detector at the SLAC B-Factory. | |
29 | // | |
30 | // Author List: | |
31 | // Sven Menke | |
32 | // | |
33 | //----------------------------------------------------------------------- | |
34 | //----------------------- | |
35 | // This Class's Header -- | |
36 | //----------------------- | |
37 | #include "EvtGenBase/EvtPatches.hh" | |
38 | #include "EvtGenBase/EvtConst.hh" | |
39 | #include "EvtGenModels/EvtVubdGamma.hh" | |
40 | #include "EvtGenBase/EvtDiLog.hh" | |
41 | ||
42 | //--------------- | |
43 | // C Headers -- | |
44 | //--------------- | |
45 | #include <math.h> | |
46 | ||
47 | //---------------- | |
48 | // Constructors -- | |
49 | //---------------- | |
50 | ||
51 | EvtVubdGamma::EvtVubdGamma(const double &alphas) | |
52 | { | |
53 | _alphas = alphas; | |
54 | ||
55 | // the range for the delta distribution in p2 is from _epsilon1 to | |
56 | // _epsilon2. It was checked with the single differential formulae | |
57 | // in the paper that these values are small enough to imitate p2 = 0 | |
58 | // for the regular terms. | |
59 | // The ()* distributions, however need further treatment. In order to | |
60 | // generate the correct spectrum in z a threshold need to be computed | |
61 | // from the desired value of the coupling alphas. The idea is that | |
62 | // for z=1 p2=0 is not allowed and therefore the part of dGamma proportional | |
63 | // to delta(p2) should go to 0 for z->1. | |
64 | // Using equation (3.1) and (3.2) it is possible to find the correct value | |
65 | // for log(_epsilon3) from this requirement. | |
66 | ||
67 | _epsilon1 = 1e-10; | |
68 | _epsilon2 = 1e-5; | |
69 | if ( alphas > 0 ) { | |
70 | double lne3 = 9./16.-2*EvtConst::pi*EvtConst::pi/3.+6*EvtConst::pi/4/alphas; if ( lne3 > 0 ) | |
71 | lne3 = -7./4. - sqrt(lne3); | |
72 | else | |
73 | lne3 = -7./4.; | |
74 | _epsilon3 = exp(lne3); | |
75 | } | |
76 | else | |
77 | _epsilon3 = 1; | |
78 | } | |
79 | ||
80 | //-------------- | |
81 | // Destructor -- | |
82 | //-------------- | |
83 | ||
84 | EvtVubdGamma::~EvtVubdGamma( ) | |
85 | { | |
86 | } | |
87 | ||
88 | //----------- | |
89 | // Methods -- | |
90 | //----------- | |
91 | ||
92 | double EvtVubdGamma::getdGdxdzdp(const double &x, const double &z, const double &p2) | |
93 | { | |
94 | // check phase space | |
95 | ||
96 | double xb = (1-x); | |
97 | ||
98 | if ( x < 0 || x > 1 || z < xb || z > (1+xb) ) | |
99 | return 0; | |
100 | ||
101 | double p2min = (0>z-1.?0:z-1.); | |
102 | double p2max = (1.-x)*(z-1.+x); | |
103 | ||
104 | if (p2 < p2min || p2 > p2max) | |
105 | return 0; | |
106 | ||
107 | // // check the phase space | |
108 | // return 1.; | |
109 | ||
110 | double dG; | |
111 | ||
112 | if ( p2 >_epsilon1 && p2< _epsilon2) { | |
113 | ||
114 | double W1 = getW1delta(x,z); | |
115 | double W4plus5 = getW4plus5delta(x,z); | |
116 | ||
117 | dG = 12. * delta(p2,p2min,p2max) * ((1.+xb-z) * (z-xb) * W1 | |
118 | + xb*(z-xb) * (W4plus5)); | |
119 | } | |
120 | else { | |
121 | ||
122 | double W1 = getW1nodelta(x,z,p2); | |
123 | double W2 = getW2nodelta(x,z,p2); | |
124 | double W3 = getW3nodelta(x,z,p2); | |
125 | double W4 = getW4nodelta(x,z,p2); | |
126 | double W5 = getW5nodelta(x,z,p2); | |
127 | ||
128 | dG = 12. * ((1.+xb-z) * (z-xb-p2) * W1 | |
129 | + (1.-z+p2) * W2 | |
130 | + (xb*(z-xb)-p2) * (W3+W4+W5)); | |
131 | } | |
132 | return dG; | |
133 | } | |
134 | ||
135 | double EvtVubdGamma::delta(const double &x, const double &xmin, const double &xmax) | |
136 | { | |
137 | if ( xmin > 0 || xmax < 0 ) return 0.; | |
138 | if ( _epsilon1 < x && x < _epsilon2 ) return 1./(_epsilon2-_epsilon1); | |
139 | return 0.0; | |
140 | } | |
141 | ||
142 | double EvtVubdGamma::getW1delta(const double &, const double &z) | |
143 | { | |
144 | double mz = 1.-z; | |
145 | ||
146 | double lz; | |
147 | if (z == 1) lz = -1.; | |
148 | else lz = log(z)/(1.-z); | |
149 | ||
150 | // ddilog_(&z) is actually the dilog of (1-z) in maple, | |
151 | // also in Neuberts paper the limit dilog(1) = pi^2/6 is used | |
152 | // this corresponds to maple's dilog(0), so | |
153 | // I take ddilog_(&mz) where mz=1-z in order to satisfy Neubert's definition | |
154 | // and to compare with Maple the argument in maple should be (1-mz) ... | |
155 | ||
156 | double dl = 4.*EvtDiLog::DiLog(mz) + 4.*pow(EvtConst::pi,2)/3.; | |
157 | ||
158 | double w = -(8.*pow(log(z),2) - 10.*log(z) + 2.*lz + dl + 5.) | |
159 | + (8.*log(z)-7.)*log(_epsilon3) - 2.*pow(log(_epsilon3),2); | |
160 | ||
161 | return (1. + w*_alphas/3./EvtConst::pi); | |
162 | } | |
163 | ||
164 | double EvtVubdGamma::getW1nodelta(const double &, const double &z, const double &p2) | |
165 | { | |
166 | ||
167 | double z2 = z*z; | |
168 | double t2 = 1.-4.*p2/z2; | |
169 | double t = sqrt(t2); | |
170 | ||
171 | double w = 0; | |
172 | if ( p2 > _epsilon2 ) | |
173 | w += 4./p2*(log((1.+t)/(1.-t))/t + log(p2/z2)) | |
174 | + 1. - (8.-z)*(2.-z)/z2/t2 | |
175 | + ((2.-z)/2./z+(8.-z)*(2.-z)/2./z2/t2)*log((1.+t)/(1.-t))/t; | |
176 | if ( p2 > _epsilon3 ) | |
177 | w += (8.*log(z)-7.)/p2 - 4.*log(p2)/p2; | |
178 | ||
179 | return w*_alphas/3./EvtConst::pi; | |
180 | } | |
181 | ||
182 | double EvtVubdGamma::getW2nodelta(const double &, const double &z, const double &p2) | |
183 | { | |
184 | ||
185 | double z2 = z*z; | |
186 | double t2 = 1.-4.*p2/z2; | |
187 | double t = sqrt(t2); | |
188 | double w11 = (32.-8.*z+z2)/4./z/t2; | |
189 | ||
190 | double w = 0; | |
191 | if ( p2 > _epsilon2 ) | |
192 | w -= (z*t2/8. + (4.-z)/4. + w11/2.)*log((1.+t)/(1.-t))/t; | |
193 | if ( p2 > _epsilon2 ) | |
194 | w += (8.-z)/4. + w11; | |
195 | ||
196 | return (w*_alphas/3./EvtConst::pi); | |
197 | } | |
198 | ||
199 | double EvtVubdGamma::getW3nodelta(const double &, const double &z, const double &p2) | |
200 | { | |
201 | double z2 = z*z; | |
202 | double t2 = 1.-4.*p2/z2; | |
203 | double t4 = t2*t2; | |
204 | double t = sqrt(t2); | |
205 | ||
206 | double w = 0; | |
207 | ||
208 | if ( p2 > _epsilon2 ) | |
209 | w += (z*t2/16. + 5.*(4.-z)/16. - (64.+56.*z-7.*z2)/16./z/t2 | |
210 | + 3.*(12.-z)/16./t4) * log((1.+t)/(1.-t))/t; | |
211 | if ( p2 > _epsilon2 ) | |
212 | w += -(8.-3.*z)/8. + (32.+22.*z-3.*z2)/4./z/t2 - 3.*(12.-z)/8./t4; | |
213 | ||
214 | return (w*_alphas/3./EvtConst::pi); | |
215 | } | |
216 | ||
217 | double EvtVubdGamma::getW4nodelta(const double &, const double &z, const double &p2) | |
218 | { | |
219 | double z2 = z*z; | |
220 | double t2 = 1.-4.*p2/z2; | |
221 | double t4 = t2*t2; | |
222 | double t = sqrt(t2); | |
223 | ||
224 | double w = 0; | |
225 | ||
226 | if ( p2 > _epsilon2 ) | |
227 | w -= ((8.-3.*z)/4./z - (22.-3.*z)/2./z/t2 + 3.*(12.-z)/4./z/t4) | |
228 | * log((1.+t)/(1.-t))/t; | |
229 | if ( p2 > _epsilon2 ) | |
230 | w += -1. - (32.-5.*z)/2./z/t2 + 3.*(12.-z)/2./z/t4 ; | |
231 | ||
232 | return w*_alphas/3./EvtConst::pi; | |
233 | } | |
234 | ||
235 | double EvtVubdGamma::getW4plus5delta(const double &, const double &z) | |
236 | { | |
237 | ||
238 | double w = 0; | |
239 | ||
240 | if ( z == 1 ) | |
241 | w = -2; | |
242 | else | |
243 | w = 2.*log(z)/(1.-z); | |
244 | ||
245 | return (w*_alphas/3./EvtConst::pi); | |
246 | } | |
247 | ||
248 | double EvtVubdGamma::getW5nodelta(const double &, const double &z, const double &p2) | |
249 | { | |
250 | double z2 = z*z; | |
251 | double t2 = 1.-4.*p2/z2; | |
252 | double t4 = t2*t2; | |
253 | double t = sqrt(t2); | |
254 | ||
255 | double w = 0; | |
256 | if ( p2 > _epsilon2 ) | |
257 | w += (1./4./z - (2.-z)/2./z2/t2 + 3.*(12.-z)/4./z2/t4) | |
258 | * log((1.+t)/(1.-t))/t; | |
259 | if ( p2 > _epsilon2 ) | |
260 | w += -(8.+z)/2./z2/t2 - 3.*(12.-z)/2./z2/t4; | |
261 | ||
262 | return (w*_alphas/3./EvtConst::pi); | |
263 | ||
264 | } | |
265 | ||
266 | ||
267 |