Line data Source code
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 0 : {
53 0 : _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 0 : _epsilon1 = 1e-10;
68 0 : _epsilon2 = 1e-5;
69 0 : if ( alphas > 0 ) {
70 0 : double lne3 = 9./16.-2*EvtConst::pi*EvtConst::pi/3.+6*EvtConst::pi/4/alphas; if ( lne3 > 0 )
71 0 : lne3 = -7./4. - sqrt(lne3);
72 : else
73 : lne3 = -7./4.;
74 0 : _epsilon3 = exp(lne3);
75 0 : }
76 : else
77 0 : _epsilon3 = 1;
78 0 : }
79 :
80 : //--------------
81 : // Destructor --
82 : //--------------
83 :
84 : EvtVubdGamma::~EvtVubdGamma( )
85 0 : {
86 0 : }
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 0 : double xb = (1-x);
97 :
98 0 : if ( x < 0 || x > 1 || z < xb || z > (1+xb) )
99 0 : return 0;
100 :
101 0 : double p2min = (0>z-1.?0:z-1.);
102 0 : double p2max = (1.-x)*(z-1.+x);
103 :
104 0 : if (p2 < p2min || p2 > p2max)
105 0 : return 0;
106 :
107 : // // check the phase space
108 : // return 1.;
109 :
110 : double dG;
111 :
112 0 : if ( p2 >_epsilon1 && p2< _epsilon2) {
113 :
114 0 : double W1 = getW1delta(x,z);
115 0 : double W4plus5 = getW4plus5delta(x,z);
116 :
117 0 : dG = 12. * delta(p2,p2min,p2max) * ((1.+xb-z) * (z-xb) * W1
118 0 : + xb*(z-xb) * (W4plus5));
119 0 : }
120 : else {
121 :
122 0 : double W1 = getW1nodelta(x,z,p2);
123 0 : double W2 = getW2nodelta(x,z,p2);
124 0 : double W3 = getW3nodelta(x,z,p2);
125 0 : double W4 = getW4nodelta(x,z,p2);
126 0 : double W5 = getW5nodelta(x,z,p2);
127 :
128 0 : dG = 12. * ((1.+xb-z) * (z-xb-p2) * W1
129 0 : + (1.-z+p2) * W2
130 0 : + (xb*(z-xb)-p2) * (W3+W4+W5));
131 : }
132 : return dG;
133 0 : }
134 :
135 : double EvtVubdGamma::delta(const double &x, const double &xmin, const double &xmax)
136 : {
137 0 : if ( xmin > 0 || xmax < 0 ) return 0.;
138 0 : if ( _epsilon1 < x && x < _epsilon2 ) return 1./(_epsilon2-_epsilon1);
139 0 : return 0.0;
140 0 : }
141 :
142 : double EvtVubdGamma::getW1delta(const double &, const double &z)
143 : {
144 0 : double mz = 1.-z;
145 :
146 : double lz;
147 0 : if (z == 1) lz = -1.;
148 0 : 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 0 : double dl = 4.*EvtDiLog::DiLog(mz) + 4.*pow(EvtConst::pi,2)/3.;
157 :
158 0 : double w = -(8.*pow(log(z),2) - 10.*log(z) + 2.*lz + dl + 5.)
159 0 : + (8.*log(z)-7.)*log(_epsilon3) - 2.*pow(log(_epsilon3),2);
160 :
161 0 : return (1. + w*_alphas/3./EvtConst::pi);
162 : }
163 :
164 : double EvtVubdGamma::getW1nodelta(const double &, const double &z, const double &p2)
165 : {
166 :
167 0 : double z2 = z*z;
168 0 : double t2 = 1.-4.*p2/z2;
169 0 : double t = sqrt(t2);
170 :
171 : double w = 0;
172 0 : if ( p2 > _epsilon2 )
173 0 : w += 4./p2*(log((1.+t)/(1.-t))/t + log(p2/z2))
174 0 : + 1. - (8.-z)*(2.-z)/z2/t2
175 0 : + ((2.-z)/2./z+(8.-z)*(2.-z)/2./z2/t2)*log((1.+t)/(1.-t))/t;
176 0 : if ( p2 > _epsilon3 )
177 0 : w += (8.*log(z)-7.)/p2 - 4.*log(p2)/p2;
178 :
179 0 : return w*_alphas/3./EvtConst::pi;
180 : }
181 :
182 : double EvtVubdGamma::getW2nodelta(const double &, const double &z, const double &p2)
183 : {
184 :
185 0 : double z2 = z*z;
186 0 : double t2 = 1.-4.*p2/z2;
187 0 : double t = sqrt(t2);
188 0 : double w11 = (32.-8.*z+z2)/4./z/t2;
189 :
190 : double w = 0;
191 0 : if ( p2 > _epsilon2 )
192 0 : w -= (z*t2/8. + (4.-z)/4. + w11/2.)*log((1.+t)/(1.-t))/t;
193 0 : if ( p2 > _epsilon2 )
194 0 : w += (8.-z)/4. + w11;
195 :
196 0 : return (w*_alphas/3./EvtConst::pi);
197 : }
198 :
199 : double EvtVubdGamma::getW3nodelta(const double &, const double &z, const double &p2)
200 : {
201 0 : double z2 = z*z;
202 0 : double t2 = 1.-4.*p2/z2;
203 0 : double t4 = t2*t2;
204 0 : double t = sqrt(t2);
205 :
206 : double w = 0;
207 :
208 0 : if ( p2 > _epsilon2 )
209 0 : w += (z*t2/16. + 5.*(4.-z)/16. - (64.+56.*z-7.*z2)/16./z/t2
210 0 : + 3.*(12.-z)/16./t4) * log((1.+t)/(1.-t))/t;
211 0 : if ( p2 > _epsilon2 )
212 0 : w += -(8.-3.*z)/8. + (32.+22.*z-3.*z2)/4./z/t2 - 3.*(12.-z)/8./t4;
213 :
214 0 : return (w*_alphas/3./EvtConst::pi);
215 : }
216 :
217 : double EvtVubdGamma::getW4nodelta(const double &, const double &z, const double &p2)
218 : {
219 0 : double z2 = z*z;
220 0 : double t2 = 1.-4.*p2/z2;
221 0 : double t4 = t2*t2;
222 0 : double t = sqrt(t2);
223 :
224 : double w = 0;
225 :
226 0 : if ( p2 > _epsilon2 )
227 0 : w -= ((8.-3.*z)/4./z - (22.-3.*z)/2./z/t2 + 3.*(12.-z)/4./z/t4)
228 0 : * log((1.+t)/(1.-t))/t;
229 0 : if ( p2 > _epsilon2 )
230 0 : w += -1. - (32.-5.*z)/2./z/t2 + 3.*(12.-z)/2./z/t4 ;
231 :
232 0 : 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 0 : if ( z == 1 )
241 0 : w = -2;
242 : else
243 0 : w = 2.*log(z)/(1.-z);
244 :
245 0 : return (w*_alphas/3./EvtConst::pi);
246 : }
247 :
248 : double EvtVubdGamma::getW5nodelta(const double &, const double &z, const double &p2)
249 : {
250 0 : double z2 = z*z;
251 0 : double t2 = 1.-4.*p2/z2;
252 0 : double t4 = t2*t2;
253 0 : double t = sqrt(t2);
254 :
255 : double w = 0;
256 0 : if ( p2 > _epsilon2 )
257 0 : w += (1./4./z - (2.-z)/2./z2/t2 + 3.*(12.-z)/4./z2/t4)
258 0 : * log((1.+t)/(1.-t))/t;
259 0 : if ( p2 > _epsilon2 )
260 0 : w += -(8.+z)/2./z2/t2 - 3.*(12.-z)/2./z2/t4;
261 :
262 0 : return (w*_alphas/3./EvtConst::pi);
263 :
264 : }
265 :
266 :
267 :
|
