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 );
+ }
+}