|
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. 59 double q2 = ( q. mass2() ); 61 double f1v, f1a, f2v, f2a; 65 q2, m_meson, &f1v, &f1a, &f2v, &f2a ); 68 p4b. set( parent-> mass(), 0.0, 0.0, 0.0 ); 159 if ( l_num == EM || l_num == MUM || l_num == TAUM ) { 165 if ( l_num == EP || l_num == MUP || l_num == TAUP ) { 172 << "Wrong lepton number" << endl; 176 amp. vertex( 0, 0, 0, l1. cont( temp_00_term1 + temp_00_term2 ) ); 177 amp. vertex( 0, 0, 1, l2. cont( temp_00_term1 + temp_00_term2 ) ); 179 amp. vertex( 0, 1, 0, l1. cont( temp_01_term1 + temp_01_term2 ) ); 180 amp. vertex( 0, 1, 1, l2. cont( temp_01_term1 + temp_01_term2 ) ); 182 amp. vertex( 1, 0, 0, l1. cont( temp_10_term1 + temp_10_term2 ) ); 183 amp. vertex( 1, 0, 1, l2. cont( temp_10_term1 + temp_10_term2 ) ); 185 amp. vertex( 1, 1, 0, l1. cont( temp_11_term1 + temp_11_term2 ) ); 186 amp. vertex( 1, 1, 1, l2. cont( temp_11_term1 + temp_11_term2 ) ); 225 dirac_part-> init( parent, p_init ); 234 listdaug[0] = baryon; 235 listdaug[1] = lepton; 236 listdaug[2] = nudaug; 238 amp. init( parent, 3, listdaug ); 241 daughter = root_part-> getDaug( 0 ); 243 trino = root_part-> getDaug( 2 ); 257 double m = root_part-> mass(); 262 double q2, elepton, plepton; 264 double erho, prho, costl; 266 double maxfoundprob = 0.0; 270 for ( massiter = 0; massiter < 3; massiter++ ) { 274 if ( massiter == 1 ) { 277 if ( massiter == 2 ) { 281 q2max = ( m - mass[0] ) * ( m - mass[0] ); 285 for ( i = 0; i < 25; i++ ) { 286 q2 = ( ( i + 0.5 ) * q2max ) / 25.0; 288 erho = ( m * m + mass[0] * mass[0] - q2 ) / ( 2.0 * m ); 290 prho = sqrt( erho * erho - mass[0] * mass[0] ); 292 p4baryon. set( erho, 0.0, 0.0, -1.0 * prho ); 293 p4w. set( m - erho, 0.0, 0.0, prho ); 296 elepton = ( q2 + mass[1] * mass[1] ) / ( 2.0 * sqrt( q2 ) ); 297 plepton = sqrt( elepton * elepton - mass[1] * mass[1] ); 301 for ( j = 0; j < 3; j++ ) { 302 costl = 0.99 * ( j - 1.0 ); 306 p4lepton. set( elepton, 0.0, plepton * sqrt( 1.0 - costl * costl ), 308 p4nu. set( plepton, 0.0, 309 -1.0 * plepton * sqrt( 1.0 - costl * costl ), 310 -1.0 * plepton * costl ); 312 EvtVector4R boost( ( m - erho ), 0.0, 0.0, 1.0 * prho ); 313 p4lepton = boostTo( p4lepton, boost ); 318 daughter-> init( baryon, p4baryon ); 319 lep-> init( lepton, p4lepton ); 320 trino-> init( nudaug, p4nu ); 322 CalcAmp( root_part, amp, FormFactors, r00, r01, r10, r11 ); 336 double a = probctl[1]; 337 double b = 0.5 * ( probctl[2] - probctl[0] ); 338 double c = 0.5 * ( probctl[2] + probctl[0] ) - probctl[1]; 341 if ( probctl[1] > prob ) 343 if ( probctl[2] > prob ) 346 if ( fabs( c ) > 1e-20 ) { 347 double ctlx = -0.5 * b / c; 348 if ( fabs( ctlx ) < 1.0 ) { 349 double probtmp = a + b * ctlx + c * ctlx * ctlx; 350 if ( probtmp > prob ) 360 if ( prob > maxfoundprob ) { 422 rho. set( 0, 0, r00 ); 423 rho. set( 0, 1, r01 ); 425 rho. set( 1, 0, r10 ); 426 rho. set( 1, 1, r11 ); 429 double pmag = vector4P. d3mag(); 430 double cosTheta = vector4P. get( 3 ) / pmag; 432 double theta = acos( cosTheta ); 433 double phi = atan2( vector4P. get( 2 ), vector4P. get( 1 ) ); 440 p4b. set( parent-> mass(), 0.0, 0.0, 0.0 ); 448 double q2 = q. mass2(); 460 double f1, f2, f3, g1, g2, g3; 461 FormFactors-> getdiracff( par_num, bar_num, q2, baryonmass, &f1, &f2, 462 &f3, &g1, &g2, &g3 ); 464 const double form_fact[6] = {f1, f2, f3, g1, g2, g3}; 469 if ( l_num == EM || l_num == MUM || l_num == TAUM ) { 475 } else if ( l_num == EP || l_num == MUP || l_num == TAUP ) { 497 if ( ( par_num == LAMB && bar_num == LAMCP ) || 498 ( par_num == LAMBB && bar_num == LAMCM ) || 499 ( par_num == LAMB && bar_num == PRO ) || 500 ( par_num == LAMBB && bar_num == PROB ) || 501 ( par_num == LAMB && bar_num == N1440 ) || 502 ( par_num == LAMBB && bar_num == N1440B ) || 503 ( par_num == LAMB && bar_num == N1710 ) || 504 ( par_num == LAMBB && bar_num == N1710B ) 508 if ( bar_num == LAMCP || bar_num == PRO || bar_num == N1440 || 511 else if ( bar_num == LAMCM || bar_num == PROB || 512 bar_num == N1440B || bar_num == N1710B ) 516 parent-> sp( 0 ), p4b, p4daught, form_fact, 519 parent-> sp( 1 ), p4b, p4daught, form_fact, 522 parent-> sp( 0 ), p4b, p4daught, form_fact, 525 parent-> sp( 1 ), p4b, p4daught, form_fact, 530 else if ( ( par_num == LAMB && bar_num == LAMC1P ) || 531 ( par_num == LAMBB && bar_num == LAMC1M ) || 532 ( par_num == LAMB && bar_num == N1535 ) || 533 ( par_num == LAMBB && bar_num == N1535B ) || 534 ( par_num == LAMB && bar_num == N1650 ) || 535 ( par_num == LAMBB && bar_num == N1650B ) ) { 537 if ( bar_num == LAMC1P || bar_num == N1535 || bar_num == N1650 ) 539 else if ( bar_num == LAMC1M || bar_num == N1535B || bar_num == N1650B ) 544 p4b, p4daught, form_fact, pflag ); 547 p4b, p4daught, form_fact, pflag ); 550 p4b, p4daught, form_fact, pflag ); 553 p4b, p4daught, form_fact, pflag ); 559 << "Dirac semilep. baryon current " 560 << "not implemented for this decay sequence." << std::endl; 564 amp. vertex( 0, 0, 0, l1 * b11 ); 565 amp. vertex( 0, 0, 1, l2 * b11 ); 567 amp. vertex( 1, 0, 0, l1 * b21 ); 568 amp. vertex( 1, 0, 1, l2 * b21 ); 570 amp. vertex( 0, 1, 0, l1 * b12 ); 571 amp. vertex( 0, 1, 1, l2 * b12 ); 573 amp. vertex( 1, 1, 0, l1 * b22 ); 574 amp. vertex( 1, 1, 1, l2 * b22 ); 583 double f1, f2, f3, f4, g1, g2, g3, g4; 584 FormFactors-> getraritaff( par_num, bar_num, q2, baryonmass, &f1, &f2, 585 &f3, &f4, &g1, &g2, &g3, &g4 ); 587 const double form_fact[8] = {f1, f2, f3, f4, g1, g2, g3, g4}; 591 EvtVector4C b11, b12, b21, b22, b13, b23, b14, b24, l1, l2; 594 if ( l_num == EM || l_num == MUM || l_num == TAUM ) { 600 } else if ( l_num == EP || l_num == MUP || l_num == TAUP ) { 617 if ( ( par_num == LAMB && bar_num == LAMC2P ) || 618 ( par_num == LAMB && bar_num == N1720 ) || 619 ( par_num == LAMB && bar_num == N1520 ) || 620 ( par_num == LAMB && bar_num == N1700 ) || 621 ( par_num == LAMB && bar_num == N1875 ) || 622 ( par_num == LAMB && bar_num == N1900 ) ) { 625 } else if ( ( par_num == LAMBB && bar_num == LAMC2M ) || 626 ( par_num == LAMBB && bar_num == N1520B ) || 627 ( par_num == LAMBB && bar_num == N1700B ) || 628 ( par_num == LAMBB && bar_num == N1875B ) ) { 633 else if ( ( par_num == LAMBB && bar_num == N1720B ) || 634 ( par_num == LAMBB && bar_num == N1900B ) ) { 638 << "Rarita-Schwinger semilep. baryon current " 639 << "not implemented for this decay sequence." << std::endl; 645 parent-> sp( 0 ), p4b, p4daught, 648 parent-> sp( 1 ), p4b, p4daught, 652 parent-> sp( 0 ), p4b, p4daught, 655 parent-> sp( 1 ), p4b, p4daught, 659 parent-> sp( 0 ), p4b, p4daught, 662 parent-> sp( 1 ), p4b, p4daught, 666 parent-> sp( 0 ), p4b, p4daught, 669 parent-> sp( 1 ), p4b, p4daught, 672 amp. vertex( 0, 0, 0, l1 * b11 ); 673 amp. vertex( 0, 0, 1, l2 * b11 ); 675 amp. vertex( 1, 0, 0, l1 * b21 ); 676 amp. vertex( 1, 0, 1, l2 * b21 ); 678 amp. vertex( 0, 1, 0, l1 * b12 ); 679 amp. vertex( 0, 1, 1, l2 * b12 ); 681 amp. vertex( 1, 1, 0, l1 * b22 ); 682 amp. vertex( 1, 1, 1, l2 * b22 ); 684 amp. vertex( 0, 2, 0, l1 * b13 ); 685 amp. vertex( 0, 2, 1, l2 * b13 ); 687 amp. vertex( 1, 2, 0, l1 * b23 ); 688 amp. vertex( 1, 2, 1, l2 * b23 ); 690 amp. vertex( 0, 3, 0, l1 * b14 ); 691 amp. vertex( 0, 3, 1, l2 * b14 ); 693 amp. vertex( 1, 3, 0, l1 * b24 ); 694 amp. vertex( 1, 3, 1, l2 * b24 ); 702 const double* ff, int pflag ) 726 else if ( pflag == 3 ) { 756 EvtVector4C current = ( ff[0] * t[0] + ff[1] * t[1] + ff[2] * t[2] - 757 ff[3] * t[3] - ff[4] * t[4] - ff[5] * t[5] ); 785 else if ( pflag == 2 ) { 795 id. setdiag( 1.0, 1.0, 1.0, 1.0 ); 798 for ( int i = 0; i < 4; i++ ) { 803 for ( int i = 0; i < 4; i++ ) { 814 ( parent / parent. mass() ); 819 ( daught / daught. mass() ); 822 t[3] = cg0 * ( id.cont2( v1 ) ); 832 ( parent / parent. mass() ); 838 ( daught / daught. mass() ); 841 t[7] = cg5 * ( id.cont2( v2 ) ); 844 EvtVector4C current = ( ff[0] * t[0] + ff[1] * t[1] + ff[2] * t[2] + 845 ff[3] * t[3] - ff[4] * t[4] - ff[5] * t[5] - 846 ff[6] * t[6] - ff[7] * t[7] );
virtual EvtRaritaSchwinger spRSParent(int) const
void setdiag(double t00, double t11, double t22, double t33)
virtual void getraritaff(EvtId parent, EvtId daught, double q2, double mass, double *f1, double *f2, double *f3, double *f4, double *g1, double *g2, double *g3, double *g4)=0
virtual EvtDiracSpinor sp(int) const
static EvtSpinType::spintype getSpinType(EvtId i)
EvtVector4R getP4Lab() const
void set(int, const EvtComplex &)
static const EvtGammaMatrix & g0()
EvtComplex EvtLeptonPCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
void set_spinor(int i, const EvtComplex &sp)
int getSpinStates() const
void makeDaughters(unsigned int ndaug, EvtId *id)
void set(int i, double d)
EvtRaritaSchwinger boostTo(const EvtRaritaSchwinger &rs, const EvtVector4R p4)
void init(EvtId part_n, const EvtVector4R &p4) override
EvtSpinDensity getSpinDensity()
virtual EvtDiracSpinor spParent(int) const
EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void vertex(const EvtComplex &)
static const EvtGammaMatrix & g1()
double CalcMaxProb(EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtSemiLeptonicFF *FormFactors, EvtComplex r00, EvtComplex r01, EvtComplex r10, EvtComplex r11)
static double getMinMass(EvtId i)
static EvtId getId(const std::string &name)
const EvtVector4R & getP4() const
static const EvtGammaMatrix & g2()
void setSpinDensityForwardHelicityBasis(const EvtSpinDensity &rho)
static double getWidth(EvtId i)
static const EvtGammaMatrix & g3()
virtual EvtDiracSpinor spParentNeutrino() const
EvtVector4C EvtBaryonVARaritaCurrent(const EvtRaritaSchwinger &Bf_vect, const EvtDiracSpinor &Bi, EvtVector4R parent, EvtVector4R daught, const double *ff, int pflag)
EvtVector4C EvtLeptonVACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void setDiagonalSpinDensity()
EvtParticle * getDaug(int i)
void CalcAmp(EvtParticle *parent, EvtAmp &, EvtSemiLeptonicFF *FormFactors) override
EvtVector4C getVector(int i) const
EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void set(int i, int j, const EvtComplex &rhoij)
EvtDiracSpinor getSpinor(int i) const
EvtVector4C EvtBaryonVACurrent(const EvtDiracSpinor &Bf, const EvtDiracSpinor &Bi, EvtVector4R parent, EvtVector4R daught, const double *ff, int pflag)
EvtComplex cont(const EvtVector4C &v4) const
static double getMass(EvtId i)
EvtComplex EvtLeptonSCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void init(EvtId p, int ndaug, EvtId *daug)
static const EvtGammaMatrix & g5()
virtual void getbaryonff(EvtId parent, EvtId daught, double t, double m_meson, double *f1v, double *f1a, double *f2v, double *f2a)=0
static double getMaxMass(EvtId i)
virtual void getdiracff(EvtId parent, EvtId daught, double q2, double mass, double *f1, double *f2, double *f3, double *g1, double *g2, double *g3)=0
double normalizedProb(const EvtSpinDensity &d)
|