|
EvtGen
2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
|
Go to the documentation of this file. 75 double CKM_A, double CKM_lambda, double CKM_barrho, double CKM_bareta, 76 double ReA7, double ImA7, double ReA10, double ImA10 ) 84 EvtComplex A10 = ReA10 * unit1 + ImA10 * uniti; 95 double q2 = q. mass2(); 97 double M1 = parent-> mass(); 111 double Relambda_qu, Imlambda_qu; 144 iddaught == EvtPDL::getId( std::string( "anti-K'_10" ) ) ) ) { 147 Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) + 148 pow( CKM_lambda, 2.0 ) * 149 ( CKM_barrho * unit1 + CKM_bareta * uniti ) / 150 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ); 151 Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq; 153 Vuq = CKM_lambda * unit1; 174 Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) * 175 ( CKM_barrho * unit1 + CKM_bareta * uniti ) / 176 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ); 177 Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq; 179 Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) - 180 0.125 * pow( CKM_lambda, 4.0 ) ); 185 << "\n\n The function EvtbTosllVectorAmpNew::CalcAmp(...)" 186 << "\n Error in the model set!" 187 << " ms = " << ms << std::endl; 191 Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda, 193 Vub = CKM_A * pow( CKM_lambda, 3.0 ) * 194 ( CKM_barrho * unit1 - CKM_bareta * uniti ) / 195 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ); 197 CKM_factor = conj( Vtq ) * Vtb; 199 lambda_qu = conj( Vuq ) * Vub / 201 Relambda_qu = real( lambda_qu ); 202 Imlambda_qu = imag( lambda_qu ); 204 double a1, a2, a0, v, t1, t2, t3; 208 q2, a1, a2, a0, v, t1, t2, t3 ); 215 ms, mc, mu, mt, Mw, ml, 216 Relambda_qu, Imlambda_qu ); 218 mb, ms, mc, mu, mt, Mw, ml, 219 Relambda_qu, Imlambda_qu ); 278 double hats = q2 / pow( M1, 2 ); 279 double hatM2 = M2 / M1; 280 double hatmb = mb / M1; 281 double hatms = ms / M1; 286 EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, c_b2q, c_barb2barq, e, f, 289 a_b2q = 2.0 * c9eff_b2q * v / ( 1.0 + hatM2 ) + 290 4.0 * ( hatmb + hatms ) * c7gam * t1 / hats; 291 a_barb2barq = 2.0 * c9eff_barb2barq * v / ( 1.0 + hatM2 ) + 292 4.0 * ( hatmb + hatms ) * c7gam * t1 / hats; 294 b_b2q = ( c9eff_b2q * a1 + 295 2.0 * ( hatmb - hatms ) * ( 1.0 - hatM2 ) * c7gam * t2 / hats ) * 297 b_barb2barq = ( c9eff_barb2barq * a1 + 2.0 * ( hatmb - hatms ) * 298 ( 1.0 - hatM2 ) * c7gam * t2 / 302 c_b2q = ( c9eff_b2q * ( 1.0 - hatM2 ) * a2 + 303 2.0 * ( hatmb - hatms ) * ( 1.0 - pow( hatM2, 2 ) ) * c7gam * t2 / 305 2.0 * ( hatmb - hatms ) * c7gam * t3 ) / 306 ( 1 - pow( hatM2, 2 ) ); 308 c_barb2barq = ( c9eff_barb2barq * ( 1.0 - hatM2 ) * a2 + 309 2.0 * ( hatmb - hatms ) * ( 1.0 - pow( hatM2, 2 ) ) * 311 2.0 * ( hatmb - hatms ) * c7gam * t3 ) / 312 ( 1 - pow( hatM2, 2 ) ); 314 e = 2.0 * c10a * v / ( 1 + hatM2 ); 316 f = ( 1.0 + hatM2 ) * c10a * a1; 318 g = c10a * a2 / ( 1 + hatM2 ); 320 h = ( ( 1.0 + hatM2 ) * a1 - ( 1.0 - hatM2 ) * a2 - 2.0 * hatM2 * a0 ) * 342 lepPlus = ( charge1 > charge2 ) ? parent-> getDaug( 1 ) : parent-> getDaug( 2 ); 343 lepMinus = ( charge1 < charge2 ) ? parent-> getDaug( 1 ) 359 EvtIdSet bmesons( "B-", "anti-B0", "anti-B_s0", "B_c-" ); 360 EvtIdSet bbarmesons( "B+", "B0", "B_s0", "B_c+" ); 364 if ( bmesons. contains( parentID ) ) { 397 for ( i = 0; i < 3; i++ ) { 403 amp. vertex( i, 0, 0, CKM_factor * ( lvc11 * E1 + lac11 * E2 ) ); 404 amp. vertex( i, 0, 1, CKM_factor * ( lvc12 * E1 + lac12 * E2 ) ); 405 amp. vertex( i, 1, 0, CKM_factor * ( lvc21 * E1 + lac21 * E2 ) ); 406 amp. vertex( i, 1, 1, CKM_factor * ( lvc22 * E1 + lac22 * E2 ) ); 410 if ( bbarmesons. contains( parentID ) ) { 414 T1 = a_barb2barq * unit1 * 417 c_barb2barq * uniti * 446 for ( i = 0; i < 3; i++ ) { 453 conj( CKM_factor ) * ( lvc11 * E3 + lac11 * E4 ) ); 455 conj( CKM_factor ) * ( lvc12 * E3 + lac12 * E4 ) ); 457 conj( CKM_factor ) * ( lvc21 * E3 + lac21 * E4 ) ); 459 conj( CKM_factor ) * ( lvc22 * E3 + lac22 * E4 ) ); 464 << "\n\n The function EvtbTosllVectorAmpNew::CalcAmp(...)" 465 << "\n Wrong B-meson number" << std::endl; 503 int Nf, int res_swch, int ias, double CKM_A, double CKM_lambda, 504 double CKM_barrho, double CKM_bareta, double ReA7, double ImA7, 505 double ReA10, double ImA10 ) 507 double maxfoundprob = -100.0; 514 if ( res_swch == 0 ) { 516 double s_min, t_for_s; 517 s_min = 4.0 * pow( ml, 2.0 ); 519 ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - 2.0 * pow( ml, 2.0 ) ); 522 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s_min ) / 524 El2 = ( s_min + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) / 528 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) ); 529 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) ); 532 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 2.0 * EV * El2 - 534 ( 2.0 * modV * modl2 ); 535 if ( ( fabs( cosVellminus ) > 1.0 ) && 536 ( fabs( cosVellminus ) <= 1.0001 ) ) { 541 cosVellminus = cosVellminus / fabs( cosVellminus ); 543 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) { 544 cosVellminus = cosVellminus / fabs( cosVellminus ); 546 << "\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):" 547 << "\n modV = " << modV << "\n modl2 = " << modl2 548 << "\n cos(theta) = " << cosVellminus 549 << "\n t_for_s = " << t_for_s << "\n s_min = " << s_min 550 << "\n EV = " << EV << "\n El2 = " << El2 551 << "\n M2 = " << M2 << "\n ml = " << ml 554 if ( fabs( cosVellminus ) > 1.0001 ) { 556 << "\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)" 557 << "\n |cos(theta)| = " << fabs( cosVellminus ) << " > 1" 558 << "\n s_min = " << s_min << "\n t_for_s = " << t_for_s 559 << "\n EV = " << EV << "\n El2 = " << El2 560 << "\n modV = " << modV << "\n modl2 = " << modl2 561 << "\n M2 = " << M2 << "\n ml = " << ml << std::endl; 566 p1. set( M1, 0.0, 0.0, 0.0 ); 567 p2. set( EV, modV, 0.0, 0.0 ); 568 k2. set( El2, modl2 * cosVellminus, 569 -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 ); 605 scalar_part-> init( parnum, p1 ); 611 listdaug[0] = mesnum; 616 amp. init( parnum, 3, listdaug ); 622 vect = root_part-> getDaug( 0 ); 623 lep1 = root_part-> getDaug( 1 ); 624 lep2 = root_part-> getDaug( 2 ); 630 vect-> init( mesnum, p2 ); 631 lep1-> init( l1num, k1 ); 632 lep2-> init( l2num, k2 ); 639 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch, ias, 640 CKM_A, CKM_lambda, CKM_barrho, CKM_bareta, ReA7, ImA7, ReA10, 654 if ( res_swch == 1 ) { 656 double t_plus, t_minus; 660 s = pow( 3.09688, 2.0 ); 662 t_plus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s; 663 t_plus = t_plus + sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) * 664 sqrt( lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) ); 667 t_minus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s; 669 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) * 670 sqrt( lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) ); 673 dt = ( t_plus - t_minus ) / 1000.0; 676 for ( k = 0; k <= 1000; k++ ) { 677 t_for_s = t_plus - dt * ( (double)k ); 679 if ( ( t_for_s < t_minus ) && ( t_for_s >= ( 0.9999 * t_minus ) ) ) { 682 if ( t_for_s < ( 0.9999 * t_minus ) ) { 684 << "\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)" 685 << "\n t_for_s = " << t_for_s << " < t_minus = " << t_minus 687 << "\n t_plus = " << t_plus << "\n dt = " << dt 688 << "\n k = " << k << "\n s = " << s 689 << "\n M1 = " << M1 << "\n M2 = " << M2 690 << "\n ml = " << ml << std::endl; 696 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s ) / 698 El2 = ( s + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) / 702 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) ); 703 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) ); 706 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 2.0 * EV * El2 - 708 ( 2.0 * modV * modl2 ); 709 if ( ( fabs( cosVellminus ) > 1.0 ) && 710 ( fabs( cosVellminus ) <= 1.0001 ) ) { 715 cosVellminus = cosVellminus / fabs( cosVellminus ); 717 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) { 718 cosVellminus = cosVellminus / fabs( cosVellminus ); 720 << "\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):" 721 << "\n modV = " << modV << "\n modl2 = " << modl2 722 << "\n cos(theta) = " << cosVellminus 723 << "\n s = " << s << "\n t_for_s = " << t_for_s 724 << "\n t_plus = " << t_plus 725 << "\n t_minus = " << t_minus << "\n dt = " << dt 726 << "\n EV = " << EV << "\n El2 = " << El2 727 << "\n M2 = " << M2 << "\n ml = " << ml 730 if ( fabs( cosVellminus ) > 1.0001 ) { 732 << "\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)" 733 << "\n |cos(theta)| = " << fabs( cosVellminus ) << " > 1" 734 << "\n s = " << s << "\n t_for_s = " << t_for_s 735 << "\n EV = " << EV << "\n El2 = " << El2 736 << "\n modV = " << modV << "\n modl2 = " << modl2 737 << "\n M2 = " << M2 << "\n ml = " << ml 743 p1. set( M1, 0.0, 0.0, 0.0 ); 744 p2. set( EV, modV, 0.0, 0.0 ); 745 k2. set( El2, modl2 * cosVellminus, 746 -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 ); 782 scalar_part-> init( parnum, p1 ); 788 listdaug[0] = mesnum; 793 amp. init( parnum, 3, listdaug ); 799 vect = root_part-> getDaug( 0 ); 800 lep1 = root_part-> getDaug( 1 ); 801 lep2 = root_part-> getDaug( 2 ); 807 vect-> init( mesnum, p2 ); 808 lep1-> init( l1num, k1 ); 809 lep2-> init( l2num, k2 ); 816 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch, 817 ias, CKM_A, CKM_lambda, CKM_barrho, CKM_bareta, ReA7, ImA7, 823 if ( nikmax > maxfoundprob ) { 824 maxfoundprob = nikmax; 827 << "\n maxfoundprob ( s =" << s << ", t = " << t_for_s 828 << " ) = " << maxfoundprob << "\n k =" << katmax 842 if ( maxfoundprob <= 0.0 ) { 844 << "\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)" 845 << "\n maxfoundprob = " << maxfoundprob << " <0 or =0!" 846 << "\n res_swch = " << res_swch << std::endl; 851 << "\n maxfoundprob (...) = " << maxfoundprob << std::endl; 853 maxfoundprob *= 1.01; 878 l = pow( a, 2.0 ) + pow( b, 2.0 ) + pow( c, 2.0 ) - 2.0 * a * b - 879 2.0 * a * c - 2.0 * b * c; EvtTensor3C directProd(const EvtVector3C &c1, const EvtVector3C &c2)
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, double ReA7, double ImA7, double ReA10, double ImA10) override
static const EvtTensor4C & g()
EvtComplex GetC7Eff(double mu, double Mw, double mt, int Nf, int ias)
EvtTensor4C dual(const EvtTensor4C &t2)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
Evt3Rank3C conj(const Evt3Rank3C &t2)
EvtVector4C cont2(const EvtVector4C &v4) const
static double getMeanMass(EvtId i)
int getSpinStates() const
void makeDaughters(unsigned int ndaug, EvtId *id)
void set(int i, double d)
EvtSpinDensity getSpinDensity()
double lambda(double a, double b, double c) override
virtual EvtDiracSpinor spParent(int) const
EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void vertex(const EvtComplex &)
virtual double getQuarkMass(int)
EvtComplex GetC9Eff(int decay_id, int res_swch, int ias, int Nf, double q2, double m2, double md, double mc, double mu, double mt, double Mw, double ml, double Relambda_qu, double Imlambda_qu)
virtual EvtVector4C epsParent(int i) const
EvtComplex GetC10Eff(double mt, double Mw)
static EvtId getId(const std::string &name)
const EvtVector4R & getP4() const
double imag(const EvtComplex &c)
void init(EvtId part_n, double e, double px, double py, double pz)
virtual void getVectorFF(EvtId, EvtId, double, double &, double &, double &, double &, double &, double &, double &)
void setDiagonalSpinDensity()
EvtParticle * getDaug(int i)
EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
int contains(const EvtId id)
EvtVector4R getP4Restframe() const
void init(EvtId p, int ndaug, EvtId *daug)
double real(const EvtComplex &c)
double normalizedProb(const EvtSpinDensity &d)
void CalcAmp(EvtParticle *parent, EvtAmp &, 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, double ReA7, double ImA7, double ReA10, double ImA10) override
|