Page MenuHomeHEPForge
Differential D84 Diff 362 src/EvtGenModels/EvtRareLbToLll.cpp

Changeset View

Standalone View

View Options

src/EvtGenModels/EvtRareLbToLll.cpp

Show First 20 Lines • Show All 63 Lines • ▼ Show 20 Lines if ( !( spin == EvtSpinType::DIRAC || spin == EvtSpinType::RARITASCHWINGER ) ) {
<< std::endl; << std::endl;
} }
// We expect that the second and third daughters // We expect that the second and third daughters
// are the ell+ and ell- // are the ell+ and ell-
checkSpinDaughter( 1, EvtSpinType::DIRAC ); checkSpinDaughter( 1, EvtSpinType::DIRAC );
checkSpinDaughter( 2, EvtSpinType::DIRAC ); checkSpinDaughter( 2, EvtSpinType::DIRAC );
// Work out whether we have electron mode
const EvtIdSet leptons{ "e-", "e+" };
if ( leptons.contains( getDaug( 1 ) ) ) {
m_electronMode = true;
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< " EvtRareLbToLll has dielectron final state" << std::endl;
}
std::string model = getArgStr( 0 ); std::string model = getArgStr( 0 );
if ( model == "Gutsche" ) { if ( model == "Gutsche" ) {
ffmodel_ = std::make_unique<EvtRareLbToLllFFGutsche>(); ffmodel_ = std::make_unique<EvtRareLbToLllFFGutsche>();
} else if ( model == "LQCD" ) { } else if ( model == "LQCD" ) {
ffmodel_ = std::make_unique<EvtRareLbToLllFFlQCD>(); ffmodel_ = std::make_unique<EvtRareLbToLllFFlQCD>();
} else if ( model == "MR" ) { } else if ( model == "MR" ) {
ffmodel_ = std::make_unique<EvtRareLbToLllFF>(); ffmodel_ = std::make_unique<EvtRareLbToLllFF>();
} else { } else {
Show All 32 Lines if ( m_maxProbability == 0 ) {
EvtParticle* lep2 = parent.getDaug( 2 ); EvtParticle* lep2 = parent.getDaug( 2 );
lambda->noLifeTime(); lambda->noLifeTime();
lep1->noLifeTime(); lep1->noLifeTime();
lep2->noLifeTime(); lep2->noLifeTime();
EvtSpinDensity rho; EvtSpinDensity rho;
rho.setDiag( parent.getSpinStates() ); rho.setDiag( parent.getSpinStates() );
double M0 = EvtPDL::getMass( getParentId() ); const double M0 = EvtPDL::getMass( getParentId() );
double mL = EvtPDL::getMass( getDaug( 0 ) ); const double mL = EvtPDL::getMass( getDaug( 0 ) );
double m1 = EvtPDL::getMass( getDaug( 1 ) ); const double m1 = EvtPDL::getMass( getDaug( 1 ) );
double m2 = EvtPDL::getMass( getDaug( 2 ) ); const double m2 = EvtPDL::getMass( getDaug( 2 ) );
double q2, pstar, elambda, theta; double q2, pstar, elambda, theta;
double q2min = ( m1 + m2 ) * ( m1 + m2 ); const double q2min = ( m1 + m2 ) * ( m1 + m2 );
tlathamUnsubmitted
Done
Inline Actions

Would be good to at least initialise these to 0 here, or better move their declarations to where they are first used so they are initialised there.

tlatham: Would be good to at least initialise these to 0 here, or better move their declarations to…
double q2max = ( M0 - mL ) * ( M0 - mL ); const double q2max = ( M0 - mL ) * ( M0 - mL );
EvtVector4R p4lambda, p4lep1, p4lep2, boost; EvtVector4R p4lambda, p4lep1, p4lep2, boost;
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< " EvtRareLbToLll is probing whole phase space ..." << std::endl; << " EvtRareLbToLll is probing whole phase space ..." << std::endl;
int i, j; int i, j;
tlathamUnsubmitted
Done
Inline Actions

These loop counters should be declared in the loop init-statement not here.

tlatham: These loop counters should be declared in the loop init-statement not here.
double prob = 0; double prob = 0;
for ( i = 0; i <= 100; i++ ) { const int nsteps = 5000;
q2 = q2min + i * ( q2max - q2min ) / 100.; for ( i = 0; i <= nsteps; i++ ) {
q2 = q2min + i * ( q2max - q2min ) / nsteps;
elambda = ( M0 * M0 + mL * mL - q2 ) / 2 / M0; elambda = ( M0 * M0 + mL * mL - q2 ) / 2 / M0;
if ( i == 0 ) if ( i == 0 )
pstar = 0; pstar = 0;
else else
pstar = sqrt( q2 - ( m1 + m2 ) * ( m1 + m2 ) ) * pstar = sqrt( q2 - ( m1 + m2 ) * ( m1 + m2 ) ) *
sqrt( q2 - ( m1 - m2 ) * ( m1 - m2 ) ) / 2 / sqrt( q2 ); sqrt( q2 - ( m1 - m2 ) * ( m1 - m2 ) ) / 2 / sqrt( q2 );
boost.set( M0 - elambda, 0, 0, +sqrt( elambda * elambda - mL * mL ) ); boost.set( M0 - elambda, 0, 0, +sqrt( elambda * elambda - mL * mL ) );
if ( i != 100 ) { if ( i != nsteps ) {
p4lambda.set( elambda, 0, 0, p4lambda.set( elambda, 0, 0,
-sqrt( elambda * elambda - mL * mL ) ); -sqrt( elambda * elambda - mL * mL ) );
} else { } else {
p4lambda.set( mL, 0, 0, 0 ); p4lambda.set( mL, 0, 0, 0 );
} }
for ( j = 0; j <= 45; j++ ) { for ( j = 0; j <= 45; j++ ) {
theta = j * EvtConst::pi / 45; theta = j * EvtConst::pi / 45;
p4lep1.set( sqrt( pstar * pstar + m1 * m1 ), 0, p4lep1.set( sqrt( pstar * pstar + m1 * m1 ), 0,
+pstar * sin( theta ), +pstar * cos( theta ) ); +pstar * sin( theta ), +pstar * cos( theta ) );
p4lep2.set( sqrt( pstar * pstar + m2 * m2 ), 0, p4lep2.set( sqrt( pstar * pstar + m2 * m2 ), 0,
-pstar * sin( theta ), -pstar * cos( theta ) ); -pstar * sin( theta ), -pstar * cos( theta ) );
//std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " pstar: " << pstar << std::endl;
if ( i != 100 ) // At maximal q2 we are already in correct frame as Lambda and W/Zvirtual are at rest if ( i != nsteps) // At maximal q2 we are already in correct frame as Lambda and W/Zvirtual are at rest
{ {
p4lep1 = boostTo( p4lep1, boost ); p4lep1 = boostTo( p4lep1, boost );
p4lep2 = boostTo( p4lep2, boost ); p4lep2 = boostTo( p4lep2, boost );
} }
lambda->init( getDaug( 0 ), p4lambda ); lambda->init( getDaug( 0 ), p4lambda );
lep1->init( getDaug( 1 ), p4lep1 ); lep1->init( getDaug( 1 ), p4lep1 );
lep2->init( getDaug( 2 ), p4lep2 ); lep2->init( getDaug( 2 ), p4lep2 );
calcAmp( amp, &parent ); calcAmp( amp, parent );
prob = rho.normalizedProb( amp.getSpinDensity() ); prob = rho.normalizedProb( amp.getSpinDensity() );
//std::cout << "q2: " << q2 << " \t theta: " << theta << " \t prob: " << prob << std::endl; // In case of electron mode add pole
//std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " q2-q2min: " << q2-(m1+m2)*(m1+m2) << std::endl; if ( m_electronMode ) {
prob /= 1.0 + m_poleSize / ( q2*q2 );
}
if ( prob > m_maxProbability ) { if ( prob > m_maxProbability ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< " - probability " << prob << " found at q2 = " << q2 << " - probability " << prob << " found at q2 = " << q2
<< " (" << 100 * ( q2 - q2min ) / ( q2max - q2min ) << " (" << nsteps * ( q2 - q2min ) / ( q2max - q2min )
<< " %) and theta = " << theta * 180 / EvtConst::pi << " %) and theta = " << theta * 180 / EvtConst::pi
<< std::endl; << std::endl;
m_maxProbability = prob; m_maxProbability = prob;
} }
} }
//::abort();
} }
//m_poleSize = 0.04*q2min; m_maxProbability *= 1.05;
m_maxProbability *= 1.2;
} }
setProbMax( m_maxProbability ); setProbMax( m_maxProbability );
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< " EvtRareLbToLll set up maximum probability to " << m_maxProbability << " EvtRareLbToLll set up maximum probability to " << m_maxProbability
<< std::endl; << std::endl;
} }
void EvtRareLbToLll::decay( EvtParticle* parent ) void EvtRareLbToLll::decay( EvtParticle* parent )
{ {
// Phase space initialization depends on what leptons are
if ( m_electronMode ) {
setWeight( parent->initializePhaseSpace( getNDaug(), getDaugs(), false,
m_poleSize, 1, 2 ) );
} else {
parent->initializePhaseSpace( getNDaug(), getDaugs() ); parent->initializePhaseSpace( getNDaug(), getDaugs() );
}
calcAmp( _amp2, parent ); calcAmp( _amp2, *parent );
} }
bool EvtRareLbToLll::isParticle( EvtParticle* parent ) const bool EvtRareLbToLll::isParticle( const EvtParticle& parent ) const
{ {
static EvtIdSet partlist( "Lambda_b0" ); const EvtIdSet partlist{ "Lambda_b0" };
return partlist.contains( parent->getId() ); return partlist.contains( parent.getId() );
} }
void EvtRareLbToLll::calcAmp( EvtAmp& amp, EvtParticle* parent ) void EvtRareLbToLll::calcAmp( EvtAmp& amp, const EvtParticle& parent )
{ {
//parent->setDiagonalSpinDensity(); //parent->setDiagonalSpinDensity();
EvtParticle* lambda = parent->getDaug( 0 ); EvtParticle* lambda = parent.getDaug( 0 );
static EvtIdSet leptons( "e-", "mu-", "tau-" ); const EvtIdSet leptons{ "e-", "mu-", "tau-" };
const bool isparticle = isParticle( parent ); const bool isparticle = isParticle( parent );
EvtParticle* lp = 0; EvtParticle* lp = 0;
EvtParticle* lm = 0; EvtParticle* lm = 0;
if ( leptons.contains( parent->getDaug( 1 )->getId() ) ) { if ( leptons.contains( parent.getDaug( 1 )->getId() ) ) {
lp = parent->getDaug( 1 ); lp = parent.getDaug( 1 );
lm = parent->getDaug( 2 ); lm = parent.getDaug( 2 );
} else { } else {
lp = parent->getDaug( 2 ); lp = parent.getDaug( 2 );
lm = parent->getDaug( 1 ); lm = parent.getDaug( 1 );
} }
EvtVector4R P; EvtVector4R P;
P.set( parent->mass(), 0.0, 0.0, 0.0 ); P.set( parent.mass(), 0.0, 0.0, 0.0 );
EvtVector4R q = lp->getP4() + lm->getP4(); EvtVector4R q = lp->getP4() + lm->getP4();
double qsq = q.mass2(); const double qsq = q.mass2();
// Leptonic currents // Leptonic currents
EvtVector4C LV[2][2]; // \bar{\ell} \gamma^{\mu} \ell EvtVector4C LV[2][2]; // \bar{\ell} \gamma^{\mu} \ell
EvtVector4C LA[2][2]; // \bar{\ell} \gamma^{\mu} \gamma^{5} \ell EvtVector4C LA[2][2]; // \bar{\ell} \gamma^{\mu} \gamma^{5} \ell
for ( int i = 0; i < 2; ++i ) { for ( int i = 0; i < 2; ++i ) {
for ( int j = 0; j < 2; ++j ) { for ( int j = 0; j < 2; ++j ) {
if ( isparticle ) { if ( isparticle ) {
LV[i][j] = EvtLeptonVCurrent( lp->spParent( i ), LV[i][j] = EvtLeptonVCurrent( lp->spParent( i ),
lm->spParent( j ) ); lm->spParent( j ) );
LA[i][j] = EvtLeptonACurrent( lp->spParent( i ), LA[i][j] = EvtLeptonACurrent( lp->spParent( i ),
lm->spParent( j ) ); lm->spParent( j ) );
} else { } else {
LV[i][j] = EvtLeptonVCurrent( lp->spParent( 1 - i ), LV[i][j] = EvtLeptonVCurrent( lp->spParent( 1 - i ),
lm->spParent( 1 - j ) ); lm->spParent( 1 - j ) );
LA[i][j] = EvtLeptonACurrent( lp->spParent( 1 - i ), LA[i][j] = EvtLeptonACurrent( lp->spParent( 1 - i ),
lm->spParent( 1 - j ) ); lm->spParent( 1 - j ) );
} }
} }
} }
EvtRareLbToLllFF::FormFactors FF; EvtRareLbToLllFF::FormFactors FF;
//F, G, FT and GT //F, G, FT and GT
ffmodel_->getFF( parent, lambda, FF ); ffmodel_->getFF( parent, *lambda, FF );
EvtComplex C7eff = wcmodel_->GetC7Eff( qsq ); EvtComplex C7eff = wcmodel_->GetC7Eff( qsq );
EvtComplex C9eff = wcmodel_->GetC9Eff( qsq ); EvtComplex C9eff = wcmodel_->GetC9Eff( qsq );
EvtComplex C10eff = wcmodel_->GetC10Eff( qsq ); EvtComplex C10eff = wcmodel_->GetC10Eff( qsq );
EvtComplex AC[4]; EvtComplex AC[4];
EvtComplex BC[4]; EvtComplex BC[4];
EvtComplex DC[4]; EvtComplex DC[4];
EvtComplex EC[4]; EvtComplex EC[4];
// check to see if particle is same or opposite parity to Lb // check to see if particle is same or opposite parity to Lb
const int parity = ffmodel_->isNatural( lambda ) ? 1 : -1; const int parity = ffmodel_->isNatural( *lambda ) ? 1 : -1;
// Lambda spin type // Lambda spin type
const EvtSpinType::spintype spin = EvtPDL::getSpinType( lambda->getId() ); const EvtSpinType::spintype spin = EvtPDL::getSpinType( lambda->getId() );
static const double mb = 5.209; static const double mb = 5.209;
jbackUnsubmitted
Done
Inline Actions

Should we remove the static here? Hopefully there are no other statics used in this model?

jback: Should we remove the static here? Hopefully there are no other statics used in this model?
// Eq. 48 + 49 // Eq. 48 + 49
for ( unsigned int i = 0; i < 4; ++i ) { for ( unsigned int i = 0; i < 4; ++i ) {
if ( parity > 0 ) { if ( parity > 0 ) {
AC[i] = -2. * mb * C7eff * FF.FT_[i] / qsq + C9eff * FF.F_[i]; AC[i] = -2. * mb * C7eff * FF.FT_[i] / qsq + C9eff * FF.F_[i];
BC[i] = -2. * mb * C7eff * FF.GT_[i] / qsq - C9eff * FF.G_[i]; BC[i] = -2. * mb * C7eff * FF.GT_[i] / qsq - C9eff * FF.G_[i];
DC[i] = C10eff * FF.F_[i]; DC[i] = C10eff * FF.F_[i];
EC[i] = -C10eff * FF.G_[i]; EC[i] = -C10eff * FF.G_[i];
Show All 14 Linesvoid EvtRareLbToLll::calcAmp( EvtAmp& amp, const EvtParticle& parent )
if ( EvtSpinType::DIRAC == spin ) { if ( EvtSpinType::DIRAC == spin ) {
EvtVector4C H1[2][2]; // vector current EvtVector4C H1[2][2]; // vector current
EvtVector4C H2[2][2]; // axial-vector EvtVector4C H2[2][2]; // axial-vector
EvtVector4C T[6]; EvtVector4C T[6];
// Hadronic currents // Hadronic currents
for ( int i = 0; i < 2; ++i ) { for ( int i = 0; i < 2; ++i ) {
for ( int j = 0; j < 2; ++j ) { for ( int j = 0; j < 2; ++j ) {
HadronicAmp( parent, lambda, T, i, j ); HadronicAmp( parent, *lambda, T, i, j );
H1[i][j] = ( cv * AC[0] * T[0] + ca * BC[0] * T[1] + H1[i][j] = ( cv * AC[0] * T[0] + ca * BC[0] * T[1] +
cs * AC[1] * T[2] + cp * BC[1] * T[3] + cs * AC[1] * T[2] + cp * BC[1] * T[3] +
cs * AC[2] * T[4] + cp * BC[2] * T[5] ); cs * AC[2] * T[4] + cp * BC[2] * T[5] );
H2[i][j] = ( cv * DC[0] * T[0] + ca * EC[0] * T[1] + H2[i][j] = ( cv * DC[0] * T[0] + ca * EC[0] * T[1] +
cs * DC[1] * T[2] + cp * EC[1] * T[3] + cs * DC[1] * T[2] + cp * EC[1] * T[3] +
cs * DC[2] * T[4] + cp * EC[2] * T[5] ); cs * DC[2] * T[4] + cp * EC[2] * T[5] );
▲ Show 20 Lines • Show All 66 Lines • ▼ Show 20 Lines if ( EvtSpinType::DIRAC == spin ) {
<< std::endl; << std::endl;
} }
return; return;
} }
// spin 1/2 daughters // spin 1/2 daughters
void EvtRareLbToLll::HadronicAmp( EvtParticle* parent, EvtParticle* lambda, void EvtRareLbToLll::HadronicAmp( const EvtParticle& parent,
EvtVector4C* T, const int i, const int j ) const EvtParticle& lambda, EvtVector4C* T,
const int i, const int j )
{ {
const EvtDiracSpinor Sfinal = lambda->spParent( j ); const EvtDiracSpinor Sfinal = lambda.spParent( j );
const EvtDiracSpinor Sinit = parent->sp( i ); const EvtDiracSpinor Sinit = parent.sp( i );
const EvtVector4R L = lambda->getP4(); const EvtVector4R L = lambda.getP4();
EvtVector4R P; EvtVector4R P;
P.set( parent->mass(), 0.0, 0.0, 0.0 ); P.set( parent.mass(), 0.0, 0.0, 0.0 );
const double Pm = parent->mass(); const double Pm = parent.mass();
const double Lm = lambda->mass(); const double Lm = lambda.mass();
// \bar{u} \gamma^{\mu} u // \bar{u} \gamma^{\mu} u
T[0] = EvtLeptonVCurrent( Sfinal, Sinit ); T[0] = EvtLeptonVCurrent( Sfinal, Sinit );
// \bar{u} \gamma^{\mu}\gamma^{5} u // \bar{u} \gamma^{\mu}\gamma^{5} u
T[1] = EvtLeptonACurrent( Sfinal, Sinit ); T[1] = EvtLeptonACurrent( Sfinal, Sinit );
// \bar{u} v^{\mu} u // \bar{u} v^{\mu} u
Show All 12 Linesvoid EvtRareLbToLll::HadronicAmp( const EvtParticle& parent,
// v = p_{\Lambda_b}/m_{\Lambda_b} // v = p_{\Lambda_b}/m_{\Lambda_b}
// v^{\prime} = p_{\Lambda}/m_{\Lambda} // v^{\prime} = p_{\Lambda}/m_{\Lambda}
return; return;
} }
// spin 3/2 daughters // spin 3/2 daughters
void EvtRareLbToLll::HadronicAmpRS( EvtParticle* parent, EvtParticle* lambda, void EvtRareLbToLll::HadronicAmpRS( const EvtParticle& parent,
EvtVector4C* T, const int i, const int j ) const EvtParticle& lambda, EvtVector4C* T,
const int i, const int j )
{ {
const EvtRaritaSchwinger Sfinal = lambda->spRSParent( j ); const EvtRaritaSchwinger Sfinal = lambda.spRSParent( j );
const EvtDiracSpinor Sinit = parent->sp( i ); const EvtDiracSpinor Sinit = parent.sp( i );
EvtVector4R P; EvtVector4R P;
P.set( parent->mass(), 0.0, 0.0, 0.0 ); P.set( parent.mass(), 0.0, 0.0, 0.0 );
const EvtVector4R L = lambda->getP4(); const EvtVector4R L = lambda.getP4();
EvtTensor4C ID; EvtTensor4C ID;
ID.setdiag( 1.0, 1.0, 1.0, 1.0 ); ID.setdiag( 1.0, 1.0, 1.0, 1.0 );
EvtDiracSpinor Sprime; EvtDiracSpinor Sprime;
for ( int i = 0; i < 4; i++ ) { for ( int i = 0; i < 4; i++ ) {
Sprime.set_spinor( i, Sfinal.getVector( i ) * P ); Sprime.set_spinor( i, Sfinal.getVector( i ) * P );
} }
const double Pmsq = P.mass2(); const double Pmsq = P.mass2();
const double Pm = parent->mass(); const double Pm = parent.mass();
const double PmLm = Pm * lambda->mass(); const double PmLm = Pm * lambda.mass();
EvtVector4C V1, V2; EvtVector4C V1, V2;
for ( int i = 0; i < 4; i++ ) { for ( int i = 0; i < 4; i++ ) {
V1.set( i, EvtLeptonSCurrent( Sfinal.getSpinor( i ), Sinit ) ); V1.set( i, EvtLeptonSCurrent( Sfinal.getSpinor( i ), Sinit ) );
V2.set( i, EvtLeptonPCurrent( Sfinal.getSpinor( i ), Sinit ) ); V2.set( i, EvtLeptonPCurrent( Sfinal.getSpinor( i ), Sinit ) );
} }
Show All 30 Lines