Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenModels/EvtRareLbToLll.cpp
Show First 20 Lines • Show All 66 Lines • ▼ Show 20 Lines | void EvtRareLbToLll::init() | ||||
// 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 | // Work out whether we have electron mode | ||||
const EvtIdSet leptons{ "e-", "e+" }; | const EvtIdSet leptons{ "e-", "e+" }; | ||||
if ( leptons.contains( getDaug( 1 ) ) ) { | if ( leptons.contains( getDaug( 1 ) ) ) { | ||||
m_electronMode = true; | m_electronMode = true; | ||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | ||||
<< " EvtRareLbToLll has dielectron final state" << std::endl; | << " EvtRareLbToLll has dielectron final state" << std::endl; | ||||
} | } | ||||
std::string model{ "LQCD" }; | std::string model{ "LQCD" }; | ||||
if ( getNArg() == 1) { | if ( getNArg() == 1 ) { | ||||
model = getArgStr( 0 ); | 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 19 Lines | void EvtRareLbToLll::initProbMax() | ||||
if ( m_maxProbability == 0 ) { | if ( m_maxProbability == 0 ) { | ||||
EvtDiracParticle parent{}; | EvtDiracParticle parent{}; | ||||
parent.noLifeTime(); | parent.noLifeTime(); | ||||
parent.init( getParentId(), | parent.init( getParentId(), | ||||
EvtVector4R( EvtPDL::getMass( getParentId() ), 0, 0, 0 ) ); | EvtVector4R( EvtPDL::getMass( getParentId() ), 0, 0, 0 ) ); | ||||
parent.setDiagonalSpinDensity(); | parent.setDiagonalSpinDensity(); | ||||
EvtAmp amp; | EvtAmp amp; | ||||
EvtId daughters[3] = {getDaug( 0 ), getDaug( 1 ), getDaug( 2 )}; | EvtId daughters[3] = { getDaug( 0 ), getDaug( 1 ), getDaug( 2 ) }; | ||||
amp.init( getParentId(), 3, daughters ); | amp.init( getParentId(), 3, daughters ); | ||||
parent.makeDaughters( 3, daughters ); | parent.makeDaughters( 3, daughters ); | ||||
EvtParticle* lambda = parent.getDaug( 0 ); | EvtParticle* lambda = parent.getDaug( 0 ); | ||||
EvtParticle* lep1 = parent.getDaug( 1 ); | EvtParticle* lep1 = parent.getDaug( 1 ); | ||||
EvtParticle* lep2 = parent.getDaug( 2 ); | EvtParticle* lep2 = parent.getDaug( 2 ); | ||||
lambda->noLifeTime(); | lambda->noLifeTime(); | ||||
lep1->noLifeTime(); | lep1->noLifeTime(); | ||||
lep2->noLifeTime(); | lep2->noLifeTime(); | ||||
Show All 14 Lines | if ( m_maxProbability == 0 ) { | ||||
EvtGenReport( EVTGEN_INFO, "EvtGen" ) | EvtGenReport( EVTGEN_INFO, "EvtGen" ) | ||||
<< " EvtRareLbToLll is probing whole phase space ..." << std::endl; | << " EvtRareLbToLll is probing whole phase space ..." << std::endl; | ||||
double prob = 0; | double prob = 0; | ||||
const int nsteps = 5000; | const int nsteps = 5000; | ||||
for ( int i = 0; i <= nsteps; i++ ) { | for ( int i = 0; i <= nsteps; i++ ) { | ||||
const double q2 = q2min + i * ( q2max - q2min ) / nsteps; | const double q2 = q2min + i * ( q2max - q2min ) / nsteps; | ||||
const double elambda = ( M0 * M0 + mL * mL - q2 ) / 2 / M0; | const double elambda = ( M0 * M0 + mL * mL - q2 ) / 2 / M0; | ||||
double pstar{0}; | double pstar{ 0 }; | ||||
if ( i != 0 ) { | 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 != nsteps ) { | 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 ( int j = 0; j <= 45; j++ ) { | for ( int j = 0; j <= 45; j++ ) { | ||||
const double 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 ) ); | ||||
if ( i != nsteps) // 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() ); | ||||
// In case of electron mode add pole | // In case of electron mode add pole | ||||
if ( m_electronMode ) { | if ( m_electronMode ) { | ||||
prob /= 1.0 + m_poleSize / ( q2*q2 ); | 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 | ||||
<< " (" << nsteps * ( 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; | ||||
} | } | ||||
▲ Show 20 Lines • Show All 316 Lines • Show Last 20 Lines |