Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenModels/EvtSSD_DirectCP.cpp
Show All 21 Lines | |||||
#include "EvtGenBase/EvtCPUtil.hh" | #include "EvtGenBase/EvtCPUtil.hh" | ||||
#include "EvtGenBase/EvtConst.hh" | #include "EvtGenBase/EvtConst.hh" | ||||
#include "EvtGenBase/EvtGenKine.hh" | #include "EvtGenBase/EvtGenKine.hh" | ||||
#include "EvtGenBase/EvtPDL.hh" | #include "EvtGenBase/EvtPDL.hh" | ||||
#include "EvtGenBase/EvtParticle.hh" | #include "EvtGenBase/EvtParticle.hh" | ||||
#include "EvtGenBase/EvtRandom.hh" | #include "EvtGenBase/EvtRandom.hh" | ||||
#include "EvtGenBase/EvtReport.hh" | #include "EvtGenBase/EvtReport.hh" | ||||
#include "EvtGenBase/EvtScalarParticle.hh" | |||||
#include "EvtGenBase/EvtTensor4C.hh" | #include "EvtGenBase/EvtTensor4C.hh" | ||||
#include "EvtGenBase/EvtVector4C.hh" | #include "EvtGenBase/EvtVector4C.hh" | ||||
#include <stdlib.h> | #include <stdlib.h> | ||||
#include <string> | #include <string> | ||||
std::string EvtSSD_DirectCP::getName() | std::string EvtSSD_DirectCP::getName() | ||||
{ | { | ||||
return "SSD_DirectCP"; | return "SSD_DirectCP"; | ||||
} | } | ||||
EvtDecayBase* EvtSSD_DirectCP::clone() | EvtDecayBase* EvtSSD_DirectCP::clone() | ||||
{ | { | ||||
return new EvtSSD_DirectCP; | return new EvtSSD_DirectCP; | ||||
} | } | ||||
void EvtSSD_DirectCP::init() | void EvtSSD_DirectCP::init() | ||||
{ | { | ||||
// check that there is 1 argument and 2-body decay | // check that there is 1 argument and 2-body decay | ||||
checkNArg( 1 ); | checkNArg( 1 ); | ||||
checkNDaug( 2 ); | checkNDaug( 2 ); | ||||
EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) ); | const EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) ); | ||||
EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) ); | const EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) ); | ||||
if ( ( !( d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR ) ) || | if ( ( !( d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR ) ) || | ||||
( !( ( d2type == EvtSpinType::SCALAR ) || | ( !( ( d2type == EvtSpinType::SCALAR ) || | ||||
( d2type == EvtSpinType::VECTOR ) || | ( d2type == EvtSpinType::VECTOR ) || | ||||
( d2type == EvtSpinType::TENSOR ) ) ) || | ( d2type == EvtSpinType::TENSOR ) ) ) || | ||||
( !( ( d1type == EvtSpinType::SCALAR ) || | ( !( ( d1type == EvtSpinType::SCALAR ) || | ||||
( d1type == EvtSpinType::VECTOR ) || | ( d1type == EvtSpinType::VECTOR ) || | ||||
( d1type == EvtSpinType::TENSOR ) ) ) ) { | ( d1type == EvtSpinType::TENSOR ) ) ) ) { | ||||
Show All 10 Lines | void EvtSSD_DirectCP::init() | ||||
_acp = getArg( 0 ); // A_CP defined as A_CP = (BR(fbar)-BR(f))/(BR(fbar)+BR(f)) | _acp = getArg( 0 ); // A_CP defined as A_CP = (BR(fbar)-BR(f))/(BR(fbar)+BR(f)) | ||||
} | } | ||||
void EvtSSD_DirectCP::initProbMax() | void EvtSSD_DirectCP::initProbMax() | ||||
{ | { | ||||
double theProbMax = 1.; | double theProbMax = 1.; | ||||
EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) ); | const EvtSpinType::spintype d2type = EvtPDL::getSpinType( getDaug( 1 ) ); | ||||
EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) ); | const EvtSpinType::spintype d1type = EvtPDL::getSpinType( getDaug( 0 ) ); | ||||
if ( d1type == EvtSpinType::TENSOR || d2type == EvtSpinType::TENSOR ) | |||||
theProbMax *= 10; | 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 ); | 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; | bool flip = false; | ||||
EvtId daugs[2]; | EvtId daugs[2]; | ||||
// decide it is B or Bbar: | // decide it is B or Bbar: | ||||
if ( EvtRandom::Flat( 0., 1. ) < ( ( 1. - _acp ) / 2. ) ) { | if ( EvtRandom::Flat( 0., 1. ) < ( ( 1. - _acp ) / 2. ) ) { | ||||
// it is a B | // it is a B | ||||
if ( EvtPDL::getStdHep( getParentId() ) < 0 ) | if ( EvtPDL::getStdHep( getParentId() ) < 0 ) | ||||
flip = true; | flip = true; | ||||
} else { | } else { | ||||
// it is a Bbar | // it is a Bbar | ||||
if ( EvtPDL::getStdHep( getParentId() ) > 0 ) | if ( EvtPDL::getStdHep( getParentId() ) > 0 ) | ||||
flip = true; | flip = true; | ||||
} | } | ||||
if ( flip ) { | if ( flip ) { | ||||
if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) { | if ( ( isB0Mixed( *parent ) ) || ( isBsMixed( *parent ) ) ) { | ||||
p->getParent()->setId( EvtPDL::chargeConj( p->getParent()->getId() ) ); | parent->getParent()->setId( | ||||
p->setId( EvtPDL::chargeConj( p->getId() ) ); | EvtPDL::chargeConj( parent->getParent()->getId() ) ); | ||||
parent->setId( EvtPDL::chargeConj( parent->getId() ) ); | |||||
} else { | } 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[0] = EvtPDL::chargeConj( getDaug( 0 ) ); | ||||
daugs[1] = EvtPDL::chargeConj( getDaug( 1 ) ); | 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 { | } else { | ||||
d = p->getDaug( 1 ); | daugs[0] = getDaug( 0 ); | ||||
momv = d->getP4(); | daugs[1] = getDaug( 1 ); | ||||
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 ) ) ); | |||||
} | } | ||||
if ( d2type == EvtSpinType::TENSOR ) { | parent->initializePhaseSpace( 2, daugs ); | ||||
double norm = d->mass() * d->mass() / | |||||
( m_parent * d->getP4().d3mag() * d->getP4().d3mag() ); | |||||
vertex( 0, norm * d->epsTensorParent( 0 ).cont1( p4_parent ) * p4_parent ); | calcAmp( *parent, _amp2 ); | ||||
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 ); | |||||
} | |||||
} | } | ||||
bool EvtSSD_DirectCP::isB0Mixed( EvtParticle* p ) | bool EvtSSD_DirectCP::isB0Mixed( const EvtParticle& p ) | ||||
{ | { | ||||
if ( !( p->getParent() ) ) | if ( !( p.getParent() ) ) | ||||
return false; | return false; | ||||
static EvtId B0 = EvtPDL::getId( "B0" ); | static const EvtId B0 = EvtPDL::getId( "B0" ); | ||||
static EvtId B0B = EvtPDL::getId( "anti-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; | return false; | ||||
if ( ( p->getParent()->getId() == B0 ) || ( p->getParent()->getId() == B0B ) ) | if ( ( p.getParent()->getId() == B0 ) || ( p.getParent()->getId() == B0B ) ) | ||||
return true; | return true; | ||||
return false; | return false; | ||||
} | } | ||||
bool EvtSSD_DirectCP::isBsMixed( EvtParticle* p ) | bool EvtSSD_DirectCP::isBsMixed( const EvtParticle& p ) | ||||
{ | { | ||||
if ( !( p->getParent() ) ) | if ( !( p.getParent() ) ) | ||||
return false; | return false; | ||||
static EvtId BS0 = EvtPDL::getId( "B_s0" ); | static const EvtId BS0 = EvtPDL::getId( "B_s0" ); | ||||
static EvtId BSB = EvtPDL::getId( "anti-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; | return false; | ||||
if ( ( p->getParent()->getId() == BS0 ) || ( p->getParent()->getId() == BSB ) ) | if ( ( p.getParent()->getId() == BS0 ) || ( p.getParent()->getId() == BSB ) ) | ||||
return true; | return true; | ||||
return false; | return false; | ||||
} | } | ||||
std::string EvtSSD_DirectCP::getParamName( int i ) | std::string EvtSSD_DirectCP::getParamName( int i ) | ||||
{ | { | ||||
switch ( i ) { | switch ( i ) { | ||||
case 0: | case 0: | ||||
return "ACP"; | return "ACP"; | ||||
default: | default: | ||||
return ""; | 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 ); | |||||
} | |||||
} |