71 double mu,
int Nf,
int res_swch,
int ias,
72 double CKM_A,
double CKM_lambda,
73 double CKM_barrho,
double CKM_bareta )
89 double q2 = q.
mass2();
91 double M1 = parent->
mass();
105 double Relambda_qu, Imlambda_qu;
138 iddaught ==
EvtPDL::getId( std::string(
"anti-K'_10" ) ) ) ) {
141 Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) +
142 pow( CKM_lambda, 2.0 ) *
143 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
144 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
145 Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq;
147 Vuq = CKM_lambda * unit1;
168 Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) *
169 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
170 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
171 Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq;
173 Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) -
174 0.125 * pow( CKM_lambda, 4.0 ) );
179 <<
"\n\n The function EvtbTosllVectorAmpNew::CalcAmp(...)"
180 <<
"\n Error in the model set!"
181 <<
" ms = " << ms << std::endl;
185 Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda,
187 Vub = CKM_A * pow( CKM_lambda, 3.0 ) *
188 ( CKM_barrho * unit1 - CKM_bareta * uniti ) /
189 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
191 CKM_factor =
conj( Vtq ) * Vtb;
193 lambda_qu =
conj( Vuq ) * Vub /
195 Relambda_qu =
real( lambda_qu );
196 Imlambda_qu =
imag( lambda_qu );
198 double a1,
a2, a0, v, t1, t2, t3;
202 q2,
a1,
a2, a0, v, t1, t2, t3 );
208 ms, mc, mu, mt, Mw, ml,
209 Relambda_qu, Imlambda_qu );
211 mb, ms, mc, mu, mt, Mw, ml,
212 Relambda_qu, Imlambda_qu );
270 double hats = q2 / pow( M1, 2 );
271 double hatM2 = M2 / M1;
272 double hatmb = mb / M1;
273 double hatms = ms / M1;
278 EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, c_b2q, c_barb2barq, e, f,
281 a_b2q = 2.0 * c9eff_b2q * v / ( 1.0 + hatM2 ) +
282 4.0 * ( hatmb + hatms ) * c7gam * t1 / hats;
283 a_barb2barq = 2.0 * c9eff_barb2barq * v / ( 1.0 + hatM2 ) +
284 4.0 * ( hatmb + hatms ) * c7gam * t1 / hats;
286 b_b2q = ( c9eff_b2q *
a1 +
287 2.0 * ( hatmb - hatms ) * ( 1.0 - hatM2 ) * c7gam * t2 / hats ) *
289 b_barb2barq = ( c9eff_barb2barq *
a1 + 2.0 * ( hatmb - hatms ) *
290 ( 1.0 - hatM2 ) * c7gam * t2 /
294 c_b2q = ( c9eff_b2q * ( 1.0 - hatM2 ) *
a2 +
295 2.0 * ( hatmb - hatms ) * ( 1.0 - pow( hatM2, 2 ) ) * c7gam * t2 /
297 2.0 * ( hatmb - hatms ) * c7gam * t3 ) /
298 ( 1 - pow( hatM2, 2 ) );
300 c_barb2barq = ( c9eff_barb2barq * ( 1.0 - hatM2 ) *
a2 +
301 2.0 * ( hatmb - hatms ) * ( 1.0 - pow( hatM2, 2 ) ) *
303 2.0 * ( hatmb - hatms ) * c7gam * t3 ) /
304 ( 1 - pow( hatM2, 2 ) );
306 e = 2.0 * c10a * v / ( 1 + hatM2 );
308 f = ( 1.0 + hatM2 ) * c10a *
a1;
310 g = c10a *
a2 / ( 1 + hatM2 );
312 h = ( ( 1.0 + hatM2 ) *
a1 - ( 1.0 - hatM2 ) *
a2 - 2.0 * hatM2 * a0 ) *
334 lepPlus = ( charge1 > charge2 ) ? parent->
getDaug( 1 ) : parent->
getDaug( 2 );
335 lepMinus = ( charge1 < charge2 ) ? parent->
getDaug( 1 )
349 EvtIdSet bmesons{
"B-",
"anti-B0",
"anti-B_s0",
"B_c-" };
350 EvtIdSet bbarmesons{
"B+",
"B0",
"B_s0",
"B_c+" };
354 if ( bmesons.
contains( parentID ) ) {
387 for ( i = 0; i < 3; i++ ) {
393 amp.
vertex( i, 0, 0, CKM_factor * ( lvc11 * E1 + lac11 * E2 ) );
394 amp.
vertex( i, 0, 1, CKM_factor * ( lvc12 * E1 + lac12 * E2 ) );
395 amp.
vertex( i, 1, 0, CKM_factor * ( lvc21 * E1 + lac21 * E2 ) );
396 amp.
vertex( i, 1, 1, CKM_factor * ( lvc22 * E1 + lac22 * E2 ) );
400 if ( bbarmesons.
contains( parentID ) ) {
404 T1 = a_barb2barq * unit1 *
407 c_barb2barq * uniti *
436 for ( i = 0; i < 3; i++ ) {
443 conj( CKM_factor ) * ( lvc11 * E3 + lac11 * E4 ) );
445 conj( CKM_factor ) * ( lvc12 * E3 + lac12 * E4 ) );
447 conj( CKM_factor ) * ( lvc21 * E3 + lac21 * E4 ) );
449 conj( CKM_factor ) * ( lvc22 * E3 + lac22 * E4 ) );
454 <<
"\n\n The function EvtbTosllVectorAmpNew::CalcAmp(...)"
455 <<
"\n Wrong B-meson number" << std::endl;
493 int Nf,
int res_swch,
int ias,
double CKM_A,
double CKM_lambda,
494 double CKM_barrho,
double CKM_bareta )
496 double maxfoundprob = -100.0;
503 if ( res_swch == 0 ) {
505 double s_min, t_for_s;
506 s_min = 4.0 * pow( ml, 2.0 );
508 ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - 2.0 * pow( ml, 2.0 ) );
511 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s_min ) /
513 El2 = ( s_min + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) /
517 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) );
518 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
521 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 2.0 * EV * El2 -
523 ( 2.0 * modV * modl2 );
524 if ( ( fabs( cosVellminus ) > 1.0 ) &&
525 ( fabs( cosVellminus ) <= 1.0001 ) ) {
530 cosVellminus = cosVellminus / fabs( cosVellminus );
532 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) {
533 cosVellminus = cosVellminus / fabs( cosVellminus );
535 <<
"\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):"
536 <<
"\n modV = " << modV <<
"\n modl2 = " << modl2
537 <<
"\n cos(theta) = " << cosVellminus
538 <<
"\n t_for_s = " << t_for_s <<
"\n s_min = " << s_min
539 <<
"\n EV = " << EV <<
"\n El2 = " << El2
540 <<
"\n M2 = " << M2 <<
"\n ml = " << ml
543 if ( fabs( cosVellminus ) > 1.0001 ) {
545 <<
"\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)"
546 <<
"\n |cos(theta)| = " << fabs( cosVellminus ) <<
" > 1"
547 <<
"\n s_min = " << s_min <<
"\n t_for_s = " << t_for_s
548 <<
"\n EV = " << EV <<
"\n El2 = " << El2
549 <<
"\n modV = " << modV <<
"\n modl2 = " << modl2
550 <<
"\n M2 = " << M2 <<
"\n ml = " << ml << std::endl;
555 p1.
set( M1, 0.0, 0.0, 0.0 );
556 p2.
set( EV, modV, 0.0, 0.0 );
557 k2.
set( El2, modl2 * cosVellminus,
558 -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 );
594 scalar_part->
init( parnum, p1 );
600 listdaug[0] = mesnum;
605 amp.
init( parnum, 3, listdaug );
611 vect = root_part->
getDaug( 0 );
612 lep1 = root_part->
getDaug( 1 );
613 lep2 = root_part->
getDaug( 2 );
619 vect->
init( mesnum, p2 );
620 lep1->
init( l1num, k1 );
621 lep2->
init( l2num, k2 );
628 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch, ias,
629 CKM_A, CKM_lambda, CKM_barrho, CKM_bareta );
642 if ( res_swch == 1 ) {
644 double t_plus, t_minus;
648 s = pow( 3.09688, 2.0 );
650 t_plus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
651 t_plus = t_plus + sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
652 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
655 t_minus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
657 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
658 sqrt(
lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
661 dt = ( t_plus - t_minus ) / 1000.0;
664 for ( k = 0; k <= 1000; k++ ) {
665 t_for_s = t_plus - dt * ( (double)k );
667 if ( ( t_for_s < t_minus ) && ( t_for_s >= ( 0.9999 * t_minus ) ) ) {
670 if ( t_for_s < ( 0.9999 * t_minus ) ) {
672 <<
"\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)"
673 <<
"\n t_for_s = " << t_for_s <<
" < t_minus = " << t_minus
675 <<
"\n t_plus = " << t_plus <<
"\n dt = " << dt
676 <<
"\n k = " << k <<
"\n s = " << s
677 <<
"\n M1 = " << M1 <<
"\n M2 = " << M2
678 <<
"\n ml = " << ml << std::endl;
684 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s ) /
686 El2 = ( s + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) /
690 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) );
691 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
694 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 2.0 * EV * El2 -
696 ( 2.0 * modV * modl2 );
697 if ( ( fabs( cosVellminus ) > 1.0 ) &&
698 ( fabs( cosVellminus ) <= 1.0001 ) ) {
703 cosVellminus = cosVellminus / fabs( cosVellminus );
705 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) {
706 cosVellminus = cosVellminus / fabs( cosVellminus );
708 <<
"\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):"
709 <<
"\n modV = " << modV <<
"\n modl2 = " << modl2
710 <<
"\n cos(theta) = " << cosVellminus
711 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
712 <<
"\n t_plus = " << t_plus
713 <<
"\n t_minus = " << t_minus <<
"\n dt = " << dt
714 <<
"\n EV = " << EV <<
"\n El2 = " << El2
715 <<
"\n M2 = " << M2 <<
"\n ml = " << ml
718 if ( fabs( cosVellminus ) > 1.0001 ) {
720 <<
"\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)"
721 <<
"\n |cos(theta)| = " << fabs( cosVellminus ) <<
" > 1"
722 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
723 <<
"\n EV = " << EV <<
"\n El2 = " << El2
724 <<
"\n modV = " << modV <<
"\n modl2 = " << modl2
725 <<
"\n M2 = " << M2 <<
"\n ml = " << ml
731 p1.
set( M1, 0.0, 0.0, 0.0 );
732 p2.
set( EV, modV, 0.0, 0.0 );
733 k2.
set( El2, modl2 * cosVellminus,
734 -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 );
770 scalar_part->
init( parnum, p1 );
776 listdaug[0] = mesnum;
781 amp.
init( parnum, 3, listdaug );
787 vect = root_part->
getDaug( 0 );
788 lep1 = root_part->
getDaug( 1 );
789 lep2 = root_part->
getDaug( 2 );
795 vect->
init( mesnum, p2 );
796 lep1->
init( l1num, k1 );
797 lep2->
init( l2num, k2 );
804 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch,
805 ias, CKM_A, CKM_lambda, CKM_barrho, CKM_bareta );
810 if ( nikmax > maxfoundprob ) {
811 maxfoundprob = nikmax;
814 <<
"\n maxfoundprob ( s =" << s <<
", t = " << t_for_s
815 <<
" ) = " << maxfoundprob <<
"\n k =" << katmax
829 if ( maxfoundprob <= 0.0 ) {
831 <<
"\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)"
832 <<
"\n maxfoundprob = " << maxfoundprob <<
" <0 or =0!"
833 <<
"\n res_swch = " << res_swch << std::endl;
838 <<
"\n maxfoundprob (...) = " << maxfoundprob << std::endl;
840 maxfoundprob *= 1.01;
virtual EvtVector4C epsParent(int i) const
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
virtual EvtDiracSpinor spParent(int) const
void makeDaughters(size_t ndaug, const EvtId *id)
void init(EvtId part_n, double e, double px, double py, double pz)
double CalcMaxProb(EvtId parnum, EvtId mesnum, EvtId l1num, EvtId l2num, EvtbTosllFFNew *formFactors, EvtbTosllWilsCoeffNLO *WilsCoeff, double mu, int Nf, int res_swch, int ias, double CKM_A, double CKM_lambda, double CKM_barrho, double CKM_bareta) override