70 double mu,
int Nf,
int res_swch,
int ias,
71 double CKM_A,
double CKM_lambda,
72 double CKM_barrho,
double CKM_bareta )
88 double q2 = q.
mass2();
90 double M1 = parent->
mass();
104 double Relambda_qu, Imlambda_qu;
132 Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) +
133 pow( CKM_lambda, 2.0 ) *
134 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
135 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
136 Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq;
138 Vuq = CKM_lambda * unit1;
166 Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) *
167 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
168 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
169 Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq;
171 Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) -
172 0.125 * pow( CKM_lambda, 4.0 ) );
177 <<
"\n\n The function EvtbTosllScalarAmpNew::CalcAmp(...)"
178 <<
"\n Error in the model set!"
179 <<
" ms = " << ms << std::endl;
183 Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda,
185 Vub = CKM_A * pow( CKM_lambda, 3.0 ) *
186 ( CKM_barrho * unit1 - CKM_bareta * uniti ) /
187 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
189 CKM_factor =
conj( Vtq ) * Vtb;
191 lambda_qu =
conj( Vuq ) * Vub /
193 Relambda_qu =
real( lambda_qu );
194 Imlambda_qu =
imag( lambda_qu );
206 ms, mc, mu, mt, Mw, ml,
207 Relambda_qu, Imlambda_qu );
209 mb, ms, mc, mu, mt, Mw, ml,
210 Relambda_qu, Imlambda_qu );
268 double hats = q2 / pow( M1, 2 );
269 double hatM2 = M2 / M1;
270 double hatmb = mb / M1;
271 double hatms = ms / M1;
274 EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, c, d;
276 a_b2q = c9eff_b2q * fp -
277 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 );
278 a_barb2barq = c9eff_barb2barq * fp -
279 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 );
281 b_b2q = ( c9eff_b2q * ( f0 - fp ) +
282 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 ) ) *
283 ( 1 - pow( hatM2, 2.0 ) ) / hats;
284 b_barb2barq = ( c9eff_barb2barq * ( f0 - fp ) +
285 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 ) ) *
286 ( 1 - pow( hatM2, 2.0 ) ) / hats;
290 d = c10a * ( 1.0 - pow( hatM2, 2 ) ) * ( f0 - fp ) / hats;
299 lepPlus = ( charge1 > charge2 ) ? parent->
getDaug( 1 ) : parent->
getDaug( 2 );
300 lepMinus = ( charge1 < charge2 ) ? parent->
getDaug( 1 )
312 EvtIdSet bmesons{
"B-",
"anti-B0",
"anti-B_s0",
"B_c-" };
313 EvtIdSet bbarmesons{
"B+",
"B0",
"B_s0",
"B_c+" };
317 if ( bmesons.
contains( parentID ) ) {
321 T1 = a_b2q * hatP + b_b2q * hatq;
323 T2 = c * hatP + d * hatq;
343 amp.
vertex( 0, 0, CKM_factor * ( lvc11 * T1 + lac11 * T2 ) );
344 amp.
vertex( 0, 1, CKM_factor * ( lvc12 * T1 + lac12 * T2 ) );
345 amp.
vertex( 1, 0, CKM_factor * ( lvc21 * T1 + lac21 * T2 ) );
346 amp.
vertex( 1, 1, CKM_factor * ( lvc22 * T1 + lac22 * T2 ) );
349 if ( bbarmesons.
contains( parentID ) ) {
353 T1 = a_barb2barq * hatP + b_barb2barq * hatq;
354 T2 = c * hatP + d * hatq;
374 amp.
vertex( 0, 0,
conj( CKM_factor ) * ( lvc11 * T1 + lac11 * T2 ) );
375 amp.
vertex( 0, 1,
conj( CKM_factor ) * ( lvc12 * T1 + lac12 * T2 ) );
376 amp.
vertex( 1, 0,
conj( CKM_factor ) * ( lvc21 * T1 + lac21 * T2 ) );
377 amp.
vertex( 1, 1,
conj( CKM_factor ) * ( lvc22 * T1 + lac22 * T2 ) );
382 <<
"\n\n The function EvtbTosllScalarAmpNew::CalcAmp(...)"
383 <<
"\n Wrong B-meson number" << std::endl;
402 int Nf,
int res_swch,
int ias,
double CKM_A,
double CKM_lambda,
403 double CKM_barrho,
double CKM_bareta )
405 double maxfoundprob = -100.0;
411 if ( res_swch == 0 ) {
414 double t_plus, t_minus;
419 s_min = 4.0 * pow( ml, 2.0 );
420 s_max = pow( ( M1 - M2 ), 2.0 );
423 ds = ( s_max - s_min ) / ( (
double)max_j );
426 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcScalarMaxProb(...)"
427 <<
"\n ds = " << ds <<
" < 0."
428 <<
"\n s_min = " << s_min <<
"\n s_max = " << s_max
429 <<
"\n M1 = " << M1 <<
"\n M2 = " << M2
430 <<
"\n ml = " << ml << std::endl;
436 for ( j = max_j / 3; j < max_j; j++ ) {
437 s = s_min + ds * ( (double)j );
439 t_plus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
441 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
442 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
445 t_minus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
447 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
448 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
452 dt = ( t_plus - t_minus ) / ( (
double)max_k );
453 if ( fabs( dt ) < 0.00001 )
456 if ( dt <= ( -0.00001 ) ) {
458 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcScalarMaxProb(...)"
459 <<
"\n dt = " << dt <<
" < 0."
460 <<
"\n s = " << s <<
"\n s_min = " << s_min
461 <<
"\n s_max = " << s_max <<
"\n ds = " << ds
462 <<
"\n j = " << j <<
"\n t_plus = " << t_plus
463 <<
"\n t_minus = " << t_minus <<
"\n M1 = " << M1
464 <<
"\n M2 = " << M2 <<
"\n ml = " << ml
470 for ( k = 0; k < max_k; k++ ) {
471 t_for_s = t_minus + dt * ( (double)k );
473 if ( ( t_for_s > t_plus ) && ( t_for_s <= ( 1.0001 * t_plus ) ) ) {
476 if ( t_for_s > ( 1.0001 * t_plus ) ) {
478 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcScalarMaxProb(...)"
479 <<
"\n t_for_s = " << t_for_s
480 <<
" > t_plus = " << t_plus <<
" ! "
481 <<
"\n t_minus = " << t_minus <<
"\n dt = " << dt
482 <<
"\n k = " << k <<
"\n s = " << s
483 <<
"\n M1 = " << M1 <<
"\n M2 = " << M2
484 <<
"\n ml = " << ml << std::endl;
490 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s ) /
495 El1 = ( pow( M1, 2.0 ) + pow( ml, 2.0 ) - t_for_s ) /
498 El1 = 1.0000001 * ml;
500 El2 = ( s + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) /
503 El2 = 1.0000001 * ml;
507 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) );
508 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
511 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) +
512 2.0 * EV * El2 - t_for_s ) /
513 ( 2.0 * modV * modl2 );
514 if ( ( fabs( cosVellminus ) > 1.0 ) &&
515 ( fabs( cosVellminus ) <= 1.0001 ) ) {
520 cosVellminus = cosVellminus / fabs( cosVellminus );
522 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) {
523 cosVellminus = cosVellminus / fabs( cosVellminus );
542 if ( fabs( cosVellminus ) > 1.0001 ) {
544 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
545 <<
"\n |cos(theta)| = " << fabs( cosVellminus ) <<
" > 1"
546 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
547 <<
"\n s_min = " << s_min <<
"\n s_max = " << s_max
548 <<
"\n t_plus = " << t_plus <<
"\n t_minus = " << t_minus
549 <<
"\n dt = " << dt <<
"\n EV = " << EV
550 <<
"\n El2 = " << El2 <<
"\n modV = " << modV
551 <<
"\n modl2 = " << modl2 <<
"\n M2 = " << M2
552 <<
"\n ml = " << ml << std::endl;
556 double sin2Vellminus = 1.0 - pow( cosVellminus, 2.0 );
557 if ( ( sin2Vellminus < 0.0 ) && ( sin2Vellminus >= -0.0001 ) ) {
560 if ( sin2Vellminus <= -0.0001 ) {
562 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
563 <<
"\n cos^2(theta) = " << sin2Vellminus <<
" < -0.001"
564 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
565 <<
"\n s_min = " << s_min <<
"\n s_max = " << s_max
566 <<
"\n t_plus = " << t_plus <<
"\n t_minus = " << t_minus
567 <<
"\n dt = " << dt <<
"\n EV = " << EV
568 <<
"\n El2 = " << El2 <<
"\n modV = " << modV
569 <<
"\n modl2 = " << modl2 <<
"\n M2 = " << M2
570 <<
"\n ml = " << ml << std::endl;
575 p1.
set( M1, 0.0, 0.0, 0.0 );
576 p2.
set( EV, modV, 0.0, 0.0 );
577 k2.
set( El2, modl2 * cosVellminus,
578 -modl2 * sqrt( sin2Vellminus ), 0.0 );
610 scalar_part->
init( parnum, p1 );
616 listdaug[0] = mesnum;
621 amp.
init( parnum, 3, listdaug );
627 vect = root_part->
getDaug( 0 );
628 lep1 = root_part->
getDaug( 1 );
629 lep2 = root_part->
getDaug( 2 );
650 vect->
init( mesnum, p2 );
651 lep1->
init( l1num, k1 );
652 lep2->
init( l2num, k2 );
665 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch,
666 ias, CKM_A, CKM_lambda, CKM_barrho, CKM_bareta );
671 if ( nikmax > maxfoundprob ) {
672 maxfoundprob = nikmax;
694 if ( res_swch == 1 ) {
696 double t_plus, t_minus;
700 s = pow( 3.09688, 2.0 );
702 t_plus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
703 t_plus = t_plus + sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
704 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
707 t_minus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
709 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
710 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
713 dt = ( t_plus - t_minus ) / 1000.0;
716 for ( k = 0; k < 1000; k++ ) {
717 t_for_s = t_minus + dt * ( (double)k );
719 if ( ( t_for_s > t_plus ) && ( t_for_s <= ( 1.0001 * t_plus ) ) ) {
722 if ( t_for_s > ( 1.0001 * t_plus ) ) {
724 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
725 <<
"\n t_for_s = " << t_for_s <<
" > t_plus = " << t_plus
727 <<
"\n t_minus = " << t_minus <<
"\n dt = " << dt
728 <<
"\n k = " << k <<
"\n s = " << s
729 <<
"\n M1 = " << M1 <<
"\n M2 = " << M2
730 <<
"\n ml = " << ml << std::endl;
736 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s ) /
738 El2 = ( s + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) /
742 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) );
743 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
746 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 2.0 * EV * El2 -
748 ( 2.0 * modV * modl2 );
749 if ( ( fabs( cosVellminus ) > 1.0 ) &&
750 ( fabs( cosVellminus ) <= 1.0001 ) ) {
755 cosVellminus = cosVellminus / fabs( cosVellminus );
757 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) {
758 cosVellminus = cosVellminus / fabs( cosVellminus );
775 if ( fabs( cosVellminus ) > 1.0001 ) {
777 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
778 <<
"\n |cos(theta)| = " << fabs( cosVellminus ) <<
" > 1"
779 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
780 <<
"\n t_plus = " << t_plus <<
"\n t_minus = " << t_minus
781 <<
"\n dt = " << dt <<
"\n EV = " << EV
782 <<
"\n El2 = " << El2 <<
"\n modV = " << modV
783 <<
"\n modl2 = " << modl2 <<
"\n M2 = " << M2
784 <<
"\n ml = " << ml << std::endl;
788 double sin2Vellminus = 1.0 - pow( cosVellminus, 2.0 );
789 if ( ( sin2Vellminus < 0.0 ) && ( sin2Vellminus >= -0.0001 ) ) {
792 if ( sin2Vellminus <= -0.0001 ) {
794 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
795 <<
"\n cos^2(theta) = " << sin2Vellminus <<
" < -0.001"
796 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
797 <<
"\n t_plus = " << t_plus <<
"\n t_minus = " << t_minus
798 <<
"\n dt = " << dt <<
"\n EV = " << EV
799 <<
"\n El2 = " << El2 <<
"\n modV = " << modV
800 <<
"\n modl2 = " << modl2 <<
"\n M2 = " << M2
801 <<
"\n ml = " << ml << std::endl;
806 p1.
set( M1, 0.0, 0.0, 0.0 );
807 p2.
set( EV, modV, 0.0, 0.0 );
808 k2.
set( El2, modl2 * cosVellminus, -modl2 * sqrt( sin2Vellminus ),
841 scalar_part->
init( parnum, p1 );
847 listdaug[0] = mesnum;
852 amp.
init( parnum, 3, listdaug );
858 vect = root_part->
getDaug( 0 );
859 lep1 = root_part->
getDaug( 1 );
860 lep2 = root_part->
getDaug( 2 );
866 vect->
init( mesnum, p2 );
867 lep1->
init( l1num, k1 );
868 lep2->
init( l2num, k2 );
875 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch,
876 ias, CKM_A, CKM_lambda, CKM_barrho, CKM_bareta );
881 if ( nikmax > maxfoundprob ) {
882 maxfoundprob = nikmax;
900 if ( maxfoundprob <= 0.0 ) {
902 <<
"\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)"
903 <<
"\n maxfoundprob = " << maxfoundprob <<
" <0 or =0!"
904 <<
"\n res_swch = " << res_swch << std::endl;
909 <<
"\n maxfoundprob (...) = " << maxfoundprob << std::endl;
911 maxfoundprob *= 1.01;
virtual void init(EvtId part_n, const EvtVector4R &p4)=0