62 int Nf,
int res_swch,
int ias,
double Egamma_min,
63 double CKM_A,
double CKM_lambda,
64 double CKM_barrho,
double CKM_bareta )
75 double M1 = parent->
mass();
91 double Relambda_qu, Imlambda_qu;
96 if ( charge1 == 0 || charge2 == 0 ) {
98 <<
"\n\n The function EvtbsTollGammaAmp::CalcAmp(...)"
99 <<
"\n Error in the leptonic charge getting!"
100 <<
"\n charge1 =" << charge1
101 <<
"\n charge2 =" << charge2 <<
"\n charge gamma ="
103 <<
"\n number of daughters =" << parent->
getNDaug() << std::endl;
110 lepPlus = ( charge1 > charge2 )
113 lepMinus = ( charge1 < charge2 )
125 if ( charge1 > charge2 ) {
138 double q2 = q.
mass2();
139 double p2 = p.
mass2();
140 double t = p_minus_p_1.
mass2();
141 double u = p_minus_p_2.
mass2();
144 double pk = 0.5 * ( p2 - q2 );
145 double p1k = 0.5 * ( pow( ml, 2.0 ) - u );
146 double p2k = 0.5 * ( pow( ml, 2.0 ) - t );
148 double hatq2 = q2 / ( M1 * M1 );
149 double Egam = 0.5 * M1 *
171 Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) +
172 pow( CKM_lambda, 2.0 ) *
173 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
174 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
175 Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq;
177 Vuq = CKM_lambda * unit1;
179 Vcq = unit1 - 0.5 * pow( CKM_lambda, 2.0 ) -
180 0.125 * pow( CKM_lambda, 4.0 ) * ( 1.0 + 4.0 * pow( CKM_A, 2.0 ) );
188 Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) *
189 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
190 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
191 Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq;
193 Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) -
194 0.125 * pow( CKM_lambda, 4.0 ) );
198 0.5 * pow( CKM_A, 2.0 ) * pow( CKM_lambda, 5.0 ) *
199 ( 1.0 - 2.0 * ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
200 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ) ) );
205 <<
"\n\n The function EvtbsTollGammaAmp::CalcAmp(..// 4-momentum of ell^+.)"
206 <<
"\n Error in the model set!"
207 <<
" mq = " << mq << std::endl;
211 Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda,
213 Vub = CKM_A * pow( CKM_lambda, 3.0 ) *
214 ( CKM_barrho * unit1 - CKM_bareta * uniti ) /
215 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
216 Vcb = unit1 * CKM_A * pow( CKM_lambda, 2.0 );
218 CKM_factor =
conj( Vtq ) * Vtb;
220 lambda_qu =
conj( Vuq ) * Vub /
222 Relambda_qu =
real( lambda_qu );
223 Imlambda_qu =
imag( lambda_qu );
225 lambda_qc =
conj( Vcq ) * Vcb /
234 if ( Egam < Egamma_min || ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ) {
239 c9eff_b2q = unit1 * 0.0;
240 c9eff_barb2barq = unit1 * 0.0;
243 c1 = WilsCoeff->
C1( mu, Mw, Nf, ias );
244 c2 = WilsCoeff->
C2( mu, Mw, Nf, ias );
245 a1 = unit1 * ( c1 + c2 / 3.0 );
246 c7gam = WilsCoeff->
GetC7Eff( mu, Mw, mt, Nf, ias );
247 c9eff_b2q = WilsCoeff->
GetC9Eff( 0, res_swch, ias, Nf, q2, mb, mq, mc, mu,
248 mt, Mw, ml, Relambda_qu, Imlambda_qu );
249 c9eff_barb2barq = WilsCoeff->
GetC9Eff( 1, res_swch, ias, Nf, q2, mb, mq,
250 mc, mu, mt, Mw, ml, Relambda_qu,
261 if ( Egam < Egamma_min || ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ) {
265 Fta_b2q = unit1 * 0.0;
266 Fta_barb2barq = unit1 * 0.0;
267 Ftv_b2q = unit1 * 0.0;
268 Ftv_barb2barq = unit1 * 0.0;
272 <<
"\n\n The function EvtbsTollGammaAmp::CalcAmp(...)"
273 <<
"\n Leptonic decay constant fb is not uninitialized in this function!"
274 <<
" fb = " << fb << std::endl;
280 a1, lambda_qu, lambda_qc, Fv, Fa, Ftv_b2q,
285 a1, lambda_qu, lambda_qc, Fv, Fa,
286 Ftv_barb2barq, Fta_barb2barq );
335 EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, e_b2q, e_barb2barq,
339 a_b2q = c9eff_b2q * Fv + 2.0 * c7gam * Ftv_b2q * mb * M1 / q2;
340 a_barb2barq = c9eff_barb2barq * Fv +
341 2.0 * c7gam * Ftv_barb2barq * mb * M1 / q2;
343 b_b2q = ( c9eff_b2q * Fa + 2.0 * c7gam * Fta_b2q * mb * M1 / q2 ) * pk /
345 b_barb2barq = ( c9eff_barb2barq * Fa +
346 2.0 * c7gam * Fta_barb2barq * mb * M1 / q2 ) *
352 f_b2q = c10a * Fa * pk / ( M1 * M1 );
356 brammT = 0.5 * c10a * ml * fb *
357 ( 1.0 / p2k + 1.0 / p1k );
378 static const EvtIdSet bmesons{
"anti-B0",
"anti-B_s0" };
379 static const EvtIdSet bbarmesons{
"B0",
"B_s0" };
383 if ( bmesons.
contains( parentID ) ) {
432 for ( i = 0; i < 2; i++ ) {
437 E1 = T1.cont2( epsG );
438 E2 = T2.cont2( epsG );
439 E3 = ( epsG * hatp ) * brammS;
442 if ( Egam < Egamma_min ||
443 ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ) {
444 CKM_factor = 0.0 * unit1;
450 ( lvc11 * E1 + lac11 * E2 + uniti * lsc11 * E3 +
451 uniti * ( ( ltc11.
cont2( hatp ) ) * epsG ) *
460 ( lvc12 * E1 + lac12 * E2 + uniti * lsc12 * E3 +
461 uniti * ( ( ltc12.
cont2( hatp ) ) * epsG ) *
470 ( lvc21 * E1 + lac21 * E2 + uniti * lsc21 * E3 +
471 uniti * ( ( ltc21.
cont2( hatp ) ) * epsG ) *
480 ( lvc22 * E1 + lac22 * E2 + uniti * lsc22 * E3 +
481 uniti * ( ( ltc22.
cont2( hatp ) ) * epsG ) *
489 if ( bbarmesons.
contains( parentID ) ) {
493 T1 = -a_barb2barq * unit1 *
497 T2 = -e_barb2barq * unit1 *
539 for ( i = 0; i < 2; i++ ) {
542 E1 = T1.cont2( barepsG );
543 E2 = T2.cont2( barepsG );
544 E3 = ( barepsG * hatp ) * brammS;
547 if ( Egam < Egamma_min ||
548 ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ) {
549 CKM_factor = 0.0 * unit1;
554 ( lvc11 * E1 + lac11 * E2 +
556 uniti * ( ( ltc11.
cont2( hatp ) ) * barepsG ) *
560 ( lvc12 * E1 + lac12 * E2 +
562 uniti * ( ( ltc12.
cont2( hatp ) ) * barepsG ) *
566 ( lvc21 * E1 + lac21 * E2 +
568 uniti * ( ( ltc21.
cont2( hatp ) ) * barepsG ) *
572 ( lvc22 * E1 + lac22 * E2 +
574 uniti * ( ( ltc22.
cont2( hatp ) ) * barepsG ) *
580 <<
"\n\n The function Evtbs2llGammaAmp::CalcAmp(...)"
581 <<
"\n Wrong B-meson number" << std::endl;
603 double mu,
int Nf,
int res_swch,
int ias,
604 double Egamma_min,
double CKM_A,
605 double CKM_lambda,
double CKM_barrho,
608 double maxfoundprob = -100.0;
617 std::string(
"omega" ) ) );
630 double list_of_max_q2_points[5];
631 list_of_max_q2_points[0] = pow( 2.0 * ml, 2.0 );
632 list_of_max_q2_points[1] = pow( Mrho, 2.0 );
633 list_of_max_q2_points[2] = pow( Momega, 2.0 );
634 list_of_max_q2_points[3] = pow( Mphi, 2.0 );
635 list_of_max_q2_points[4] =
636 pow( M1, 2.0 ) - 2.0 * M1 * Egamma_min;
649 if ( Egamma_min > Mrho ) {
651 <<
"\n\n In the function EvtbsTollGammaAmp::CalcScalarMaxProb(...)"
652 <<
"\n Bad photon energy cut: Egamma_min > M_rho0 in the rest frame of B-meson."
653 <<
"\n Mrho = " << Mrho <<
"\n Egamma_min = " << Egamma_min
658 if ( Egamma_min <= 0 ) {
660 <<
"\n\n In the function EvtbsTollGammaAmp::CalcScalarMaxProb(...)"
661 <<
"\n Bad photon energy cut: Egamma_min <= 0 in the rest frame of B-meson."
662 <<
"\n Egamma_min = " << Egamma_min << std::endl;
666 if ( res_swch == 0 || res_swch == 1 ) {
668 for ( i_list = 0; i_list <= 4; i_list++ ) {
676 s = list_of_max_q2_points[i_list];
678 t_plus = pow( M1, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
679 t_plus = t_plus + sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
680 ( pow( M1, 2.0 ) - s );
683 t_minus = pow( M1, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
684 t_minus = t_minus - sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
685 ( pow( M1, 2.0 ) - s );
688 if ( fabs( t_plus - t_minus ) < 0.000001 )
692 double dt = ( t_plus - t_minus ) / ( (
double)max_ijk );
693 if ( fabs( dt ) < 0.00001 )
698 <<
"\n\n In the function EvtbsTollGammaAmp::CalcScalarMaxProb(...)"
699 <<
"\n dt = " << dt <<
" < 0."
700 <<
"\n s = " << s <<
"\n t_plus = " << t_plus
701 <<
"\n t_minus = " << t_minus <<
"\n M1 = " << M1
702 <<
"\n ml = " << ml << std::endl;
706 for ( ijk = 0; ijk <= max_ijk; ijk++ ) {
707 t_for_s = t_minus + dt * ( (double)ijk );
711 Eg = ( pow( M1, 2.0 ) - s ) / ( 2.0 * M1 );
712 El2 = ( s + t_for_s - pow( ml, 2.0 ) ) /
716 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
719 cosBellminus = ( pow( ml, 2.0 ) + 2.0 * Eg * El2 - t_for_s ) /
720 ( 2.0 * Eg * modl2 );
722 if ( ( fabs( cosBellminus ) > 1.0 ) &&
723 ( fabs( cosBellminus ) <= 1.0001 ) ) {
725 <<
"\n Debug in the function EvtbsTollGammaAmp::CalcMaxProb(...):"
726 <<
"\n cos(theta) = " << cosBellminus << std::endl;
727 cosBellminus = cosBellminus / fabs( cosBellminus );
729 if ( fabs( cosBellminus ) > 1.0001 ) {
731 <<
"\n\n In the function EvtbsTollGammaAmp::CalcMaxProb(...)"
732 <<
"\n |cos(theta)| = " << fabs( cosBellminus ) <<
" > 1"
733 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
734 <<
"\n t_plus = " << t_plus <<
"\n t_minus = " << t_minus
735 <<
"\n dt = " << dt <<
"\n Eg = " << Eg
736 <<
"\n El2 = " << El2 <<
"\n modl2 = " << modl2
737 <<
"\n ml = " << ml << std::endl;
742 p.
set( M1, 0.0, 0.0, 0.0 );
743 k.
set( Eg, Eg, 0.0, 0.0 );
744 p2.
set( El2, modl2 * cosBellminus,
745 -modl2 * sqrt( 1.0 - pow( cosBellminus, 2.0 ) ), 0.0 );
754 scalar_part->
init( parnum, p );
760 listdaug[0] = photnum;
765 amp.
init( parnum, 3, listdaug );
771 gamm = root_part->
getDaug( 0 );
772 lep1 = root_part->
getDaug( 1 );
773 lep2 = root_part->
getDaug( 2 );
779 gamm->
init( photnum, k );
780 lep1->
init( l1num, p1 );
781 lep2->
init( l2num, p2 );
788 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf,
789 res_swch, ias, Egamma_min, CKM_A, CKM_lambda,
790 CKM_barrho, CKM_bareta );
795 if ( nikmax > maxfoundprob ) {
796 maxfoundprob = nikmax;
798 <<
"\n maxfoundprob ( s =" << s <<
", t = " << t_for_s
799 <<
" ) = " << maxfoundprob <<
"\n ijk =" << ijk
814 <<
"\n\n In the function EvtbsTollVectorAmpNew::CalcMaxProb(...)"
815 <<
"\n Unexpected value of the variable res_swch !!!"
816 <<
"\n res_swch = " << res_swch << std::endl;
820 if ( maxfoundprob <= 0.0 ) {
822 <<
"\n\n In the function EvtbsTollVectorAmpNew::CalcMaxProb(...)"
823 <<
"\n maxfoundprob = " << maxfoundprob <<
" <0 or =0!"
824 <<
"\n res_swch = " << res_swch << std::endl;
828 maxfoundprob *= 1.01;
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
void init(EvtId part_n, double e, double px, double py, double pz)