66 double mu,
int Nf,
int sr,
int res_swch,
67 int ias,
double Egamma_min,
double CKM_A,
68 double CKM_lambda,
double CKM_barrho,
69 double CKM_bareta,
double mumumass_min )
80 double M1 = parent->
mass();
96 double Relambda_qu, Imlambda_qu;
101 if ( charge1 == 0 || charge2 == 0 ) {
103 <<
"\n\n The function EvtbsTollGammaAmp::CalcAmp(...)"
104 <<
"\n Error in the leptonic charge getting!"
105 <<
"\n charge1 =" << charge1
106 <<
"\n charge2 =" << charge2 <<
"\n charge gamma ="
108 <<
"\n number of daughters =" << parent->
getNDaug() << std::endl;
115 lepPlus = ( charge1 > charge2 )
118 lepMinus = ( charge1 < charge2 )
130 if ( charge1 > charge2 ) {
143 double q2 = q.
mass2();
144 double p2 = p.
mass2();
145 double t = p_minus_p_1.
mass2();
146 double u = p_minus_p_2.
mass2();
149 double pk = 0.5 * ( p2 - q2 );
150 double p1k = 0.5 * ( pow( ml, 2.0 ) - u );
151 double p2k = 0.5 * ( pow( ml, 2.0 ) - t );
153 double hatq2 = q2 / ( M1 * M1 );
154 double Egam = 0.5 * M1 *
176 Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) +
177 pow( CKM_lambda, 2.0 ) *
178 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
179 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
180 Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq;
182 Vuq = CKM_lambda * unit1;
184 Vcq = unit1 - 0.5 * pow( CKM_lambda, 2.0 ) -
185 0.125 * pow( CKM_lambda, 4.0 ) * ( 1.0 + 4.0 * pow( CKM_A, 2.0 ) );
193 Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) *
194 ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
195 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
196 Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq;
198 Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) -
199 0.125 * pow( CKM_lambda, 4.0 ) );
203 0.5 * pow( CKM_A, 2.0 ) * pow( CKM_lambda, 5.0 ) *
204 ( 1.0 - 2.0 * ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
205 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ) ) );
210 <<
"\n\n The function EvtbsTollGammaISRFSRAmp::CalcAmp(..// 4-momentum of ell^+.)"
211 <<
"\n Error in the model set!"
212 <<
" mq = " << mq << std::endl;
216 Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda,
218 Vub = CKM_A * pow( CKM_lambda, 3.0 ) *
219 ( CKM_barrho * unit1 - CKM_bareta * uniti ) /
220 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
221 Vcb = unit1 * CKM_A * pow( CKM_lambda, 2.0 );
223 CKM_factor =
conj( Vtq ) * Vtb;
225 lambda_qu =
conj( Vuq ) * Vub /
227 Relambda_qu =
real( lambda_qu );
228 Imlambda_qu =
imag( lambda_qu );
230 lambda_qc =
conj( Vcq ) * Vcb /
239 if ( Egam < Egamma_min || ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ||
240 ( q2 <= mumumass_min * mumumass_min ) ) {
245 c9eff_b2q = unit1 * 0.0;
246 c9eff_barb2barq = unit1 * 0.0;
249 c1 = WilsCoeff->
C1( mu, Mw, Nf, ias );
250 c2 = WilsCoeff->
C2( mu, Mw, Nf, ias );
251 a1 = unit1 * ( c1 + c2 / 3.0 );
252 c7gam = WilsCoeff->
GetC7Eff( mu, Mw, mt, Nf, ias );
253 c9eff_b2q = WilsCoeff->
GetC9Eff( 0, res_swch, ias, Nf, q2, mb, mq, mc, mu,
254 mt, Mw, ml, Relambda_qu, Imlambda_qu );
255 c9eff_barb2barq = WilsCoeff->
GetC9Eff( 1, res_swch, ias, Nf, q2, mb, mq,
256 mc, mu, mt, Mw, ml, Relambda_qu,
267 if ( Egam < Egamma_min || ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ||
268 ( q2 <= mumumass_min * mumumass_min ) ) {
272 Fta_b2q = unit1 * 0.0;
273 Fta_barb2barq = unit1 * 0.0;
274 Ftv_b2q = unit1 * 0.0;
275 Ftv_barb2barq = unit1 * 0.0;
279 <<
"\n\n The function EvtbsTollGammaISRFSRAmp::CalcAmp(...)"
280 <<
"\n Leptonic decay constant fb is not uninitialized in this function!"
281 <<
" fb = " << fb << std::endl;
287 a1, lambda_qu, lambda_qc, Fv, Fa, Ftv_b2q,
292 a1, lambda_qu, lambda_qc, Fv, Fa,
293 Ftv_barb2barq, Fta_barb2barq );
343 EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, e_b2q, e_barb2barq,
347 a_b2q = c9eff_b2q * Fv + 2.0 * c7gam * Ftv_b2q * mb * M1 / q2;
348 a_barb2barq = c9eff_barb2barq * Fv +
349 2.0 * c7gam * Ftv_barb2barq * mb * M1 / q2;
351 b_b2q = ( c9eff_b2q * Fa + 2.0 * c7gam * Fta_b2q * mb * M1 / q2 ) * pk /
353 b_barb2barq = ( c9eff_barb2barq * Fa +
354 2.0 * c7gam * Fta_barb2barq * mb * M1 / q2 ) *
360 f_b2q = c10a * Fa * pk / ( M1 * M1 );
364 brammT = 0.5 * c10a * ml * fb *
365 ( 1.0 / p2k + 1.0 / p1k );
402 static const EvtIdSet bmesons{
"anti-B0",
"anti-B_s0" };
403 static const EvtIdSet bbarmesons{
"B0",
"B_s0" };
407 if ( bmesons.
contains( parentID ) ) {
456 for ( i = 0; i < 2; i++ ) {
461 E1 = T1.
cont2( epsG );
462 E2 = T2.
cont2( epsG );
463 E3 = ( epsG * hatp ) * brammS;
466 if ( Egam < Egamma_min ||
467 ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ||
468 ( q2 <= mumumass_min * mumumass_min ) ) {
469 CKM_factor = 0.0 * unit1;
475 ( lvc11 * E1 + lac11 * E2 + uniti * lsc11 * E3 +
476 uniti * ( ( ltc11.
cont2( hatp ) ) * epsG ) *
485 ( lvc12 * E1 + lac12 * E2 + uniti * lsc12 * E3 +
486 uniti * ( ( ltc12.
cont2( hatp ) ) * epsG ) *
495 ( lvc21 * E1 + lac21 * E2 + uniti * lsc21 * E3 +
496 uniti * ( ( ltc21.
cont2( hatp ) ) * epsG ) *
505 ( lvc22 * E1 + lac22 * E2 + uniti * lsc22 * E3 +
506 uniti * ( ( ltc22.
cont2( hatp ) ) * epsG ) *
514 if ( bbarmesons.
contains( parentID ) ) {
518 T1 = -a_barb2barq * unit1 *
522 T2 = -e_barb2barq * unit1 *
564 for ( i = 0; i < 2; i++ ) {
567 E1 = T1.
cont2( barepsG );
568 E2 = T2.
cont2( barepsG );
569 E3 = ( barepsG * hatp ) * brammS;
572 if ( Egam < Egamma_min ||
573 ( res_swch == 1 && q2 >= 9.199 && q2 <= 15.333 ) ||
574 ( q2 <= mumumass_min * mumumass_min ) ) {
575 CKM_factor = 0.0 * unit1;
580 ( lvc11 * E1 + lac11 * E2 +
582 uniti * ( ( ltc11.
cont2( hatp ) ) * barepsG ) *
586 ( lvc12 * E1 + lac12 * E2 +
588 uniti * ( ( ltc12.
cont2( hatp ) ) * barepsG ) *
592 ( lvc21 * E1 + lac21 * E2 +
594 uniti * ( ( ltc21.
cont2( hatp ) ) * barepsG ) *
598 ( lvc22 * E1 + lac22 * E2 +
600 uniti * ( ( ltc22.
cont2( hatp ) ) * barepsG ) *
606 <<
"\n\n The function Evtbs2llGammaISRFSRAmp::CalcAmp(...)"
607 <<
"\n Wrong B-meson number" << std::endl;
629 int Nf,
int sr,
int res_swch,
int ias,
double Egamma_min,
double CKM_A,
630 double CKM_lambda,
double CKM_barrho,
double CKM_bareta,
double mumumass_min )
632 double maxfoundprob = -100.0;
640 std::string(
"omega" ) ) );
653 double list_of_max_q2_points[5];
654 list_of_max_q2_points[0] = pow( 2.0 * ml, 2.0 );
655 list_of_max_q2_points[1] = pow( Mrho, 2.0 );
656 list_of_max_q2_points[2] = pow( Momega, 2.0 );
657 list_of_max_q2_points[3] = pow( Mphi, 2.0 );
658 list_of_max_q2_points[4] =
659 pow( M1, 2.0 ) - 2.0 * M1 * Egamma_min;
672 if ( Egamma_min > Mrho ) {
674 <<
"\n\n In the function Evtbs2llGammaISRFSRAmp::CalcMaxProb(...)"
675 <<
"\n Bad photon energy cut: Egamma_min > M_rho0 in the rest frame of B-meson."
676 <<
"\n Mrho = " << Mrho <<
"\n Egamma_min = " << Egamma_min
681 if ( Egamma_min <= 0 ) {
683 <<
"\n\n In the function Evtbs2llGammaISRFSRAmp::CalcMaxProb(...)"
684 <<
"\n Bad photon energy cut: Egamma_min <= 0 in the rest frame of B-meson."
685 <<
"\n Egamma_min = " << Egamma_min << std::endl;
689 if ( res_swch == 0 || res_swch == 1 ) {
691 for ( i_list = 0; i_list <= 4; i_list++ ) {
699 s = list_of_max_q2_points[i_list];
701 t_plus = pow( M1, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
702 t_plus = t_plus + sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
703 ( pow( M1, 2.0 ) - s );
706 t_minus = pow( M1, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
707 t_minus = t_minus - sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
708 ( pow( M1, 2.0 ) - s );
711 if ( fabs( t_plus - t_minus ) < 0.000001 )
715 double dt = ( t_plus - t_minus ) / ( (
double)max_ijk );
716 if ( fabs( dt ) < 0.00001 )
721 <<
"\n\n In the function EvtbsTollGammaISRFSRAmp::CalcScalarMaxProb(...)"
722 <<
"\n dt = " << dt <<
" < 0."
723 <<
"\n s = " << s <<
"\n t_plus = " << t_plus
724 <<
"\n t_minus = " << t_minus <<
"\n M1 = " << M1
725 <<
"\n ml = " << ml << std::endl;
729 for ( ijk = 0; ijk <= max_ijk; ijk++ ) {
730 t_for_s = t_minus + dt * ( (double)ijk );
734 Eg = ( pow( M1, 2.0 ) - s ) / ( 2.0 * M1 );
735 El2 = ( s + t_for_s - pow( ml, 2.0 ) ) /
739 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
742 cosBellminus = ( pow( ml, 2.0 ) + 2.0 * Eg * El2 - t_for_s ) /
743 ( 2.0 * Eg * modl2 );
745 if ( ( fabs( cosBellminus ) > 1.0 ) &&
746 ( fabs( cosBellminus ) <= 1.0001 ) ) {
748 <<
"\n Debug in the function EvtbsTollGammaISRFSRAmp::CalcMaxProb(...):"
749 <<
"\n cos(theta) = " << cosBellminus << std::endl;
750 cosBellminus = cosBellminus / fabs( cosBellminus );
752 if ( fabs( cosBellminus ) > 1.0001 ) {
754 <<
"\n\n In the function EvtbsTollGammaISRFSRAmp::CalcMaxProb(...)"
755 <<
"\n |cos(theta)| = " << fabs( cosBellminus ) <<
" > 1"
756 <<
"\n s = " << s <<
"\n t_for_s = " << t_for_s
757 <<
"\n t_plus = " << t_plus <<
"\n t_minus = " << t_minus
758 <<
"\n dt = " << dt <<
"\n Eg = " << Eg
759 <<
"\n El2 = " << El2 <<
"\n modl2 = " << modl2
760 <<
"\n ml = " << ml << std::endl;
765 p.
set( M1, 0.0, 0.0, 0.0 );
766 k.
set( Eg, Eg, 0.0, 0.0 );
767 p2.
set( El2, modl2 * cosBellminus,
768 -modl2 * sqrt( 1.0 - pow( cosBellminus, 2.0 ) ), 0.0 );
777 scalar_part->
init( parnum, p );
783 listdaug[0] = photnum;
788 amp.
init( parnum, 3, listdaug );
794 gamm = root_part->
getDaug( 0 );
795 lep1 = root_part->
getDaug( 1 );
796 lep2 = root_part->
getDaug( 2 );
802 gamm->
init( photnum, k );
803 lep1->
init( l1num, p1 );
804 lep2->
init( l2num, p2 );
811 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, sr,
812 res_swch, ias, Egamma_min, CKM_A, CKM_lambda,
813 CKM_barrho, CKM_bareta, mumumass_min );
818 if ( nikmax > maxfoundprob ) {
819 double maxfoundprob_old;
820 maxfoundprob_old = maxfoundprob;
821 maxfoundprob = nikmax;
823 <<
"\n maxfoundprob ( s =" << s <<
", t = " << t_for_s
824 <<
" ) = " << maxfoundprob
825 <<
"\n maxfoundprob_old = " << maxfoundprob_old
826 <<
"\n ijk =" << ijk << std::endl;
840 <<
"\n\n In the function Evtbs2llGammaISRFSRAmp::CalcMaxProb(...)"
841 <<
"\n Unexpected value of the variable res_swch !!!"
842 <<
"\n res_swch = " << res_swch << std::endl;
846 if ( maxfoundprob < 0.0 ) {
848 <<
"\n\n In the function Evtbs2llGammaISRFSRAmp::CalcMaxProb(...)"
849 <<
"\n maxfoundprob = " << maxfoundprob <<
" <0 or =0!"
850 <<
"\n mu =" << mu <<
" Nf =" << Nf <<
"\n sr =" << sr
851 <<
" res_swch =" << res_swch <<
" ias =" << ias
852 <<
"\n Egamma_min =" << Egamma_min <<
"\n CKM_A = " << CKM_A
853 <<
" CKM_lambda = " << CKM_lambda <<
"\n CKM_barrho = " << CKM_barrho
854 <<
" CKM_bareta = " << CKM_bareta << std::endl;
858 if ( maxfoundprob == 0.0 ) {
860 <<
"\n\n In the function Evtbs2llGammaISRFSRAmp::CalcMaxProb(...)"
861 <<
"\n maxfoundprob = " << maxfoundprob <<
" <0 or =0!"
862 <<
"\n mu =" << mu <<
" Nf =" << Nf <<
"\n sr =" << sr
863 <<
" res_swch =" << res_swch <<
" ias =" << ias
864 <<
"\n Egamma_min =" << Egamma_min <<
"\n CKM_A = " << CKM_A
865 <<
" CKM_lambda = " << CKM_lambda <<
"\n CKM_barrho = " << CKM_barrho
866 <<
" CKM_bareta = " << CKM_bareta << std::endl;
867 maxfoundprob = 0.00000001;
870 maxfoundprob *= 1.01;
873 <<
"\n **********************************************************************"
874 <<
"\n The function Evtbs2llGammaISRFSRAmp::CalcMaxProb(...) passed with arguments:"
875 <<
"\n mu =" << mu <<
" Nf =" << Nf <<
"\n sr =" << sr
876 <<
" res_swch =" << res_swch <<
" ias =" << ias
877 <<
"\n CKM_A = " << CKM_A <<
" CKM_lambda = " << CKM_lambda
878 <<
"\n Egamma_min =" << Egamma_min <<
"\n CKM_barrho = " << CKM_barrho
879 <<
" CKM_bareta = " << CKM_bareta
880 <<
"\n The distribution maximum maxfoundprob =" << maxfoundprob
881 <<
"\n **********************************************************************"
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
void init(EvtId part_n, double e, double px, double py, double pz)