|
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. 71 double CKM_A, double CKM_lambda, double CKM_barrho, double CKM_bareta, 72 double ReA7, double ImA7, double ReA10, double ImA10 ) 80 EvtComplex A10 = ReA10 * unit1 + ImA10 * uniti; 91 double q2 = q. mass2(); 93 double M1 = parent-> mass(); 107 double Relambda_qu, Imlambda_qu; 135 Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) + 136 pow( CKM_lambda, 2.0 ) * 137 ( CKM_barrho * unit1 + CKM_bareta * uniti ) / 138 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ); 139 Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq; 141 Vuq = CKM_lambda * unit1; 162 Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) * 163 ( CKM_barrho * unit1 + CKM_bareta * uniti ) / 164 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ); 165 Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq; 167 Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) - 168 0.125 * pow( CKM_lambda, 4.0 ) ); 173 << "\n\n The function EvtbTosllScalarAmpNew::CalcAmp(...)" 174 << "\n Error in the model set!" 175 << " ms = " << ms << std::endl; 179 Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda, 181 Vub = CKM_A * pow( CKM_lambda, 3.0 ) * 182 ( CKM_barrho * unit1 - CKM_bareta * uniti ) / 183 sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ); 185 CKM_factor = conj( Vtq ) * Vtb; 187 lambda_qu = conj( Vuq ) * Vub / 189 Relambda_qu = real( lambda_qu ); 190 Imlambda_qu = imag( lambda_qu ); 203 ms, mc, mu, mt, Mw, ml, 204 Relambda_qu, Imlambda_qu ); 206 mb, ms, mc, mu, mt, Mw, ml, 207 Relambda_qu, Imlambda_qu ); 266 double hats = q2 / pow( M1, 2 ); 267 double hatM2 = M2 / M1; 268 double hatmb = mb / M1; 269 double hatms = ms / M1; 272 EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, c, d; 274 a_b2q = c9eff_b2q * fp - 275 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 ); 276 a_barb2barq = c9eff_barb2barq * fp - 277 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 ); 279 b_b2q = ( c9eff_b2q * ( f0 - fp ) + 280 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 ) ) * 281 ( 1 - pow( hatM2, 2.0 ) ) / hats; 282 b_barb2barq = ( c9eff_barb2barq * ( f0 - fp ) + 283 2.0 * c7gam * ( hatmb + hatms ) * ft / ( 1.0 + hatM2 ) ) * 284 ( 1 - pow( hatM2, 2.0 ) ) / hats; 288 d = c10a * ( 1.0 - pow( hatM2, 2 ) ) * ( f0 - fp ) / hats; 297 lepPlus = ( charge1 > charge2 ) ? parent-> getDaug( 1 ) : parent-> getDaug( 2 ); 298 lepMinus = ( charge1 < charge2 ) ? parent-> getDaug( 1 ) 310 EvtIdSet bmesons( "B-", "anti-B0", "anti-B_s0", "B_c-" ); 311 EvtIdSet bbarmesons( "B+", "B0", "B_s0", "B_c+" ); 315 if ( bmesons. contains( parentID ) ) { 319 T1 = a_b2q * hatP + b_b2q * hatq; 321 T2 = c * hatP + d * hatq; 341 amp. vertex( 0, 0, CKM_factor * ( lvc11 * T1 + lac11 * T2 ) ); 342 amp. vertex( 0, 1, CKM_factor * ( lvc12 * T1 + lac12 * T2 ) ); 343 amp. vertex( 1, 0, CKM_factor * ( lvc21 * T1 + lac21 * T2 ) ); 344 amp. vertex( 1, 1, CKM_factor * ( lvc22 * T1 + lac22 * T2 ) ); 347 if ( bbarmesons. contains( parentID ) ) { 351 T1 = a_barb2barq * hatP + b_barb2barq * hatq; 352 T2 = c * hatP + d * hatq; 372 amp. vertex( 0, 0, conj( CKM_factor ) * ( lvc11 * T1 + lac11 * T2 ) ); 373 amp. vertex( 0, 1, conj( CKM_factor ) * ( lvc12 * T1 + lac12 * T2 ) ); 374 amp. vertex( 1, 0, conj( CKM_factor ) * ( lvc21 * T1 + lac21 * T2 ) ); 375 amp. vertex( 1, 1, conj( CKM_factor ) * ( lvc22 * T1 + lac22 * T2 ) ); 380 << "\n\n The function EvtbTosllScalarAmpNew::CalcAmp(...)" 381 << "\n Wrong B-meson number" << std::endl; 400 int Nf, int res_swch, int ias, double CKM_A, double CKM_lambda, 401 double CKM_barrho, double CKM_bareta, double ReA7, double ImA7, 402 double ReA10, double ImA10 ) 404 double maxfoundprob = -100.0; 410 if ( res_swch == 0 ) { 413 double t_plus, t_minus; 418 s_min = 4.0 * pow( ml, 2.0 ); 419 s_max = pow( ( M1 - M2 ), 2.0 ); 422 ds = ( s_max - s_min ) / ( ( double)max_j ); 425 << "\n\n In the function EvtbTosllScalarAmpNew::CalcScalarMaxProb(...)" 426 << "\n ds = " << ds << " < 0." 427 << "\n s_min = " << s_min << "\n s_max = " << s_max 428 << "\n M1 = " << M1 << "\n M2 = " << M2 429 << "\n ml = " << ml << std::endl; 435 for ( j = max_j / 3; j < max_j; j++ ) { 436 s = s_min + ds * ( (double)j ); 438 t_plus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s; 440 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) * 441 sqrt( lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) ); 444 t_minus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s; 446 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) * 447 sqrt( lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) ); 451 dt = ( t_plus - t_minus ) / ( ( double)max_k ); 452 if ( fabs( dt ) < 0.00001 ) 455 if ( dt <= ( -0.00001 ) ) { 457 << "\n\n In the function EvtbTosllScalarAmpNew::CalcScalarMaxProb(...)" 458 << "\n dt = " << dt << " < 0." 459 << "\n s = " << s << "\n s_min = " << s_min 460 << "\n s_max = " << s_max << "\n ds = " << ds 461 << "\n j = " << j << "\n t_plus = " << t_plus 462 << "\n t_minus = " << t_minus << "\n M1 = " << M1 463 << "\n M2 = " << M2 << "\n ml = " << ml 469 for ( k = 0; k < max_k; k++ ) { 470 t_for_s = t_minus + dt * ( (double)k ); 472 if ( ( t_for_s > t_plus ) && ( t_for_s <= ( 1.0001 * t_plus ) ) ) { 475 if ( t_for_s > ( 1.0001 * t_plus ) ) { 477 << "\n\n In the function EvtbTosllScalarAmpNew::CalcScalarMaxProb(...)" 478 << "\n t_for_s = " << t_for_s 479 << " > t_plus = " << t_plus << " ! " 480 << "\n t_minus = " << t_minus << "\n dt = " << dt 481 << "\n k = " << k << "\n s = " << s 482 << "\n M1 = " << M1 << "\n M2 = " << M2 483 << "\n ml = " << ml << std::endl; 489 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s ) / 491 El2 = ( s + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) / 495 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) ); 496 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) ); 499 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 500 2.0 * EV * El2 - t_for_s ) / 501 ( 2.0 * modV * modl2 ); 502 if ( ( fabs( cosVellminus ) > 1.0 ) && 503 ( fabs( cosVellminus ) <= 1.0001 ) ) { 508 cosVellminus = cosVellminus / fabs( cosVellminus ); 510 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) { 511 cosVellminus = cosVellminus / fabs( cosVellminus ); 530 if ( fabs( cosVellminus ) > 1.0001 ) { 532 << "\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)" 533 << "\n |cos(theta)| = " << fabs( cosVellminus ) << " > 1" 534 << "\n s = " << s << "\n t_for_s = " << t_for_s 535 << "\n s_min = " << s_min << "\n s_max = " << s_max 536 << "\n t_plus = " << t_plus << "\n t_minus = " << t_minus 537 << "\n dt = " << dt << "\n EV = " << EV 538 << "\n El2 = " << El2 << "\n modV = " << modV 539 << "\n modl2 = " << modl2 << "\n M2 = " << M2 540 << "\n ml = " << ml << std::endl; 545 p1. set( M1, 0.0, 0.0, 0.0 ); 546 p2. set( EV, modV, 0.0, 0.0 ); 547 k2. set( El2, modl2 * cosVellminus, 548 -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 ); 580 scalar_part-> init( parnum, p1 ); 586 listdaug[0] = mesnum; 591 amp. init( parnum, 3, listdaug ); 597 vect = root_part-> getDaug( 0 ); 598 lep1 = root_part-> getDaug( 1 ); 599 lep2 = root_part-> getDaug( 2 ); 605 vect-> init( mesnum, p2 ); 606 lep1-> init( l1num, k1 ); 607 lep2-> init( l2num, k2 ); 614 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, 615 res_swch, ias, CKM_A, CKM_lambda, CKM_barrho, 616 CKM_bareta, ReA7, ImA7, ReA10, ImA10 ); 621 if ( nikmax > maxfoundprob ) { 622 maxfoundprob = nikmax; 644 if ( res_swch == 1 ) { 646 double t_plus, t_minus; 650 s = pow( 3.09688, 2.0 ); 652 t_plus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s; 653 t_plus = t_plus + sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) * 654 sqrt( lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) ); 657 t_minus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s; 659 sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) * 660 sqrt( lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) ); 663 dt = ( t_plus - t_minus ) / 1000.0; 666 for ( k = 0; k < 1000; k++ ) { 667 t_for_s = t_minus + dt * ( (double)k ); 669 if ( ( t_for_s > t_plus ) && ( t_for_s <= ( 1.0001 * t_plus ) ) ) { 672 if ( t_for_s > ( 1.0001 * t_plus ) ) { 674 << "\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)" 675 << "\n t_for_s = " << t_for_s << " > t_plus = " << t_plus 677 << "\n t_minus = " << t_minus << "\n dt = " << dt 678 << "\n k = " << k << "\n s = " << s 679 << "\n M1 = " << M1 << "\n M2 = " << M2 680 << "\n ml = " << ml << std::endl; 686 EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s ) / 688 El2 = ( s + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) / 692 modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) ); 693 modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) ); 696 cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 2.0 * EV * El2 - 698 ( 2.0 * modV * modl2 ); 699 if ( ( fabs( cosVellminus ) > 1.0 ) && 700 ( fabs( cosVellminus ) <= 1.0001 ) ) { 705 cosVellminus = cosVellminus / fabs( cosVellminus ); 707 if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) { 708 cosVellminus = cosVellminus / fabs( cosVellminus ); 725 if ( fabs( cosVellminus ) > 1.0001 ) { 727 << "\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)" 728 << "\n |cos(theta)| = " << fabs( cosVellminus ) << " > 1" 729 << "\n s = " << s << "\n t_for_s = " << t_for_s 730 << "\n t_plus = " << t_plus << "\n t_minus = " << t_minus 731 << "\n dt = " << dt << "\n EV = " << EV 732 << "\n El2 = " << El2 << "\n modV = " << modV 733 << "\n modl2 = " << modl2 << "\n M2 = " << M2 734 << "\n ml = " << ml << std::endl; 739 p1. set( M1, 0.0, 0.0, 0.0 ); 740 p2. set( EV, modV, 0.0, 0.0 ); 741 k2. set( El2, modl2 * cosVellminus, 742 -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 ); 774 scalar_part-> init( parnum, p1 ); 780 listdaug[0] = mesnum; 785 amp. init( parnum, 3, listdaug ); 791 vect = root_part-> getDaug( 0 ); 792 lep1 = root_part-> getDaug( 1 ); 793 lep2 = root_part-> getDaug( 2 ); 799 vect-> init( mesnum, p2 ); 800 lep1-> init( l1num, k1 ); 801 lep2-> init( l2num, k2 ); 808 CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch, 809 ias, CKM_A, CKM_lambda, CKM_barrho, CKM_bareta, ReA7, ImA7, 815 if ( nikmax > maxfoundprob ) { 816 maxfoundprob = nikmax; 834 if ( maxfoundprob <= 0.0 ) { 836 << "\n\n In the function EvtbTosllScalarAmpNew::CalcMaxProb(...)" 837 << "\n maxfoundprob = " << maxfoundprob << " <0 or =0!" 838 << "\n res_swch = " << res_swch << std::endl; 843 << "\n maxfoundprob (...) = " << maxfoundprob << std::endl; 845 maxfoundprob *= 1.01; 867 l = pow( a, 2.0 ) + pow( b, 2.0 ) + pow( c, 2.0 ) - 2.0 * a * b - 868 2.0 * a * c - 2.0 * b * c;
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
EvtComplex GetC7Eff(double mu, double Mw, double mt, int Nf, int ias)
virtual void getScalarFF(EvtId, EvtId, double, double &, double &, double &)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
Evt3Rank3C conj(const Evt3Rank3C &t2)
static double getMeanMass(EvtId i)
int getSpinStates() const
void makeDaughters(unsigned int ndaug, EvtId *id)
void set(int i, double d)
EvtSpinDensity getSpinDensity()
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)
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
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)
void setDiagonalSpinDensity()
EvtParticle * getDaug(int i)
EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
int contains(const EvtId id)
double lambda(double a, double b, double c) override
EvtVector4R getP4Restframe() const
void init(EvtId p, int ndaug, EvtId *daug)
double real(const EvtComplex &c)
double normalizedProb(const EvtSpinDensity &d)
|