Page MenuHomeHEPForge
Differential D92 Diff 392 src/EvtGenModels/EvtSSD_DirectCP.cpp

Changeset View

Standalone View

View Options

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