diff --git a/EvtGenModels/EvtSSD_DirectCP.hh b/EvtGenModels/EvtSSD_DirectCP.hh --- a/EvtGenModels/EvtSSD_DirectCP.hh +++ b/EvtGenModels/EvtSSD_DirectCP.hh @@ -38,9 +38,12 @@ void decay( EvtParticle* p ) override; std::string getParamName( int i ) override; + protected: + void calcAmp( const EvtParticle& parent, EvtAmp& amp ) const; + private: - bool isB0Mixed( EvtParticle* p ); - bool isBsMixed( EvtParticle* p ); + bool isB0Mixed( const EvtParticle& p ); + bool isBsMixed( const EvtParticle& p ); //Arguments diff --git a/History.md b/History.md --- a/History.md +++ b/History.md @@ -11,6 +11,11 @@ === ## R02-02-0X +1 March 2023 Fernando Abudinen +* D92: Bugfix probmax issue for TENSOR daughter in EvtSSD_DirectCP model. + Calculation of amplitude in EvtSSD_DirectCP model moved to dedicated calcAmp function. + Got rid of a few static variables along the way. + 3 Feb 2023 John Back * D91: Check for non-zero momentum for EvtSLBaryonAmp parent spin density matrix. Print out integrals of JSON test histograms. diff --git a/src/EvtGenModels/EvtSSD_DirectCP.cpp b/src/EvtGenModels/EvtSSD_DirectCP.cpp --- a/src/EvtGenModels/EvtSSD_DirectCP.cpp +++ b/src/EvtGenModels/EvtSSD_DirectCP.cpp @@ -27,6 +27,7 @@ #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtRandom.hh" #include "EvtGenBase/EvtReport.hh" +#include "EvtGenBase/EvtScalarParticle.hh" #include "EvtGenBase/EvtTensor4C.hh" #include "EvtGenBase/EvtVector4C.hh" @@ -50,8 +51,8 @@ checkNArg( 1 ); checkNDaug( 2 ); - EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) ); - EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) ); + const EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) ); + const EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) ); if ( ( !( d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR ) ) || ( !( ( d2type == EvtSpinType::SCALAR ) || @@ -78,15 +79,89 @@ { double theProbMax = 1.; - EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) ); - EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) ); - if ( d1type == EvtSpinType::TENSOR || d2type == EvtSpinType::TENSOR ) - theProbMax *= 10; + const EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) ); + const EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) ); + + if ( d1type != EvtSpinType::SCALAR || d2type != EvtSpinType::SCALAR ) { + // Create a scalar parent at rest and initialize it + // Use noLifeTime() cludge to avoid generating random numbers + + EvtScalarParticle parent{}; + parent.noLifeTime(); + parent.init( getParentId(), + EvtVector4R( EvtPDL::getMass( getParentId() ), 0, 0, 0 ) ); + parent.setDiagonalSpinDensity(); + + // Create daughters and initialize amplitude + EvtAmp amp; + EvtId daughters[2] = { getDaug( 0 ), getDaug( 1 ) }; + amp.init( getParentId(), 2, daughters ); + parent.makeDaughters( 2, daughters ); + + const int scalarDaughterIndex = d1type == EvtSpinType::SCALAR ? 0 : 1; + const int nonScalarDaughterIndex = d1type == EvtSpinType::SCALAR ? 1 : 0; + + EvtParticle* scalarDaughter = parent.getDaug( scalarDaughterIndex ); + EvtParticle* nonScalarDaughter = parent.getDaug( nonScalarDaughterIndex ); + + scalarDaughter->noLifeTime(); + nonScalarDaughter->noLifeTime(); + + EvtSpinDensity rho; + rho.setDiag( parent.getSpinStates() ); + + // Momentum of daughters in parent's frame + const double m_parent = EvtPDL::getMass( getParentId() ); + const double m_sd = EvtPDL::getMass( getDaug( scalarDaughterIndex ) ); + const double m_nd = EvtPDL::getMass( getDaug( nonScalarDaughterIndex ) ); + const double pstar = + sqrt( pow( m_parent, 2 ) - pow( ( m_sd + m_nd ), 2 ) ) * + sqrt( pow( m_parent, 2 ) - pow( ( m_nd - m_sd ), 2 ) ) / + ( 2 * m_parent ); + + EvtVector4R p4_sd, p4_nd; + + const int nsteps = 16; + + double prob_max = 0; + double theta_max = 0; + + for ( int i = 0; i <= nsteps; i++ ) { + const double theta = i * EvtConst::pi / nsteps; + + p4_sd.set( sqrt( pow( pstar, 2 ) + pow( m_sd, 2 ) ), 0, + +pstar * sin( theta ), +pstar * cos( theta ) ); + + p4_nd.set( sqrt( pow( pstar, 2 ) + pow( m_nd, 2 ) ), 0, + -pstar * sin( theta ), -pstar * cos( theta ) ); + + scalarDaughter->init( getDaug( scalarDaughterIndex ), p4_sd ); + nonScalarDaughter->init( getDaug( nonScalarDaughterIndex ), p4_nd ); + + calcAmp( parent, amp ); + + const double i_prob = rho.normalizedProb( amp.getSpinDensity() ); + + if ( i_prob > prob_max ) { + prob_max = i_prob; + theta_max = theta; + } + } + + EvtGenReport( EVTGEN_INFO, "EvtGen" ) + << " - probability " << prob_max + << " found at theta = " << theta_max << std::endl; + theProbMax *= 1.01 * prob_max; + } setProbMax( theProbMax ); + + EvtGenReport( EVTGEN_INFO, "EvtGen" ) + << " EvtSSD_DirectCP - set up maximum probability to " << theProbMax + << std::endl; } -void EvtSSD_DirectCP::decay( EvtParticle* p ) +void EvtSSD_DirectCP::decay( EvtParticle* parent ) { bool flip = false; EvtId daugs[2]; @@ -103,97 +178,54 @@ } if ( flip ) { - if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) { - p->getParent()->setId( EvtPDL::chargeConj( p->getParent()->getId() ) ); - p->setId( EvtPDL::chargeConj( p->getId() ) ); + if ( ( isB0Mixed( *parent ) ) || ( isBsMixed( *parent ) ) ) { + parent->getParent()->setId( + EvtPDL::chargeConj( parent->getParent()->getId() ) ); + parent->setId( EvtPDL::chargeConj( parent->getId() ) ); } else { - p->setId( EvtPDL::chargeConj( p->getId() ) ); + parent->setId( EvtPDL::chargeConj( parent->getId() ) ); } - } - - if ( !flip ) { - daugs[0] = getDaug( 0 ); - daugs[1] = getDaug( 1 ); - } else { daugs[0] = EvtPDL::chargeConj( getDaug( 0 ) ); daugs[1] = EvtPDL::chargeConj( getDaug( 1 ) ); - } - - EvtParticle* d; - p->initializePhaseSpace( 2, daugs ); - - EvtVector4R p4_parent = p->getP4Restframe(); - double m_parent = p4_parent.mass(); - - EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) ); - - EvtVector4R momv; - EvtVector4R moms; - - if ( d2type == EvtSpinType::SCALAR ) { - d2type = EvtPDL::getSpinType( getDaug( 0 ) ); - d = p->getDaug( 0 ); - momv = d->getP4(); - moms = p->getDaug( 1 )->getP4(); } else { - d = p->getDaug( 1 ); - momv = d->getP4(); - moms = p->getDaug( 0 )->getP4(); - } - - if ( d2type == EvtSpinType::SCALAR ) { - vertex( 1. ); - } - - if ( d2type == EvtSpinType::VECTOR ) { - double norm = momv.mass() / ( momv.d3mag() * p->mass() ); - - vertex( 0, norm * p4_parent * ( d->epsParent( 0 ) ) ); - vertex( 1, norm * p4_parent * ( d->epsParent( 1 ) ) ); - vertex( 2, norm * p4_parent * ( d->epsParent( 2 ) ) ); + daugs[0] = getDaug( 0 ); + daugs[1] = getDaug( 1 ); } - if ( d2type == EvtSpinType::TENSOR ) { - double norm = d->mass() * d->mass() / - ( m_parent * d->getP4().d3mag() * d->getP4().d3mag() ); + parent->initializePhaseSpace( 2, daugs ); - vertex( 0, norm * d->epsTensorParent( 0 ).cont1( p4_parent ) * p4_parent ); - vertex( 1, norm * d->epsTensorParent( 1 ).cont1( p4_parent ) * p4_parent ); - vertex( 2, norm * d->epsTensorParent( 2 ).cont1( p4_parent ) * p4_parent ); - vertex( 3, norm * d->epsTensorParent( 3 ).cont1( p4_parent ) * p4_parent ); - vertex( 4, norm * d->epsTensorParent( 4 ).cont1( p4_parent ) * p4_parent ); - } + calcAmp( *parent, _amp2 ); } -bool EvtSSD_DirectCP::isB0Mixed( EvtParticle* p ) +bool EvtSSD_DirectCP::isB0Mixed( const EvtParticle& p ) { - if ( !( p->getParent() ) ) + if ( !( p.getParent() ) ) return false; - static EvtId B0 = EvtPDL::getId( "B0" ); - static EvtId B0B = EvtPDL::getId( "anti-B0" ); + static const EvtId B0 = EvtPDL::getId( "B0" ); + static const EvtId B0B = EvtPDL::getId( "anti-B0" ); - if ( ( p->getId() != B0 ) && ( p->getId() != B0B ) ) + if ( ( p.getId() != B0 ) && ( p.getId() != B0B ) ) return false; - if ( ( p->getParent()->getId() == B0 ) || ( p->getParent()->getId() == B0B ) ) + if ( ( p.getParent()->getId() == B0 ) || ( p.getParent()->getId() == B0B ) ) return true; return false; } -bool EvtSSD_DirectCP::isBsMixed( EvtParticle* p ) +bool EvtSSD_DirectCP::isBsMixed( const EvtParticle& p ) { - if ( !( p->getParent() ) ) + if ( !( p.getParent() ) ) return false; - static EvtId BS0 = EvtPDL::getId( "B_s0" ); - static EvtId BSB = EvtPDL::getId( "anti-B_s0" ); + static const EvtId BS0 = EvtPDL::getId( "B_s0" ); + static const EvtId BSB = EvtPDL::getId( "anti-B_s0" ); - if ( ( p->getId() != BS0 ) && ( p->getId() != BSB ) ) + if ( ( p.getId() != BS0 ) && ( p.getId() != BSB ) ) return false; - if ( ( p->getParent()->getId() == BS0 ) || ( p->getParent()->getId() == BSB ) ) + if ( ( p.getParent()->getId() == BS0 ) || ( p.getParent()->getId() == BSB ) ) return true; return false; @@ -208,3 +240,53 @@ return ""; } } + +void EvtSSD_DirectCP::calcAmp( const EvtParticle& parent, EvtAmp& amp ) const +{ + const EvtSpinType::spintype d1type = EvtPDL::getSpinType( + parent.getDaug( 0 )->getId() ); + const EvtSpinType::spintype d2type = EvtPDL::getSpinType( + parent.getDaug( 1 )->getId() ); + + if ( d1type == EvtSpinType::SCALAR && d2type == EvtSpinType::SCALAR ) { + amp.vertex( 1. ); + return; + } + + const EvtVector4R p4_parent = parent.getP4Restframe(); + + const int nonScalarDaughterIndex = d1type == EvtSpinType::SCALAR ? 1 : 0; + + const EvtParticle& daughter = *parent.getDaug( nonScalarDaughterIndex ); + + const EvtSpinType::spintype nonScalarType = EvtPDL::getSpinType( + daughter.getId() ); + + if ( nonScalarType == EvtSpinType::VECTOR ) { + const EvtVector4R momv = daughter.getP4(); + const double norm = momv.mass() / ( momv.d3mag() * parent.mass() ); + + amp.vertex( 0, norm * p4_parent * ( daughter.epsParent( 0 ) ) ); + amp.vertex( 1, norm * p4_parent * ( daughter.epsParent( 1 ) ) ); + amp.vertex( 2, norm * p4_parent * ( daughter.epsParent( 2 ) ) ); + + } + + // This is for the EvtSpinType::TENSOR case. + else { + const double norm = daughter.mass() * daughter.mass() / + ( p4_parent.mass() * daughter.getP4().d3mag() * + daughter.getP4().d3mag() ); + + amp.vertex( 0, norm * daughter.epsTensorParent( 0 ).cont1( p4_parent ) * + p4_parent ); + amp.vertex( 1, norm * daughter.epsTensorParent( 1 ).cont1( p4_parent ) * + p4_parent ); + amp.vertex( 2, norm * daughter.epsTensorParent( 2 ).cont1( p4_parent ) * + p4_parent ); + amp.vertex( 3, norm * daughter.epsTensorParent( 3 ).cont1( p4_parent ) * + p4_parent ); + amp.vertex( 4, norm * daughter.epsTensorParent( 4 ).cont1( p4_parent ) * + p4_parent ); + } +}