Changeset View
Changeset View
Standalone View
Standalone View
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; | const double q2min = ( m1 + m2 ) * ( m1 + m2 ); | ||||
tlatham: Would be good to at least initialise these to 0 here, or better move their declarations to… | |||||
double q2min = ( m1 + m2 ) * ( m1 + m2 ); | const double q2max = ( M0 - mL ) * ( M0 - mL ); | ||||
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; | |||||
double prob = 0; | double prob = 0; | ||||
Done Inline ActionsThese 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. | |||||
for ( i = 0; i <= 100; i++ ) { | const int nsteps = 5000; | ||||
q2 = q2min + i * ( q2max - q2min ) / 100.; | for ( int i = 0; i <= nsteps; i++ ) { | ||||
elambda = ( M0 * M0 + mL * mL - q2 ) / 2 / M0; | const double q2 = q2min + i * ( q2max - q2min ) / nsteps; | ||||
if ( i == 0 ) | const double elambda = ( M0 * M0 + mL * mL - q2 ) / 2 / M0; | ||||
pstar = 0; | double pstar{0}; | ||||
else | if ( i != 0 ) { | ||||
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 ( int j = 0; j <= 45; j++ ) { | ||||
theta = j * EvtConst::pi / 45; | const double 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 ); | const 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; | const EvtParticle* lp = 0; | ||||
EvtParticle* lm = 0; | const 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; | const double mb = 5.209; | ||||
Done Inline ActionsShould 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 Lines | void 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 All 24 Lines | if ( EvtSpinType::DIRAC == spin ) { | ||||
EvtVector4C T[8]; | EvtVector4C T[8]; | ||||
EvtVector4C H1[2][4]; // vector current // swaped | EvtVector4C H1[2][4]; // vector current // swaped | ||||
EvtVector4C H2[2][4]; // axial-vector | EvtVector4C H2[2][4]; // axial-vector | ||||
// Build hadronic amplitude | // Build hadronic amplitude | ||||
for ( int i = 0; i < 2; ++i ) { | for ( int i = 0; i < 2; ++i ) { | ||||
for ( int j = 0; j < 4; ++j ) { | for ( int j = 0; j < 4; ++j ) { | ||||
HadronicAmpRS( 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] + | ||||
cs * AC[3] * T[6] + cp * BC[3] * T[7] ); | cs * AC[3] * T[6] + cp * BC[3] * T[7] ); | ||||
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] + | ||||
cs * DC[3] * T[6] + cp * EC[3] * T[7] ); | cs * DC[3] * T[6] + cp * EC[3] * T[7] ); | ||||
Show All 26 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 Lines | void 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 |
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.