70 double CKM_A,
double CKM_lambda,
double CKM_barrho,
double CKM_bareta,
71 double ReA7,
double ImA7,
double ReA10,
double ImA10 )
79 EvtComplex A10 = ReA10 * unit1 + ImA10 * uniti;
90 double q2 = q.
mass2();
92 double M1 = parent->
mass();
106 double Relambda_qu, Imlambda_qu;
134 Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) +
135 pow( CKM_lambda, 2.0 ) *
136 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
137 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
138 Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq;
140 Vuq = CKM_lambda * unit1;
161 Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) *
162 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
163 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
164 Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq;
166 Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) -
167 0.125 * pow( CKM_lambda, 4.0 ) );
172 <<
"\n\n The function EvtbTosllScalarAmpNew::CalcAmp(...)"
173 <<
"\n Error in the model set!"
174 <<
" ms = " << ms << std::endl;
178 Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda,
180 Vub = CKM_A * pow( CKM_lambda, 3.0 ) *
181 ( CKM_barrho * unit1 - CKM_bareta * uniti ) /
182 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
184 CKM_factor =
conj( Vtq ) * Vtb;
186 lambda_qu =
conj( Vuq ) * Vub /
188 Relambda_qu =
real( lambda_qu );
189 Imlambda_qu =
imag( lambda_qu );
202 ms, mc, mu, mt, Mw, ml,
203 Relambda_qu, Imlambda_qu );
205 mb, ms, mc, mu, mt, Mw, ml,
206 Relambda_qu, Imlambda_qu );
265 double hats = q2 / pow( M1, 2 );
266 double hatM2 = M2 / M1;
267 double hatmb = mb / M1;
268 double hatms = ms / M1;
271 EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, c, d;
273 a_b2q = c9eff_b2q * fp -
274 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 );
275 a_barb2barq = c9eff_barb2barq * fp -
276 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 );
278 b_b2q = ( c9eff_b2q * ( f0 - fp ) +
279 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 ) ) *
280 ( 1 - pow( hatM2, 2.0 ) ) / hats;
281 b_barb2barq = ( c9eff_barb2barq * ( f0 - fp ) +
282 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 ) ) *
283 ( 1 - pow( hatM2, 2.0 ) ) / hats;
287 d = c10a * ( 1.0 - pow( hatM2, 2 ) ) * ( f0 - fp ) / hats;
296 lepPlus = ( charge1 > charge2 ) ? parent->
getDaug( 1 ) : parent->
getDaug( 2 );
297 lepMinus = ( charge1 < charge2 ) ? parent->
getDaug( 1 )
309 EvtIdSet bmesons{
"B-",
"anti-B0",
"anti-B_s0",
"B_c-" };
310 EvtIdSet bbarmesons{
"B+",
"B0",
"B_s0",
"B_c+" };
314 if ( bmesons.
contains( parentID ) ) {
318 T1 = a_b2q * hatP + b_b2q * hatq;
320 T2 = c * hatP + d * hatq;
340 amp.
vertex( 0, 0, CKM_factor * ( lvc11 * T1 + lac11 * T2 ) );
341 amp.
vertex( 0, 1, CKM_factor * ( lvc12 * T1 + lac12 * T2 ) );
342 amp.
vertex( 1, 0, CKM_factor * ( lvc21 * T1 + lac21 * T2 ) );
343 amp.
vertex( 1, 1, CKM_factor * ( lvc22 * T1 + lac22 * T2 ) );
346 if ( bbarmesons.
contains( parentID ) ) {
350 T1 = a_barb2barq * hatP + b_barb2barq * hatq;
351 T2 = c * hatP + d * hatq;
371 amp.
vertex( 0, 0,
conj( CKM_factor ) * ( lvc11 * T1 + lac11 * T2 ) );
372 amp.
vertex( 0, 1,
conj( CKM_factor ) * ( lvc12 * T1 + lac12 * T2 ) );
373 amp.
vertex( 1, 0,
conj( CKM_factor ) * ( lvc21 * T1 + lac21 * T2 ) );
374 amp.
vertex( 1, 1,
conj( CKM_factor ) * ( lvc22 * T1 + lac22 * T2 ) );
379 <<
"\n\n The function EvtbTosllScalarAmpNew::CalcAmp(...)"
380 <<
"\n Wrong B-meson number" << std::endl;
399 int Nf,
int res_swch,
int ias,
double CKM_A,
double CKM_lambda,
400 double CKM_barrho,
double CKM_bareta,
double ReA7,
double ImA7,
401 double ReA10,
double ImA10 )
403 double maxfoundprob = -100.0;
409 if ( res_swch == 0 ) {
412 double t_plus, t_minus;
417 s_min = 4.0 * pow( ml, 2.0 );
418 s_max = pow( ( M1 - M2 ), 2.0 );
421 ds = ( s_max - s_min ) / ( (
double)max_j );
424 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcScalarMaxProb(...)"
425 <<
"\n ds = " << ds <<
" < 0."
426 <<
"\n s_min = " << s_min <<
"\n s_max = " << s_max
427 <<
"\n M1 = " << M1 <<
"\n M2 = " << M2
428 <<
"\n ml = " << ml << std::endl;
434 for ( j = max_j / 3; j < max_j; j++ ) {
435 s = s_min + ds * ( (double)j );
437 t_plus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
439 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
440 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
443 t_minus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
445 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
446 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
450 dt = ( t_plus - t_minus ) / ( (
double)max_k );
451 if ( fabs( dt ) < 0.00001 )
454 if ( dt <= ( -0.00001 ) ) {
456 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcScalarMaxProb(...)"
457 <<
"\n dt = " << dt <<
" < 0."
458 <<
"\n s = " << s <<
"\n s_min = " << s_min
459 <<
"\n s_max = " << s_max <<
"\n ds = " << ds
460 <<
"\n j = " << j <<
"\n t_plus = " << t_plus
461 <<
"\n t_minus = " << t_minus <<
"\n M1 = " << M1
462 <<
"\n M2 = " << M2 <<
"\n ml = " << ml
468 for ( k = 0; k < max_k; k++ ) {
469 t_for_s = t_minus + dt * ( (double)k );
471 if ( ( t_for_s > t_plus ) && ( t_for_s <= ( 1.0001 * t_plus ) ) ) {
474 if ( t_for_s > ( 1.0001 * t_plus ) ) {
476 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcScalarMaxProb(...)"
477 <<
"\n t_for_s = " << t_for_s
478 <<
" > t_plus = " << t_plus <<
" ! "
479 <<
"\n t_minus = " << t_minus <<
"\n dt = " << dt
480 <<
"\n k = " << k <<
"\n s = " << s
481 <<
"\n M1 = " << M1 <<
"\n M2 = " << M2
482 <<
"\n ml = " << ml << std::endl;
488 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s ) /
490 El2 = ( s + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) /
494 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) );
495 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
498 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) +
499 2.0 * EV * El2 - t_for_s ) /
500 ( 2.0 * modV * modl2 );
501 if ( ( fabs( cosVellminus ) > 1.0 ) &&
502 ( fabs( cosVellminus ) <= 1.0001 ) ) {
507 cosVellminus = cosVellminus / fabs( cosVellminus );
509 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) {
510 cosVellminus = cosVellminus / fabs( cosVellminus );
529 if ( fabs( cosVellminus ) > 1.0001 ) {
531 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
532 <<
"\n |cos(theta)| = " << fabs( cosVellminus ) <<
" > 1"
533 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
534 <<
"\n s_min = " << s_min <<
"\n s_max = " << s_max
535 <<
"\n t_plus = " << t_plus <<
"\n t_minus = " << t_minus
536 <<
"\n dt = " << dt <<
"\n EV = " << EV
537 <<
"\n El2 = " << El2 <<
"\n modV = " << modV
538 <<
"\n modl2 = " << modl2 <<
"\n M2 = " << M2
539 <<
"\n ml = " << ml << std::endl;
544 p1.
set( M1, 0.0, 0.0, 0.0 );
545 p2.
set( EV, modV, 0.0, 0.0 );
546 k2.
set( El2, modl2 * cosVellminus,
547 -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 );
579 scalar_part->
init( parnum, p1 );
585 listdaug[0] = mesnum;
590 amp.
init( parnum, 3, listdaug );
596 vect = root_part->
getDaug( 0 );
597 lep1 = root_part->
getDaug( 1 );
598 lep2 = root_part->
getDaug( 2 );
604 vect->
init( mesnum, p2 );
605 lep1->
init( l1num, k1 );
606 lep2->
init( l2num, k2 );
613 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf,
614 res_swch, ias, CKM_A, CKM_lambda, CKM_barrho,
615 CKM_bareta, ReA7, ImA7, ReA10, ImA10 );
620 if ( nikmax > maxfoundprob ) {
621 maxfoundprob = nikmax;
643 if ( res_swch == 1 ) {
645 double t_plus, t_minus;
649 s = pow( 3.09688, 2.0 );
651 t_plus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
652 t_plus = t_plus + sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
653 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
656 t_minus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
658 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
659 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
662 dt = ( t_plus - t_minus ) / 1000.0;
665 for ( k = 0; k < 1000; k++ ) {
666 t_for_s = t_minus + dt * ( (double)k );
668 if ( ( t_for_s > t_plus ) && ( t_for_s <= ( 1.0001 * t_plus ) ) ) {
671 if ( t_for_s > ( 1.0001 * t_plus ) ) {
673 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
674 <<
"\n t_for_s = " << t_for_s <<
" > t_plus = " << t_plus
676 <<
"\n t_minus = " << t_minus <<
"\n dt = " << dt
677 <<
"\n k = " << k <<
"\n s = " << s
678 <<
"\n M1 = " << M1 <<
"\n M2 = " << M2
679 <<
"\n ml = " << ml << std::endl;
685 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s ) /
687 El2 = ( s + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) /
691 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) );
692 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
695 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 2.0 * EV * El2 -
697 ( 2.0 * modV * modl2 );
698 if ( ( fabs( cosVellminus ) > 1.0 ) &&
699 ( fabs( cosVellminus ) <= 1.0001 ) ) {
704 cosVellminus = cosVellminus / fabs( cosVellminus );
706 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) {
707 cosVellminus = cosVellminus / fabs( cosVellminus );
724 if ( fabs( cosVellminus ) > 1.0001 ) {
726 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
727 <<
"\n |cos(theta)| = " << fabs( cosVellminus ) <<
" > 1"
728 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
729 <<
"\n t_plus = " << t_plus <<
"\n t_minus = " << t_minus
730 <<
"\n dt = " << dt <<
"\n EV = " << EV
731 <<
"\n El2 = " << El2 <<
"\n modV = " << modV
732 <<
"\n modl2 = " << modl2 <<
"\n M2 = " << M2
733 <<
"\n ml = " << ml << std::endl;
738 p1.
set( M1, 0.0, 0.0, 0.0 );
739 p2.
set( EV, modV, 0.0, 0.0 );
740 k2.
set( El2, modl2 * cosVellminus,
741 -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 );
773 scalar_part->
init( parnum, p1 );
779 listdaug[0] = mesnum;
784 amp.
init( parnum, 3, listdaug );
790 vect = root_part->
getDaug( 0 );
791 lep1 = root_part->
getDaug( 1 );
792 lep2 = root_part->
getDaug( 2 );
798 vect->
init( mesnum, p2 );
799 lep1->
init( l1num, k1 );
800 lep2->
init( l2num, k2 );
807 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch,
808 ias, CKM_A, CKM_lambda, CKM_barrho, CKM_bareta, ReA7, ImA7,
814 if ( nikmax > maxfoundprob ) {
815 maxfoundprob = nikmax;
833 if ( maxfoundprob <= 0.0 ) {
835 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
836 <<
"\n maxfoundprob = " << maxfoundprob <<
" <0 or =0!"
837 <<
"\n res_swch = " << res_swch << std::endl;
842 <<
"\n maxfoundprob (...) = " << maxfoundprob << std::endl;
844 maxfoundprob *= 1.01;
virtual void init(EvtId part_n, const EvtVector4R &p4)=0