BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtbTosllAmp.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 1998 Caltech, UCSB
10//
11// Module: EvtSemiLeptonicAmp.cc
12//
13// Description: Routine to implement semileptonic decays to pseudo-scalar
14// mesons.
15//
16// Modification history:
17//
18// DJL April 17,1998 Module created
19//
20//------------------------------------------------------------------------
21//
22#include "EvtbTosllAmp.hh"
23#include "../EvtGenBase/EvtAmp.hh"
24#include "../EvtGenBase/EvtDiracSpinor.hh"
25#include "../EvtGenBase/EvtGenKine.hh"
26#include "../EvtGenBase/EvtId.hh"
27#include "../EvtGenBase/EvtPDL.hh"
28#include "../EvtGenBase/EvtParticle.hh"
29#include "../EvtGenBase/EvtPatches.hh"
30#include "../EvtGenBase/EvtReport.hh"
31#include "../EvtGenBase/EvtScalarParticle.hh"
32#include "../EvtGenBase/EvtTensor4C.hh"
33#include "../EvtGenBase/EvtVector4C.hh"
34#include "../EvtGenBase/EvtVectorParticle.hh"
35
36double EvtbTosllAmp::CalcMaxProb( EvtId parent, EvtId meson, EvtId lepton1, EvtId lepton2,
37 EvtbTosllFF* FormFactors, double& poleSize ) {
38
39 // This routine takes the arguements parent, meson, and lepton
40 // number, and a form factor model, and returns a maximum
41 // probability for this semileptonic form factor model. A
42 // brute force method is used. The 2D cos theta lepton and
43 // q2 phase space is probed.
44
45 // Start by declaring a particle at rest.
46
47 // It only makes sense to have a scalar parent. For now.
48 // This should be generalized later.
49
50 EvtScalarParticle* scalar_part;
51 EvtParticle* root_part;
52
53 scalar_part = new EvtScalarParticle;
54
55 // cludge to avoid generating random numbers!
56 scalar_part->noLifeTime();
57
58 EvtVector4R p_init;
59
60 p_init.set( EvtPDL::getMass( parent ), 0.0, 0.0, 0.0 );
61 scalar_part->init( parent, p_init );
62 root_part = (EvtParticle*)scalar_part;
63 root_part->setDiagonalSpinDensity();
64
65 EvtParticle *daughter, *lep1, *lep2;
66
67 EvtAmp amp;
68
69 EvtId listdaug[3];
70 listdaug[0] = meson;
71 listdaug[1] = lepton1;
72 listdaug[2] = lepton2;
73
74 amp.init( parent, 3, listdaug );
75
76 root_part->makeDaughters( 3, listdaug );
77 daughter = root_part->getDaug( 0 );
78 lep1 = root_part->getDaug( 1 );
79 lep2 = root_part->getDaug( 2 );
80
81 // cludge to avoid generating random numbers!
82 daughter->noLifeTime();
83 lep1->noLifeTime();
84 lep2->noLifeTime();
85
86 // Initial particle is unpolarized, well it is a scalar so it is
87 // trivial
88 EvtSpinDensity rho;
89 rho.SetDiag( root_part->getSpinStates() );
90
91 double mass[3];
92
93 double m = root_part->mass();
94
95 EvtVector4R p4meson, p4lepton1, p4lepton2, p4w;
96 double q2max;
97
98 double q2, elepton, plepton;
99 int i, j;
100 double erho, prho, costl;
101
102 double maxfoundprob = 0.0;
103 double prob = -10.0;
104 int massiter;
105
106 double maxpole = 0;
107
108 for ( massiter = 0; massiter < 3; massiter++ )
109 {
110
111 mass[0] = EvtPDL::getMeanMass( meson );
112 mass[1] = EvtPDL::getMeanMass( lepton1 );
113 mass[2] = EvtPDL::getMeanMass( lepton2 );
114 if ( massiter == 1 ) { mass[0] = EvtPDL::getMinMass( meson ); }
115 if ( massiter == 2 )
116 {
117 mass[0] = EvtPDL::getMaxMass( meson );
118 if ( ( mass[0] + mass[1] + mass[2] ) > m ) mass[0] = m - mass[1] - mass[2] - 0.00001;
119 }
120
121 q2max = ( m - mass[0] ) * ( m - mass[0] );
122
123 // loop over q2
124 // cout << "m " << m << "mass[0] " << mass[0] << " q2max "<< q2max << endl;
125 for ( i = 0; i < 25; i++ )
126 {
127 // want to avoid picking up the tail of the photon propagator
128 q2 = ( ( i + 1.5 ) * q2max ) / 26.0;
129
130 if ( i == 0 ) q2 = 4 * ( mass[1] * mass[1] );
131
132 erho = ( m * m + mass[0] * mass[0] - q2 ) / ( 2.0 * m );
133
134 prho = sqrt( erho * erho - mass[0] * mass[0] );
135
136 p4meson.set( erho, 0.0, 0.0, -1.0 * prho );
137 p4w.set( m - erho, 0.0, 0.0, prho );
138
139 // This is in the W rest frame
140 elepton = ( q2 + mass[1] * mass[1] ) / ( 2.0 * sqrt( q2 ) );
141 plepton = sqrt( elepton * elepton - mass[1] * mass[1] );
142
143 double probctl[3];
144
145 for ( j = 0; j < 3; j++ )
146 {
147
148 costl = 0.99 * ( j - 1.0 );
149
150 // These are in the W rest frame. Need to boost out into
151 // the B frame.
152 p4lepton1.set( elepton, 0.0, plepton * sqrt( 1.0 - costl * costl ), plepton * costl );
153 p4lepton2.set( elepton, 0.0, -1.0 * plepton * sqrt( 1.0 - costl * costl ),
154 -1.0 * plepton * costl );
155
156 EvtVector4R boost( ( m - erho ), 0.0, 0.0, 1.0 * prho );
157 p4lepton1 = boostTo( p4lepton1, boost );
158 p4lepton2 = boostTo( p4lepton2, boost );
159
160 // Now initialize the daughters...
161
162 daughter->init( meson, p4meson );
163 lep1->init( lepton1, p4lepton1 );
164 lep2->init( lepton2, p4lepton2 );
165
166 CalcAmp( root_part, amp, FormFactors );
167
168 // Now find the probability at this q2 and cos theta lepton point
169 // and compare to maxfoundprob.
170
171 // Do a little magic to get the probability!!
172
173 // cout <<"amp:"<<amp.getSpinDensity()<<endl;
174
175 prob = rho.NormalizedProb( amp.getSpinDensity() );
176
177 // cout << "prob:"<<q2<<" "<<costl<<" "<<prob<<endl;
178
179 probctl[j] = prob;
180 }
181
182 // probclt contains prob at ctl=-1,0,1.
183 // prob=a+b*ctl+c*ctl^2
184
185 double a = probctl[1];
186 double b = 0.5 * ( probctl[2] - probctl[0] );
187 double c = 0.5 * ( probctl[2] + probctl[0] ) - probctl[1];
188
189 prob = probctl[0];
190 if ( probctl[1] > prob ) prob = probctl[1];
191 if ( probctl[2] > prob ) prob = probctl[2];
192
193 if ( fabs( c ) > 1e-20 )
194 {
195 double ctlx = -0.5 * b / c;
196 if ( fabs( ctlx ) < 1.0 )
197 {
198 double probtmp = a + b * ctlx + c * ctlx * ctlx;
199 if ( probtmp > prob ) prob = probtmp;
200 }
201 }
202
203 // report(DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
204 // << probctl[0]<<" "
205 // << probctl[1]<<" "
206 // << probctl[2]<<endl;
207
208 if ( i == 0 )
209 {
210 maxpole = prob;
211 continue;
212 }
213
214 if ( prob > maxfoundprob ) { maxfoundprob = prob; }
215
216 // cout << "q2,maxfoundprob:"<<q2<<" "<<maxfoundprob<<endl;
217 }
218 if ( EvtPDL::getWidth( meson ) <= 0.0 )
219 {
220 // if the particle is narrow dont bother with changing the mass.
221 massiter = 4;
222 }
223 }
224
225 root_part->deleteTree();
226
227 poleSize = 0.04 * ( maxpole / maxfoundprob ) * 4 * ( mass[1] * mass[1] );
228
229 // poleSize=0.002;
230
231 // cout <<"maxfoundprob,maxpole,poleSize:"<<maxfoundprob<<" "
232 // <<maxpole<<" "<<poleSize<<endl;
233
234 maxfoundprob *= 1.15;
235
236 return maxfoundprob;
237}
238
239EvtComplex EvtbTosllAmp::GetC7Eff( double q2, bool nnlo ) {
240
241 if ( !nnlo ) return -0.313;
242 double mbeff = 4.8;
243 double shat = q2 / mbeff / mbeff;
244 double logshat;
245 logshat = log( shat );
246
247 double muscale;
248 muscale = 2.5;
249 double alphas;
250 alphas = 0.267;
251 double A7;
252 A7 = -0.353 + 0.023;
253 double A8;
254 A8 = -0.164;
255 double A9;
256 A9 = 4.287 + ( -0.218 );
257 double A10;
258 A10 = -4.592 + 0.379;
259 double C1;
260 C1 = -0.697;
261 double C2;
262 C2 = 1.046;
263 double T9;
264 T9 = 0.114 + 0.280;
265 double U9;
266 U9 = 0.045 + 0.023;
267 double W9;
268 W9 = 0.044 + 0.016;
269
270 double Lmu;
271 Lmu = log( muscale / mbeff );
272
273 EvtComplex uniti( 0.0, 1.0 );
274
275 EvtComplex c7eff;
276 if ( shat > 0.25 )
277 {
278 c7eff = A7;
279 return c7eff;
280 }
281
282 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
283 muscale = 5.0;
284 alphas = 0.215;
285 A7 = -0.312 + 0.008;
286 A8 = -0.148;
287 A9 = 4.174 + ( -0.035 );
288 A10 = -4.592 + 0.379;
289 C1 = -0.487;
290 C2 = 1.024;
291 T9 = 0.374 + 0.252;
292 U9 = 0.033 + 0.015;
293 W9 = 0.032 + 0.012;
294 Lmu = log( muscale / mbeff );
295
296 EvtComplex F71;
297 EvtComplex f71;
298 EvtComplex k7100( -0.68192, -0.074998 );
299 EvtComplex k7101( 0.0, 0.0 );
300 EvtComplex k7110( -0.23935, -0.12289 );
301 EvtComplex k7111( 0.0027424, 0.019676 );
302 EvtComplex k7120( -0.0018555, -0.175 );
303 EvtComplex k7121( 0.022864, 0.011456 );
304 EvtComplex k7130( 0.28248, -0.12783 );
305 EvtComplex k7131( 0.029027, -0.0082265 );
306 f71 = k7100 + k7101 * logshat + shat * ( k7110 + k7111 * logshat ) +
307 shat * shat * ( k7120 + k7121 * logshat ) +
308 shat * shat * shat * ( k7130 + k7131 * logshat );
309 F71 = ( -208.0 / 243.0 ) * Lmu + f71;
310
311 EvtComplex F72;
312 EvtComplex f72;
313 EvtComplex k7200( 4.0915, 0.44999 );
314 EvtComplex k7201( 0.0, 0.0 );
315 EvtComplex k7210( 1.4361, 0.73732 );
316 EvtComplex k7211( -0.016454, -0.11806 );
317 EvtComplex k7220( 0.011133, 1.05 );
318 EvtComplex k7221( -0.13718, -0.068733 );
319 EvtComplex k7230( -1.6949, 0.76698 );
320 EvtComplex k7231( -0.17416, 0.049359 );
321 f72 = k7200 + k7201 * logshat + shat * ( k7210 + k7211 * logshat ) +
322 shat * shat * ( k7220 + k7221 * logshat ) +
323 shat * shat * shat * ( k7230 + k7231 * logshat );
324 F72 = ( 416.0 / 81.0 ) * Lmu + f72;
325
326 EvtComplex F78;
327 F78 = ( -32.0 / 9.0 ) * Lmu + 8.0 * EvtConst::pi * EvtConst::pi / 27.0 + ( -44.0 / 9.0 ) +
328 ( -8.0 * EvtConst::pi / 9.0 ) * uniti +
329 ( 4.0 / 3.0 * EvtConst::pi * EvtConst::pi - 40.0 / 3.0 ) * shat +
330 ( 32.0 * EvtConst::pi * EvtConst::pi / 9.0 - 316.0 / 9.0 ) * shat * shat +
331 ( 200.0 * EvtConst::pi * EvtConst::pi / 27.0 - 658.0 / 9.0 ) * shat * shat * shat +
332 ( -8.0 * logshat / 9.0 ) * ( shat + shat * shat + shat * shat * shat );
333
334 c7eff = A7 - alphas / ( 4.0 * EvtConst::pi ) * ( C1 * F71 + C2 * F72 + A8 * F78 );
335
336 return c7eff;
337}
338
339EvtComplex EvtbTosllAmp::GetC9Eff( double q2, bool nnlo, bool btod ) {
340
341 if ( !nnlo ) return 4.344;
342 double mbeff = 4.8;
343 double shat = q2 / mbeff / mbeff;
344 double logshat;
345 logshat = log( shat );
346 double mchat = 0.29;
347
348 double muscale;
349 muscale = 2.5;
350 double alphas;
351 alphas = 0.267;
352 double A7;
353 A7 = -0.353 + 0.023;
354 double A8;
355 A8 = -0.164;
356 double A9;
357 A9 = 4.287 + ( -0.218 );
358 double A10;
359 A10 = -4.592 + 0.379;
360 double C1;
361 C1 = -0.697;
362 double C2;
363 C2 = 1.046;
364 double T9;
365 T9 = 0.114 + 0.280;
366 double U9;
367 U9 = 0.045 + 0.023;
368 double W9;
369 W9 = 0.044 + 0.016;
370
371 double Lmu;
372 Lmu = log( muscale / mbeff );
373
374 EvtComplex uniti( 0.0, 1.0 );
375
376 EvtComplex hc;
377 double xarg;
378 xarg = 4.0 * mchat / shat;
379 hc = -4.0 / 9.0 * log( mchat * mchat ) + 8.0 / 27.0 + 4.0 * xarg / 9.0;
380
381 if ( xarg < 1.0 )
382 {
383 hc =
384 hc - 2.0 / 9.0 * ( 2.0 + xarg ) * sqrt( fabs( 1.0 - xarg ) ) *
385 ( log( fabs( ( sqrt( 1.0 - xarg ) + 1.0 ) / ( sqrt( 1.0 - xarg ) - 1.0 ) ) ) -
386 uniti * EvtConst::pi );
387 }
388 else
389 {
390 hc = hc - 2.0 / 9.0 * ( 2.0 + xarg ) * sqrt( fabs( 1.0 - xarg ) ) * 2.0 *
391 atan( 1.0 / sqrt( xarg - 1.0 ) );
392 }
393
394 EvtComplex h1;
395 xarg = 4.0 / shat;
396 h1 = 8.0 / 27.0 + 4.0 * xarg / 9.0;
397 if ( xarg < 1.0 )
398 {
399 h1 =
400 h1 - 2.0 / 9.0 * ( 2.0 + xarg ) * sqrt( fabs( 1.0 - xarg ) ) *
401 ( log( fabs( ( sqrt( 1.0 - xarg ) + 1.0 ) / ( sqrt( 1.0 - xarg ) - 1.0 ) ) ) -
402 uniti * EvtConst::pi );
403 }
404 else
405 {
406 h1 = h1 - 2.0 / 9.0 * ( 2.0 + xarg ) * sqrt( fabs( 1.0 - xarg ) ) * 2.0 *
407 atan( 1.0 / sqrt( xarg - 1.0 ) );
408 }
409
410 EvtComplex h0;
411 h0 = 8.0 / 27.0 - 4.0 * log( 2.0 ) / 9.0 + 4.0 * uniti * EvtConst::pi / 9.0;
412
413 // X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
414 // h(\hat m_u^2, hat s))
415 EvtComplex Vudstar( 1.0 - 0.2279 * 0.2279 / 2.0, 0.0 );
416 EvtComplex Vub( ( 0.118 + 0.273 ) / 2.0, -1.0 * ( 0.305 + 0.393 ) / 2.0 );
417 EvtComplex Vtdstar( 1.0 - ( 0.118 + 0.273 ) / 2.0, ( 0.305 + 0.393 ) / 2.0 );
418 EvtComplex Vtb( 1.0, 0.0 );
419
420 EvtComplex Xd;
421 Xd = ( Vudstar * Vub / Vtdstar * Vtb ) * ( 4.0 / 3.0 * C1 + C2 ) * ( hc - h0 );
422
423 EvtComplex c9eff = 4.344;
424 if ( shat > 0.25 )
425 {
426 c9eff = A9 + T9 * hc + U9 * h1 + W9 * h0;
427 if ( btod ) { c9eff += Xd; }
428
429 return c9eff;
430 }
431
432 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
433 muscale = 5.0;
434 alphas = 0.215;
435 A9 = 4.174 + ( -0.035 );
436 C1 = -0.487;
437 C2 = 1.024;
438 A8 = -0.148;
439 T9 = 0.374 + 0.252;
440 U9 = 0.033 + 0.015;
441 W9 = 0.032 + 0.012;
442 Lmu = log( muscale / mbeff );
443
444 EvtComplex F91;
445 EvtComplex f91;
446 EvtComplex k9100( -11.973, 0.16371 );
447 EvtComplex k9101( -0.081271, -0.059691 );
448 EvtComplex k9110( -28.432, -0.25044 );
449 EvtComplex k9111( -0.040243, 0.016442 );
450 EvtComplex k9120( -57.114, -0.86486 );
451 EvtComplex k9121( -0.035191, 0.027909 );
452 EvtComplex k9130( -128.8, -2.5243 );
453 EvtComplex k9131( -0.017587, 0.050639 );
454 f91 = k9100 + k9101 * logshat + shat * ( k9110 + k9111 * logshat ) +
455 shat * shat * ( k9120 + k9121 * logshat ) +
456 shat * shat * shat * ( k9130 + k9131 * logshat );
457 F91 =
458 ( -1424.0 / 729.0 + 16.0 * uniti * EvtConst::pi / 243.0 + 64.0 / 27.0 * log( mchat ) ) *
459 Lmu -
460 16.0 * Lmu * logshat / 243.0 +
461 ( 16.0 / 1215.0 - 32.0 / 135.0 / mchat / mchat ) * Lmu * shat +
462 ( 4.0 / 2835.0 - 8.0 / 315.0 / mchat / mchat / mchat / mchat ) * Lmu * shat * shat +
463 ( 16.0 / 76545.0 - 32.0 / 8505.0 / mchat / mchat / mchat / mchat / mchat / mchat ) *
464 Lmu * shat * shat * shat -
465 256.0 * Lmu * Lmu / 243.0 + f91;
466
467 EvtComplex F92;
468 EvtComplex f92;
469 EvtComplex k9200( 6.6338, -0.98225 );
470 EvtComplex k9201( 0.48763, 0.35815 );
471 EvtComplex k9210( 3.3585, 1.5026 );
472 EvtComplex k9211( 0.24146, -0.098649 );
473 EvtComplex k9220( -1.1906, 5.1892 );
474 EvtComplex k9221( 0.21115, -0.16745 );
475 EvtComplex k9230( -17.12, 15.146 );
476 EvtComplex k9231( 0.10552, -0.30383 );
477 f92 = k9200 + k9201 * logshat + shat * ( k9210 + k9211 * logshat ) +
478 shat * shat * ( k9220 + k9221 * logshat ) +
479 shat * shat * shat * ( k9230 + k9231 * logshat );
480 F92 = ( 256.0 / 243.0 - 32.0 * uniti * EvtConst::pi / 81.0 - 128.0 / 9.0 * log( mchat ) ) *
481 Lmu +
482 32.0 * Lmu * logshat / 81.0 +
483 ( -32.0 / 405.0 + 64.0 / 45.0 / mchat / mchat ) * Lmu * shat +
484 ( -8.0 / 945.0 + 16.0 / 105.0 / mchat / mchat / mchat / mchat ) * Lmu * shat * shat +
485 ( -32.0 / 25515.0 + 64.0 / 2835.0 / mchat / mchat / mchat / mchat / mchat / mchat ) *
486 Lmu * shat * shat * shat +
487 512.0 * Lmu * Lmu / 81.0 + f92;
488
489 EvtComplex F98;
490 F98 =
491 104.0 / 9.0 - 32.0 * EvtConst::pi * EvtConst::pi / 27.0 +
492 ( 1184.0 / 27.0 - 40.0 * EvtConst::pi * EvtConst::pi / 9.0 ) * shat +
493 ( 14212.0 / 135.0 - 32.0 * EvtConst::pi * EvtConst::pi / 3.0 ) * shat * shat +
494 ( 193444.0 / 945.0 - 560.0 * EvtConst::pi * EvtConst::pi / 27.0 ) * shat * shat * shat +
495 16.0 * logshat / 9.0 * ( 1.0 + shat + shat * shat + shat * shat * shat );
496
497 Xd = ( Vudstar * Vub / Vtdstar * Vtb ) * ( 4.0 / 3.0 * C1 + C2 ) * ( hc - h0 );
498
499 c9eff = A9 + T9 * hc + U9 * h1 + W9 * h0 -
500 alphas / ( 4.0 * EvtConst::pi ) * ( C1 * F91 + C2 * F92 + A8 * F98 );
501 if ( btod ) { c9eff += Xd; }
502
503 return c9eff;
504}
505
506EvtComplex EvtbTosllAmp::GetC10Eff( double q2, bool nnlo ) {
507
508 if ( !nnlo ) return -4.669;
509 double A10;
510 A10 = -4.592 + 0.379;
511
512 EvtComplex c10eff;
513 c10eff = A10;
514
515 return c10eff;
516}
517
518double EvtbTosllAmp::dGdsProb( double mb, double ms, double ml, double s ) {
519 // Compute the decay probability density function given a value of s
520 // according to Ali's paper
521
522 double delta, lambda, prob;
523 double f1, f2, f3, f4;
524 double msh, mlh, sh;
525
526 mlh = ml / mb;
527 msh = ms / mb;
528 sh = s / ( mb * mb );
529
530 EvtComplex c9eff = EvtbTosllAmp::GetC9Eff( sh * mb );
531 EvtComplex c7eff = EvtbTosllAmp::GetC7Eff( sh * mb );
532 EvtComplex c10eff = EvtbTosllAmp::GetC10Eff( sh * mb );
533
534 double alphas = 0.119 / ( 1 + 0.119 * log( pow( 4.8, 2 ) / pow( 91.1867, 2 ) ) * 23.0 /
535 12.0 / EvtConst::pi );
536 double omega9 =
537 -2.0 / 9.0 * EvtConst::pi * EvtConst::pi - 4.0 / 3.0 * ddilog_( sh ) -
538 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
539 ( 5.0 + 4.0 * sh ) / ( 3.0 * ( 1.0 + 2.0 * sh ) ) * log( 1.0 - sh ) -
540 2.0 * sh * ( 1.0 + sh ) * ( 1.0 - 2.0 * sh ) /
541 ( 3.0 * pow( 1.0 - sh, 2 ) * ( 1.0 + 2.0 * sh ) ) * log( sh ) +
542 ( 5.0 + 9.0 * sh - 6.0 * sh * sh ) / ( 6.0 * ( 1.0 - sh ) * ( 1.0 + 2.0 * sh ) );
543 double eta9 = 1.0 + alphas * omega9 / EvtConst::pi;
544 double omega7 = -8.0 / 3.0 * log( 4.8 / mb ) - 4.0 / 3.0 * ddilog_( sh ) -
545 2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
546 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
547 log( 1 - sh ) * ( 8.0 + sh ) / ( 2.0 + sh ) / 3.0 -
548 2.0 / 3.0 * sh * ( 2.0 - 2.0 * sh - sh * sh ) * log( sh ) /
549 pow( ( 1.0 - sh ), 2 ) / ( 2.0 + sh ) -
550 ( 16.0 - 11.0 * sh - 17.0 * sh * sh ) / 18.0 / ( 2.0 + sh ) / ( 1.0 - sh );
551 double eta7 = 1.0 + alphas * omega7 / EvtConst::pi;
552
553 double omega79 = -4.0 / 3.0 * log( 4.8 / mb ) - 4.0 / 3.0 * ddilog_( sh ) -
554 2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
555 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
556 1.0 / 9.0 * ( 2.0 + 7.0 * sh ) * log( 1.0 - sh ) / sh -
557 2.0 / 9.0 * sh * ( 3.0 - 2.0 * sh ) * log( sh ) / pow( ( 1.0 - sh ), 2 ) +
558 1.0 / 18.0 * ( 5.0 - 9.0 * sh ) / ( 1.0 - sh );
559 double eta79 = 1.0 + alphas * omega79 / EvtConst::pi;
560
561 double c7c9 = abs( c7eff ) * real( c9eff );
562 c7c9 *= pow( eta79, 2 );
563 double c7c7 = pow( abs( c7eff ), 2 );
564 c7c7 *= pow( eta7, 2 );
565
566 double c9c9plusc10c10 = pow( abs( c9eff ), 2 ) + pow( abs( c10eff ), 2 );
567 c9c9plusc10c10 *= pow( eta9, 2 );
568 double c9c9minusc10c10 = pow( abs( c9eff ), 2 ) - pow( abs( c10eff ), 2 );
569 c9c9minusc10c10 *= pow( eta9, 2 );
570
571 lambda = 1.0 + sh * sh + pow( msh, 4 ) - 2.0 * ( sh + sh * msh * msh + msh * msh );
572
573 f1 = pow( 1.0 - msh * msh, 2 ) - sh * ( 1.0 + msh * msh );
574 f2 = 2.0 * ( 1.0 + msh * msh ) * pow( 1.0 - msh * msh, 2 ) -
575 sh * ( 1.0 + 14.0 * msh * msh + pow( msh, 4 ) ) - sh * sh * ( 1.0 + msh * msh );
576 f3 = pow( 1.0 - msh * msh, 2 ) + sh * ( 1.0 + msh * msh ) - 2.0 * sh * sh +
577 lambda * 2.0 * mlh * mlh / sh;
578 f4 = 1.0 - sh + msh * msh;
579
580 delta = ( 12.0 * c7c9 * f1 + 4.0 * c7c7 * f2 / sh ) * ( 1.0 + 2.0 * mlh * mlh / sh ) +
581 c9c9plusc10c10 * f3 + 6.0 * mlh * mlh * c9c9minusc10c10 * f4;
582
583 prob = sqrt( lambda * ( 1.0 - 4.0 * mlh * mlh / sh ) ) * delta;
584
585 return prob;
586}
587
588double EvtbTosllAmp::dGdsdupProb( double mb, double ms, double ml, double s, double u ) {
589 // Compute the decay probability density function given a value of s and u
590 // according to Ali's paper
591
592 double prob;
593 double f1sp, f2sp, f3sp;
594
595 double sh = s / ( mb * mb );
596
597 EvtComplex c9eff = EvtbTosllAmp::GetC9Eff( sh * mb );
598 EvtComplex c7eff = EvtbTosllAmp::GetC7Eff( sh * mb );
599 EvtComplex c10eff = EvtbTosllAmp::GetC10Eff( sh * mb );
600
601 double alphas = 0.119 / ( 1 + 0.119 * log( pow( 4.8, 2 ) / pow( 91.1867, 2 ) ) * 23.0 /
602 12.0 / EvtConst::pi );
603 double omega9 =
604 -2.0 / 9.0 * EvtConst::pi * EvtConst::pi - 4.0 / 3.0 * ddilog_( sh ) -
605 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
606 ( 5.0 + 4.0 * sh ) / ( 3.0 * ( 1.0 + 2.0 * sh ) ) * log( 1.0 - sh ) -
607 2.0 * sh * ( 1.0 + sh ) * ( 1.0 - 2.0 * sh ) /
608 ( 3.0 * pow( 1.0 - sh, 2 ) * ( 1.0 + 2.0 * sh ) ) * log( sh ) +
609 ( 5.0 + 9.0 * sh - 6.0 * sh * sh ) / ( 6.0 * ( 1.0 - sh ) * ( 1.0 + 2.0 * sh ) );
610 double eta9 = 1.0 + alphas * omega9 / EvtConst::pi;
611 double omega7 = -8.0 / 3.0 * log( 4.8 / mb ) - 4.0 / 3.0 * ddilog_( sh ) -
612 2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
613 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
614 log( 1 - sh ) * ( 8.0 + sh ) / ( 2.0 + sh ) / 3.0 -
615 2.0 / 3.0 * sh * ( 2.0 - 2.0 * sh - sh * sh ) * log( sh ) /
616 pow( ( 1.0 - sh ), 2 ) / ( 2.0 + sh ) -
617 ( 16.0 - 11.0 * sh - 17.0 * sh * sh ) / 18.0 / ( 2.0 + sh ) / ( 1.0 - sh );
618 double eta7 = 1.0 + alphas * omega7 / EvtConst::pi;
619
620 double omega79 = -4.0 / 3.0 * log( 4.8 / mb ) - 4.0 / 3.0 * ddilog_( sh ) -
621 2.0 / 9.0 * EvtConst::pi * EvtConst::pi -
622 2.0 / 3.0 * log( sh ) * log( 1.0 - sh ) -
623 1.0 / 9.0 * ( 2.0 + 7.0 * sh ) * log( 1.0 - sh ) / sh -
624 2.0 / 9.0 * sh * ( 3.0 - 2.0 * sh ) * log( sh ) / pow( ( 1.0 - sh ), 2 ) +
625 1.0 / 18.0 * ( 5.0 - 9.0 * sh ) / ( 1.0 - sh );
626 double eta79 = 1.0 + alphas * omega79 / EvtConst::pi;
627
628 double c7c9 = abs( c7eff ) * real( c9eff );
629 c7c9 *= pow( eta79, 2 );
630 double c7c7 = pow( abs( c7eff ), 2 );
631 c7c7 *= pow( eta7, 2 );
632
633 double c9c9plusc10c10 = pow( abs( c9eff ), 2 ) + pow( abs( c10eff ), 2 );
634 c9c9plusc10c10 *= pow( eta9, 2 );
635 double c9c9minusc10c10 = pow( abs( c9eff ), 2 ) - pow( abs( c10eff ), 2 );
636 c9c9minusc10c10 *= pow( eta9, 2 );
637 double c7c10 = abs( c7eff ) * real( c10eff );
638 c7c10 *= eta7;
639 c7c10 *= eta9;
640 double c9c10 = real( c9eff ) * real( c10eff );
641 c9c10 *= pow( eta9, 2 );
642
643 f1sp =
644 ( pow( mb * mb - ms * ms, 2 ) - s * s ) * c9c9plusc10c10 +
645 4.0 *
646 ( pow( mb, 4 ) - ms * ms * mb * mb - pow( ms, 4 ) * ( 1.0 - ms * ms / ( mb * mb ) ) -
647 8.0 * s * ms * ms - s * s * ( 1.0 + ms * ms / ( mb * mb ) ) ) *
648 mb * mb * c7c7 /
649 s
650 // kludged mass term
651 * ( 1.0 + 2.0 * ml * ml / s ) -
652 8.0 * ( s * ( mb * mb + ms * ms ) - pow( mb * mb - ms * ms, 2 ) ) *
653 c7c9
654 // kludged mass term
655 * ( 1.0 + 2.0 * ml * ml / s );
656
657 f2sp = 4.0 * s * c9c10 + 8.0 * ( mb * mb + ms * ms ) * c7c10;
658 f3sp = -( c9c9plusc10c10 ) + 4.0 * ( 1.0 + pow( ms / mb, 4 ) ) * mb * mb * c7c7 /
659 s
660 // kludged mass term
661 * ( 1.0 + 2.0 * ml * ml / s );
662
663 prob = ( f1sp + f2sp * u + f3sp * u * u ) / pow( mb, 3 );
664
665 return prob;
666}
mass
double mass
Definition CosmicGenerator.cxx:128
f1
TFile * f1
Definition DataBase/tau_mode.c:4
ml
character *LEPTONflag integer iresonances real zeta5 real adp3 real large_3 real zeta5 common params adp3 common switch large_3 common lepton LEPTONflag common RESFIT IRESON common RES iresonances common alpgmu era0 common physparams ERMW common leptomass ml
Definition Ekhara/src/ekhara/common.h:15
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition EvtDiracSpinor.cc:67
ddilog_
double ddilog_(const double &sh)
s
XmlRpcServer s
Definition HelloServer.cpp:11
EvtAmp::init
void init(EvtId p, int ndaug, EvtId *daug)
Definition EvtAmp.cc:65
EvtAmp::getSpinDensity
EvtSpinDensity getSpinDensity()
Definition EvtAmp.cc:135
EvtConst::pi
static const double pi
Definition EvtConst.hh:27
EvtId
Definition EvtId.hh:27
EvtPDL::getWidth
static double getWidth(EvtId i)
Definition EvtPDL.hh:59
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition EvtPDL.hh:43
EvtPDL::getMinMass
static double getMinMass(EvtId i)
Definition EvtPDL.hh:56
EvtPDL::getMaxMass
static double getMaxMass(EvtId i)
Definition EvtPDL.hh:55
EvtPDL::getMass
static double getMass(EvtId i)
Definition EvtPDL.hh:44
EvtParticle::noLifeTime
void noLifeTime()
Definition EvtParticle.hh:358
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition EvtParticle.cc:1150
EvtParticle::init
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtParticle::getSpinStates
int getSpinStates() const
Definition EvtParticle.cc:104
EvtParticle::setDiagonalSpinDensity
void setDiagonalSpinDensity()
Definition EvtParticle.cc:116
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::mass
double mass() const
Definition EvtParticle.cc:114
EvtParticle::deleteTree
void deleteTree()
Definition EvtParticle.cc:575
EvtScalarParticle::init
void init(EvtId part_n, double e, double px, double py, double pz)
Definition EvtScalarParticle.cc:27
EvtSpinDensity::NormalizedProb
double NormalizedProb(const EvtSpinDensity &d)
Definition EvtSpinDensity.cc:112
EvtSpinDensity::SetDiag
void SetDiag(int n)
Definition EvtSpinDensity.cc:101
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:173
EvtbTosllAmp::GetC7Eff
EvtComplex GetC7Eff(double q2, bool nnlo=true)
Definition EvtbTosllAmp.cc:239
EvtbTosllAmp::GetC10Eff
EvtComplex GetC10Eff(double q2, bool nnlo=true)
Definition EvtbTosllAmp.cc:506
EvtbTosllAmp::dGdsProb
double dGdsProb(double mb, double ms, double ml, double s)
Definition EvtbTosllAmp.cc:518
EvtbTosllAmp::CalcAmp
virtual void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtbTosllFF *formFactors)=0
EvtbTosllAmp::dGdsdupProb
double dGdsdupProb(double mb, double ms, double ml, double s, double u)
Definition EvtbTosllAmp.cc:588
EvtbTosllAmp::GetC9Eff
EvtComplex GetC9Eff(double q2, bool nnlo=true, bool btod=false)
Definition EvtbTosllAmp.cc:339
EvtbTosllAmp::CalcMaxProb
double CalcMaxProb(EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtbTosllFF *formFactors, double &poleSize)
Definition EvtbTosllAmp.cc:36