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