EvtGen 2.2.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
Loading...
Searching...
No Matches
EvtBtoXsllUtil.cpp
Go to the documentation of this file.
1
2/***********************************************************************
3* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4* *
5* This file is part of EvtGen. *
6* *
7* EvtGen is free software: you can redistribute it and/or modify *
8* it under the terms of the GNU General Public License as published by *
9* the Free Software Foundation, either version 3 of the License, or *
10* (at your option) any later version. *
11* *
12* EvtGen is distributed in the hope that it will be useful, *
13* but WITHOUT ANY WARRANTY; without even the implied warranty of *
14* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15* GNU General Public License for more details. *
16* *
17* You should have received a copy of the GNU General Public License *
18* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19***********************************************************************/
20
21#include "EvtGenModels/EvtBtoXsllUtil.hh"
22
23#include "EvtGenBase/EvtComplex.hh"
24#include "EvtGenBase/EvtConst.hh"
25#include "EvtGenBase/EvtDiLog.hh"
26#include "EvtGenBase/EvtGenKine.hh"
27#include "EvtGenBase/EvtPDL.hh"
28#include "EvtGenBase/EvtParticle.hh"
29#include "EvtGenBase/EvtRandom.hh"
30#include "EvtGenBase/EvtReport.hh"
31
32#include <stdlib.h>
33
34EvtComplex EvtBtoXsllUtil::GetC7Eff0( double sh, bool nnlo )
35{
36 // This function returns the zeroth-order alpha_s part of C7
37
38 if ( !nnlo )
39 return -0.313;
40
41 double A7;
42
43 // use energy scale of 2.5 GeV as a computational trick (G.Hiller)
44 // at least for shat > 0.25
45 A7 = -0.353 + 0.023;
46
47 EvtComplex c7eff;
48 if ( sh > 0.25 ) {
49 c7eff = A7;
50 return c7eff;
51 }
52
53 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
54 A7 = -0.312 + 0.008;
55 c7eff = A7;
56
57 return c7eff;
58}
59
60EvtComplex EvtBtoXsllUtil::GetC7Eff1( double sh, double mbeff, bool nnlo )
61{
62 // This function returns the first-order alpha_s part of C7
63
64 if ( !nnlo )
65 return 0.0;
66 double logsh;
67 logsh = log( sh );
68
69 EvtComplex uniti( 0.0, 1.0 );
70
71 EvtComplex c7eff = 0.0;
72 if ( sh > 0.25 ) {
73 return c7eff;
74 }
75
76 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
77 double muscale = 5.0;
78 double alphas = 0.215;
79 //double A7 = -0.312 + 0.008;
80 double A8 = -0.148;
81 //double A9 = 4.174 + (-0.035);
82 //double A10 = -4.592 + 0.379;
83 double C1 = -0.487;
84 double C2 = 1.024;
85 //double T9 = 0.374 + 0.252;
86 //double U9 = 0.033 + 0.015;
87 //double W9 = 0.032 + 0.012;
88 double Lmu = log( muscale / mbeff );
89
90 EvtComplex F71;
91 EvtComplex f71;
92 EvtComplex k7100( -0.68192, -0.074998 );
93 EvtComplex k7101( 0.0, 0.0 );
94 EvtComplex k7110( -0.23935, -0.12289 );
95 EvtComplex k7111( 0.0027424, 0.019676 );
96 EvtComplex k7120( -0.0018555, -0.175 );
97 EvtComplex k7121( 0.022864, 0.011456 );
98 EvtComplex k7130( 0.28248, -0.12783 );
99 EvtComplex k7131( 0.029027, -0.0082265 );
100 f71 = k7100 + k7101 * logsh + sh * ( k7110 + k7111 * logsh ) +
101 sh * sh * ( k7120 + k7121 * logsh ) +
102 sh * sh * sh * ( k7130 + k7131 * logsh );
103 F71 = ( -208.0 / 243.0 ) * Lmu + f71;
104
105 EvtComplex F72;
106 EvtComplex f72;
107 EvtComplex k7200( 4.0915, 0.44999 );
108 EvtComplex k7201( 0.0, 0.0 );
109 EvtComplex k7210( 1.4361, 0.73732 );
110 EvtComplex k7211( -0.016454, -0.11806 );
111 EvtComplex k7220( 0.011133, 1.05 );
112 EvtComplex k7221( -0.13718, -0.068733 );
113 EvtComplex k7230( -1.6949, 0.76698 );
114 EvtComplex k7231( -0.17416, 0.049359 );
115 f72 = k7200 + k7201 * logsh + sh * ( k7210 + k7211 * logsh ) +
116 sh * sh * ( k7220 + k7221 * logsh ) +
117 sh * sh * sh * ( k7230 + k7231 * logsh );
118 F72 = ( 416.0 / 81.0 ) * Lmu + f72;
119
120 EvtComplex F78;
121 F78 = ( -32.0 / 9.0 ) * Lmu + 8.0 * EvtConst::pi * EvtConst::pi / 27.0 +
122 ( -44.0 / 9.0 ) + ( -8.0 * EvtConst::pi / 9.0 ) * uniti +
123 ( 4.0 / 3.0 * EvtConst::pi * EvtConst::pi - 40.0 / 3.0 ) * sh +
124 ( 32.0 * EvtConst::pi * EvtConst::pi / 9.0 - 316.0 / 9.0 ) * sh * sh +
125 ( 200.0 * EvtConst::pi * EvtConst::pi / 27.0 - 658.0 / 9.0 ) * sh *
126 sh * sh +
127 ( -8.0 * logsh / 9.0 ) * ( sh + sh * sh + sh * sh * sh );
128
129 c7eff = -alphas / ( 4.0 * EvtConst::pi ) * ( C1 * F71 + C2 * F72 + A8 * F78 );
130
131 return c7eff;
132}
133
134EvtComplex EvtBtoXsllUtil::GetC9Eff0( double sh, double /* mbeff */, bool nnlo,
135 bool btod )
136{
137 // This function returns the zeroth-order alpha_s part of C9
138
139 if ( !nnlo )
140 return 4.344;
141 double mch = 0.29;
142
143 double A9;
144 A9 = 4.287 + ( -0.218 );
145 double C1;
146 C1 = -0.697;
147 double C2;
148 C2 = 1.046;
149 double T9;
150 T9 = 0.114 + 0.280;
151 double U9;
152 U9 = 0.045 + 0.023;
153 double W9;
154 W9 = 0.044 + 0.016;
155
156 EvtComplex uniti( 0.0, 1.0 );
157
158 EvtComplex hc;
159 double xarg;
160 xarg = 4.0 * mch / sh;
161
162 hc = -4.0 / 9.0 * log( mch * mch ) + 8.0 / 27.0 + 4.0 * xarg / 9.0;
163 if ( xarg < 1.0 ) {
164 hc = hc - 2.0 / 9.0 * ( 2.0 + xarg ) * sqrt( fabs( 1.0 - xarg ) ) *
165 ( log( ( sqrt( 1.0 - xarg ) + 1.0 ) /
166 ( sqrt( 1.0 - xarg ) - 1.0 ) ) -
167 uniti * EvtConst::pi );
168 } else {
169 hc = hc - 2.0 / 9.0 * ( 2.0 + xarg ) * sqrt( fabs( 1.0 - xarg ) ) *
170 2.0 * atan( 1.0 / sqrt( xarg - 1.0 ) );
171 }
172
173 EvtComplex h1;
174 xarg = 4.0 / sh;
175 h1 = 8.0 / 27.0 + 4.0 * xarg / 9.0;
176 if ( xarg < 1.0 ) {
177 h1 = h1 - 2.0 / 9.0 * ( 2.0 + xarg ) * sqrt( fabs( 1.0 - xarg ) ) *
178 ( log( ( sqrt( 1.0 - xarg ) + 1.0 ) /
179 ( sqrt( 1.0 - xarg ) - 1.0 ) ) -
180 uniti * EvtConst::pi );
181 } else {
182 h1 = h1 - 2.0 / 9.0 * ( 2.0 + xarg ) * sqrt( fabs( 1.0 - xarg ) ) *
183 2.0 * atan( 1.0 / sqrt( xarg - 1.0 ) );
184 }
185
186 EvtComplex h0;
187 h0 = 8.0 / 27.0 - 4.0 * log( 2.0 ) / 9.0 + 4.0 * uniti * EvtConst::pi / 9.0;
188
189 // X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
190 // h(\hat m_u^2, hat s))
191 EvtComplex Vudstar( 1.0 - 0.2279 * 0.2279 / 2.0, 0.0 );
192 EvtComplex Vub( ( 0.118 + 0.273 ) / 2.0, -1.0 * ( 0.305 + 0.393 ) / 2.0 );
193 EvtComplex Vtdstar( 1.0 - ( 0.118 + 0.273 ) / 2.0, ( 0.305 + 0.393 ) / 2.0 );
194 EvtComplex Vtb( 1.0, 0.0 );
195
196 EvtComplex Xd;
197 Xd = ( Vudstar * Vub / Vtdstar * Vtb ) * ( 4.0 / 3.0 * C1 + C2 ) *
198 ( hc - h0 );
199
200 EvtComplex c9eff = 4.344;
201 if ( sh > 0.25 ) {
202 c9eff = A9 + T9 * hc + U9 * h1 + W9 * h0;
203 if ( btod ) {
204 c9eff += Xd;
205 }
206 return c9eff;
207 }
208
209 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
210 A9 = 4.174 + ( -0.035 );
211 C1 = -0.487;
212 C2 = 1.024;
213 T9 = 0.374 + 0.252;
214 U9 = 0.033 + 0.015;
215 W9 = 0.032 + 0.012;
216
217 Xd = ( Vudstar * Vub / Vtdstar * Vtb ) * ( 4.0 / 3.0 * C1 + C2 ) *
218 ( hc - h0 );
219
220 c9eff = A9 + T9 * hc + U9 * h1 + W9 * h0;
221
222 if ( btod ) {
223 c9eff += Xd;
224 }
225
226 return c9eff;
227}
228
229EvtComplex EvtBtoXsllUtil::GetC9Eff1( double sh, double mbeff, bool nnlo,
230 bool /*btod*/ )
231{
232 // This function returns the first-order alpha_s part of C9
233
234 if ( !nnlo )
235 return 0.0;
236 double logsh;
237 logsh = log( sh );
238 double mch = 0.29;
239
240 EvtComplex uniti( 0.0, 1.0 );
241
242 EvtComplex c9eff = 0.0;
243 if ( sh > 0.25 ) {
244 return c9eff;
245 }
246
247 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
248 double muscale = 5.0;
249 double alphas = 0.215;
250 double C1 = -0.487;
251 double C2 = 1.024;
252 double A8 = -0.148;
253 double Lmu = log( muscale / mbeff );
254
255 EvtComplex F91;
256 EvtComplex f91;
257 EvtComplex k9100( -11.973, 0.16371 );
258 EvtComplex k9101( -0.081271, -0.059691 );
259 EvtComplex k9110( -28.432, -0.25044 );
260 EvtComplex k9111( -0.040243, 0.016442 );
261 EvtComplex k9120( -57.114, -0.86486 );
262 EvtComplex k9121( -0.035191, 0.027909 );
263 EvtComplex k9130( -128.8, -2.5243 );
264 EvtComplex k9131( -0.017587, 0.050639 );
265 f91 = k9100 + k9101 * logsh + sh * ( k9110 + k9111 * logsh ) +
266 sh * sh * ( k9120 + k9121 * logsh ) +
267 sh * sh * sh * ( k9130 + k9131 * logsh );
268 F91 = ( -1424.0 / 729.0 + 16.0 * uniti * EvtConst::pi / 243.0 +
269 64.0 / 27.0 * log( mch ) ) *
270 Lmu -
271 16.0 * Lmu * logsh / 243.0 +
272 ( 16.0 / 1215.0 - 32.0 / 135.0 / mch / mch ) * Lmu * sh +
273 ( 4.0 / 2835.0 - 8.0 / 315.0 / mch / mch / mch / mch ) * Lmu * sh * sh +
274 ( 16.0 / 76545.0 - 32.0 / 8505.0 / mch / mch / mch / mch / mch / mch ) *
275 Lmu * sh * sh * sh -
276 256.0 * Lmu * Lmu / 243.0 + f91;
277
278 EvtComplex F92;
279 EvtComplex f92;
280 EvtComplex k9200( 6.6338, -0.98225 );
281 EvtComplex k9201( 0.48763, 0.35815 );
282 EvtComplex k9210( 3.3585, 1.5026 );
283 EvtComplex k9211( 0.24146, -0.098649 );
284 EvtComplex k9220( -1.1906, 5.1892 );
285 EvtComplex k9221( 0.21115, -0.16745 );
286 EvtComplex k9230( -17.12, 15.146 );
287 EvtComplex k9231( 0.10552, -0.30383 );
288 f92 = k9200 + k9201 * logsh + sh * ( k9210 + k9211 * logsh ) +
289 sh * sh * ( k9220 + k9221 * logsh ) +
290 sh * sh * sh * ( k9230 + k9231 * logsh );
291 F92 = ( 256.0 / 243.0 - 32.0 * uniti * EvtConst::pi / 81.0 -
292 128.0 / 9.0 * log( mch ) ) *
293 Lmu +
294 32.0 * Lmu * logsh / 81.0 +
295 ( -32.0 / 405.0 + 64.0 / 45.0 / mch / mch ) * Lmu * sh +
296 ( -8.0 / 945.0 + 16.0 / 105.0 / mch / mch / mch / mch ) * Lmu * sh * sh +
297 ( -32.0 / 25515.0 + 64.0 / 2835.0 / mch / mch / mch / mch / mch / mch ) *
298 Lmu * sh * sh * sh +
299 512.0 * Lmu * Lmu / 81.0 + f92;
300
301 EvtComplex F98;
302 F98 = 104.0 / 9.0 - 32.0 * EvtConst::pi * EvtConst::pi / 27.0 +
303 ( 1184.0 / 27.0 - 40.0 * EvtConst::pi * EvtConst::pi / 9.0 ) * sh +
304 ( 14212.0 / 135.0 - 32.0 * EvtConst::pi * EvtConst::pi / 3.0 ) * sh *
305 sh +
306 ( 193444.0 / 945.0 - 560.0 * EvtConst::pi * EvtConst::pi / 27.0 ) *
307 sh * sh * sh +
308 16.0 * logsh / 9.0 * ( 1.0 + sh + sh * sh + sh * sh * sh );
309
310 c9eff = -alphas / ( 4.0 * EvtConst::pi ) * ( C1 * F91 + C2 * F92 + A8 * F98 );
311
312 return c9eff;
313}
314
315EvtComplex EvtBtoXsllUtil::GetC10Eff( double /*sh*/, bool nnlo )
316{
317 if ( !nnlo )
318 return -4.669;
319 double A10;
320 A10 = -4.592 + 0.379;
321
322 EvtComplex c10eff;
323 c10eff = A10;
324
325 return c10eff;
326}
327
328double EvtBtoXsllUtil::dGdsProb( double mb, double ms, double ml, double s )
329{
330 // Compute the decay probability density function given a value of s
331 // according to Ali-Lunghi-Greub-Hiller's 2002 paper
332 // Note that the form given below is taken from
333 // F.Kruger and L.M.Sehgal, Phys. Lett. B380, 199 (1996)
334 // but the differential rate as a function of dilepton mass
335 // in this latter paper reduces to Eq.(12) in ALGH's 2002 paper
336 // for ml = 0 and ms = 0.
337
338 bool btod = false;
339 bool nnlo = true;
340
341 double delta, lambda, prob;
342 double f1, f2, f3, f4;
343 double msh, mlh, sh;
344 double mbeff = 4.8;
345
346 mlh = ml / mb;
347 msh = ms / mb;
348 // set lepton and strange-quark masses to 0 if need to
349 // be in strict agreement with ALGH 2002 paper
350 // mlh = 0.0; msh = 0.0;
351 // sh = s / (mb*mb);
352 sh = s / ( mbeff * mbeff );
353
354 // if sh >1.0 code will return a nan. so just skip it
355 if ( sh > 1.0 )
356 return 0.0;
357
358 EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0( sh, nnlo );
359 EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1( sh, mbeff, nnlo );
360 EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0( sh, mbeff, nnlo, btod );
361 EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1( sh, mbeff, nnlo, btod );
362 EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff( sh, nnlo );
363
364 double alphas = 0.119 /
365 ( 1 + 0.119 * log( pow( 4.8, 2 ) / pow( 91.1867, 2 ) ) *
366 23.0 / 12.0 / EvtConst::pi );
367
368 double omega7 = -8.0 / 3.0 * log( 4.8 / mb ) -
369 4.0 / 3.0 * EvtDiLog::DiLog( sh ) -
370 2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
371 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
372 log( 1 - sh ) * ( 8.0 + sh ) / ( 2.0 + sh ) / 3.0 -
373 2.0 / 3.0 * sh * ( 2.0 - 2.0 * sh - sh * sh ) * log( sh ) /
374 pow( ( 1.0 - sh ), 2 ) / ( 2.0 + sh ) -
375 ( 16.0 - 11.0 * sh - 17.0 * sh * sh ) / 18.0 /
376 ( 2.0 + sh ) / ( 1.0 - sh );
377 double eta7 = 1.0 + alphas * omega7 / EvtConst::pi;
378
379 double omega79 = -4.0 / 3.0 * log( 4.8 / mb ) -
380 4.0 / 3.0 * EvtDiLog::DiLog( sh ) -
381 2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
382 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
383 1.0 / 9.0 * ( 2.0 + 7.0 * sh ) * log( 1.0 - sh ) / sh -
384 2.0 / 9.0 * sh * ( 3.0 - 2.0 * sh ) * log( sh ) /
385 pow( ( 1.0 - sh ), 2 ) +
386 1.0 / 18.0 * ( 5.0 - 9.0 * sh ) / ( 1.0 - sh );
387 double eta79 = 1.0 + alphas * omega79 / EvtConst::pi;
388
389 double omega9 = -2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
390 4.0 / 3.0 * EvtDiLog::DiLog( sh ) -
391 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
392 ( 5.0 + 4.0 * sh ) / ( 3.0 * ( 1.0 + 2.0 * sh ) ) *
393 log( 1.0 - sh ) -
394 2.0 * sh * ( 1.0 + sh ) * ( 1.0 - 2.0 * sh ) /
395 ( 3.0 * pow( 1.0 - sh, 2 ) * ( 1.0 + 2.0 * sh ) ) *
396 log( sh ) +
397 ( 5.0 + 9.0 * sh - 6.0 * sh * sh ) /
398 ( 6.0 * ( 1.0 - sh ) * ( 1.0 + 2.0 * sh ) );
399 double eta9 = 1.0 + alphas * omega9 / EvtConst::pi;
400
401 EvtComplex c7eff = eta7 * c7eff0 + c7eff1;
402 EvtComplex c9eff = eta9 * c9eff0 + c9eff1;
403 c10eff *= eta9;
404
405 double c7c7 = abs2( c7eff );
406 double c7c9 = real( ( eta79 * c7eff0 + c7eff1 ) *
407 conj( eta79 * c9eff0 + c9eff1 ) );
408 double c9c9plusc10c10 = abs2( c9eff ) + abs2( c10eff );
409 double c9c9minusc10c10 = abs2( c9eff ) - abs2( c10eff );
410
411 // Power corrections according to ALGH 2002
412 double lambda_1 = -0.2;
413 double lambda_2 = 0.12;
414 double C1 = -0.487;
415 double C2 = 1.024;
416 double mc = 0.29 * mb;
417
418 EvtComplex F;
419 double r = s / ( 4.0 * mc * mc );
420 EvtComplex uniti( 0.0, 1.0 );
421 F = 3.0 / ( 2.0 * r );
422 if ( r < 1 ) {
423 F *= 1.0 / sqrt( r * ( 1.0 - r ) ) * atan( sqrt( r / ( 1.0 - r ) ) ) -
424 1.0;
425 } else {
426 F *= 0.5 / sqrt( r * ( r - 1.0 ) ) *
427 ( log( ( 1.0 - sqrt( 1.0 - 1.0 / r ) ) /
428 ( 1.0 + sqrt( 1.0 - 1.0 / r ) ) ) +
429 uniti * EvtConst::pi ) -
430 1.0;
431 }
432
433 double G1 = 1.0 + lambda_1 / ( 2.0 * mb * mb ) +
434 3.0 * ( 1.0 - 15.0 * sh * sh + 10.0 * sh * sh * sh ) /
435 ( ( 1.0 - sh ) * ( 1.0 - sh ) * ( 1.0 + 2.0 * sh ) ) *
436 lambda_2 / ( 2.0 * mb * mb );
437 double G2 = 1.0 + lambda_1 / ( 2.0 * mb * mb ) -
438 3.0 * ( 6.0 + 3.0 * sh - 5.0 * sh * sh * sh ) /
439 ( ( 1.0 - sh ) * ( 1.0 - sh ) * ( 2.0 + sh ) ) * lambda_2 /
440 ( 2.0 * mb * mb );
441 double G3 = 1.0 + lambda_1 / ( 2.0 * mb * mb ) -
442 ( 5.0 + 6.0 * sh - 7.0 * sh * sh ) /
443 ( ( 1.0 - sh ) * ( 1.0 - sh ) ) * lambda_2 /
444 ( 2.0 * mb * mb );
445 double Gc = -8.0 / 9.0 * ( C2 - C1 / 6.0 ) * lambda_2 / ( mc * mc ) *
446 real( F * ( conj( c9eff ) * ( 2.0 + sh ) +
447 conj( c7eff ) * ( 1.0 + 6.0 * sh - sh * sh ) / sh ) );
448
449 // end of power corrections section
450 // now back to Kruger & Sehgal expressions
451
452 double msh2 = msh * msh;
453 lambda = 1.0 + sh * sh + msh2 * msh2 - 2.0 * ( sh + sh * msh2 + msh2 );
454 // negative lambda screw up sqrt below!
455 if ( lambda < 0.0 )
456 return 0.0;
457
458 f1 = pow( 1.0 - msh2, 2 ) - sh * ( 1.0 + msh2 );
459 f2 = 2.0 * ( 1.0 + msh2 ) * pow( 1.0 - msh2, 2 ) -
460 sh * ( 1.0 + 14.0 * msh2 + pow( msh, 4 ) ) - sh * sh * ( 1.0 + msh2 );
461 f3 = pow( 1.0 - msh2, 2 ) + sh * ( 1.0 + msh2 ) - 2.0 * sh * sh +
462 lambda * 2.0 * mlh * mlh / sh;
463 f4 = 1.0 - sh + msh2;
464
465 delta = ( 12.0 * c7c9 * f1 * G3 + 4.0 * c7c7 * f2 * G2 / sh ) *
466 ( 1.0 + 2.0 * mlh * mlh / sh ) +
467 c9c9plusc10c10 * f3 * G1 + 6.0 * mlh * mlh * c9c9minusc10c10 * f4 +
468 Gc;
469
470 // avoid negative probs
471 if ( delta < 0.0 )
472 delta = 0.;
473 // negative when sh < 4*mlh*mlh
474 // s < 4*ml*ml
476 prob = sqrt( lambda * ( 1.0 - 4.0 * ml * ml / s ) ) * delta;
477
478 // if ( !(prob>=0.0) && !(prob<=0.0) ) {
479 //nan
480 // std::cout << lambda << " " << mlh << " " << sh << " " << delta << " " << mb << " " << mbeff << std::endl;
481 // std::cout << 4.0*mlh*mlh/sh << " " << 4.0*ml*ml/s << " " << s-4.0*ml*ml << " " << ml << std::endl;
482 // std::cout << sh << " " << sh*sh << " " << msh2*msh2 << " " << msh << std::endl;
483 //std::cout << ( 12.0*c7c9*f1*G3 + 4.0*c7c7*f2*G2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
484 // <<" " << c9c9plusc10c10*f3*G1
485 // << " "<< 6.0*mlh*mlh*c9c9minusc10c10*f4
486 // << " "<< Gc << std::endl;
487 //std::cout << C2 << " " << C1 << " "<< lambda_2 << " " << mc << " " << real(F*(conj(c9eff)*(2.0+sh)+conj(c7eff)*(1.0 + 6.0*sh - sh*sh)/sh)) << " " << sh << " " << r << std::endl;
488 //std::cout << c9eff << " " << eta9 << " " <<c9eff0 << " " << c9eff1 << " " << alphas << " " << omega9 << " " << sh << std::endl;
489
490 //}
491 // else{
492 // if ( sh > 1.0) std::cout << "not a nan \n";
493 // }
494 return prob;
495}
496
497double EvtBtoXsllUtil::dGdsdupProb( double mb, double ms, double ml, double s,
498 double u )
499{
500 // Compute the decay probability density function given a value of s and u
501 // according to Ali-Hiller-Handoko-Morozumi's 1997 paper
502 // see Appendix E
503
504 bool btod = false;
505 bool nnlo = true;
506
507 double prob;
508 double f1sp, f2sp, f3sp;
509 double mbeff = 4.8;
510
511 // double sh = s / (mb*mb);
512 double sh = s / ( mbeff * mbeff );
513
514 // if sh >1.0 code will return a nan. so just skip it
515 if ( sh > 1.0 )
516 return 0.0;
517
518 EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0( sh, nnlo );
519 EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1( sh, mbeff, nnlo );
520 EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0( sh, mbeff, nnlo, btod );
521 EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1( sh, mbeff, nnlo, btod );
522 EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff( sh, nnlo );
523
524 double alphas = 0.119 /
525 ( 1 + 0.119 * log( pow( 4.8, 2 ) / pow( 91.1867, 2 ) ) *
526 23.0 / 12.0 / EvtConst::pi );
527
528 double omega7 = -8.0 / 3.0 * log( 4.8 / mb ) -
529 4.0 / 3.0 * EvtDiLog::DiLog( sh ) -
530 2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
531 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
532 log( 1 - sh ) * ( 8.0 + sh ) / ( 2.0 + sh ) / 3.0 -
533 2.0 / 3.0 * sh * ( 2.0 - 2.0 * sh - sh * sh ) * log( sh ) /
534 pow( ( 1.0 - sh ), 2 ) / ( 2.0 + sh ) -
535 ( 16.0 - 11.0 * sh - 17.0 * sh * sh ) / 18.0 /
536 ( 2.0 + sh ) / ( 1.0 - sh );
537 double eta7 = 1.0 + alphas * omega7 / EvtConst::pi;
538
539 double omega79 = -4.0 / 3.0 * log( 4.8 / mb ) -
540 4.0 / 3.0 * EvtDiLog::DiLog( sh ) -
541 2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
542 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
543 1.0 / 9.0 * ( 2.0 + 7.0 * sh ) * log( 1.0 - sh ) / sh -
544 2.0 / 9.0 * sh * ( 3.0 - 2.0 * sh ) * log( sh ) /
545 pow( ( 1.0 - sh ), 2 ) +
546 1.0 / 18.0 * ( 5.0 - 9.0 * sh ) / ( 1.0 - sh );
547 double eta79 = 1.0 + alphas * omega79 / EvtConst::pi;
548
549 double omega9 = -2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
550 4.0 / 3.0 * EvtDiLog::DiLog( sh ) -
551 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
552 ( 5.0 + 4.0 * sh ) / ( 3.0 * ( 1.0 + 2.0 * sh ) ) *
553 log( 1.0 - sh ) -
554 2.0 * sh * ( 1.0 + sh ) * ( 1.0 - 2.0 * sh ) /
555 ( 3.0 * pow( 1.0 - sh, 2 ) * ( 1.0 + 2.0 * sh ) ) *
556 log( sh ) +
557 ( 5.0 + 9.0 * sh - 6.0 * sh * sh ) /
558 ( 6.0 * ( 1.0 - sh ) * ( 1.0 + 2.0 * sh ) );
559 double eta9 = 1.0 + alphas * omega9 / EvtConst::pi;
560
561 EvtComplex c7eff = eta7 * c7eff0 + c7eff1;
562 EvtComplex c9eff = eta9 * c9eff0 + c9eff1;
563 c10eff *= eta9;
564
565 double c7c7 = abs2( c7eff );
566 double c7c9 = real( ( eta79 * c7eff0 + c7eff1 ) *
567 conj( eta79 * c9eff0 + c9eff1 ) );
568 double c7c10 = real( ( eta79 * c7eff0 + c7eff1 ) * conj( eta9 * c10eff ) );
569 double c9c10 = real( ( eta9 * c9eff0 + c9eff1 ) * conj( eta9 * c10eff ) );
570 double c9c9plusc10c10 = abs2( c9eff ) + abs2( c10eff );
571
572 f1sp = ( pow( mb * mb - ms * ms, 2 ) - s * s ) * c9c9plusc10c10 +
573 4.0 *
574 ( pow( mb, 4 ) - ms * ms * mb * mb -
575 pow( ms, 4 ) * ( 1.0 - ms * ms / ( mb * mb ) ) -
576 8.0 * s * ms * ms - s * s * ( 1.0 + ms * ms / ( mb * mb ) ) ) *
577 mb * mb * c7c7 /
578 s
579 // kludged mass term
580 * ( 1.0 + 2.0 * ml * ml / s ) -
581 8.0 * ( s * ( mb * mb + ms * ms ) - pow( mb * mb - ms * ms, 2 ) ) *
582 c7c9
583 // kludged mass term
584 * ( 1.0 + 2.0 * ml * ml / s );
585
586 f2sp = 4.0 * s * c9c10 + 8.0 * ( mb * mb + ms * ms ) * c7c10;
587 f3sp = -( c9c9plusc10c10 ) + 4.0 * ( 1.0 + pow( ms / mb, 4 ) ) * mb * mb *
588 c7c7 /
589 s
590 // kludged mass term
591 * ( 1.0 + 2.0 * ml * ml / s );
592
593 prob = ( f1sp + f2sp * u + f3sp * u * u ) / pow( mb, 3 );
594 if ( prob < 0.0 )
595 prob = 0.;
596
597 return prob;
598}
599
600double EvtBtoXsllUtil::FermiMomentum( double pf )
601{
602 // Pick a value for the b-quark Fermi motion momentum
603 // according to Ali's Gaussian model
604
605 double pb, pbmax, xbox, ybox;
606 pb = 0.0;
607 pbmax = 5.0 * pf;
608
609 while ( pb == 0.0 ) {
610 xbox = EvtRandom::Flat( pbmax );
611 ybox = EvtRandom::Flat();
612 if ( ybox < FermiMomentumProb( xbox, pf ) ) {
613 pb = xbox;
614 }
615 }
616
617 return pb;
618}
619
620double EvtBtoXsllUtil::FermiMomentumProb( double pb, double pf )
621{
622 // Compute probability according to Ali's Gaussian model
623 // the function chosen has a convenient maximum value of 1 for pb = pf
624
625 double prsq = ( pb * pb ) / ( pf * pf );
626 double prob = prsq * exp( 1.0 - prsq );
627
628 return prob;
629}
conj
EvtComplex conj(const EvtComplex &c)
Definition EvtComplex.hh:198
real
double real(const EvtComplex &c)
Definition EvtComplex.hh:232
abs2
double abs2(const EvtComplex &c)
Definition EvtComplex.hh:211
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:242
lambda
double lambda(double q, double m1, double m2)
Definition EvtFlatQ2.cpp:30
EvtBtoXsllUtil::GetC10Eff
EvtComplex GetC10Eff(double sh, bool nnlo=true)
Definition EvtBtoXsllUtil.cpp:315
EvtBtoXsllUtil::dGdsProb
double dGdsProb(double mb, double ms, double ml, double s)
Definition EvtBtoXsllUtil.cpp:328
EvtBtoXsllUtil::FermiMomentumProb
double FermiMomentumProb(double pb, double pf)
Definition EvtBtoXsllUtil.cpp:620
EvtBtoXsllUtil::FermiMomentum
double FermiMomentum(double pf)
Definition EvtBtoXsllUtil.cpp:600
EvtBtoXsllUtil::GetC9Eff0
EvtComplex GetC9Eff0(double sh, double mb, bool nnlo=true, bool btod=false)
Definition EvtBtoXsllUtil.cpp:134
EvtBtoXsllUtil::GetC7Eff1
EvtComplex GetC7Eff1(double sh, double mb, bool nnlo=true)
Definition EvtBtoXsllUtil.cpp:60
EvtBtoXsllUtil::GetC7Eff0
EvtComplex GetC7Eff0(double sh, bool nnlo=true)
Definition EvtBtoXsllUtil.cpp:34
EvtBtoXsllUtil::GetC9Eff1
EvtComplex GetC9Eff1(double sh, double mb, bool nnlo=true, bool btod=false)
Definition EvtBtoXsllUtil.cpp:229
EvtBtoXsllUtil::dGdsdupProb
double dGdsdupProb(double mb, double ms, double ml, double s, double u)
Definition EvtBtoXsllUtil.cpp:497
EvtConst::pi
static const double pi
Definition EvtConst.hh:26
EvtRandom::Flat
static double Flat()
Definition EvtRandom.cpp:95
EvtDiLog::DiLog
double DiLog(double x)
Definition EvtDiLog.cpp:31
Generated by doxygen 1.16.1