Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenBase/EvtDecayAmp.cpp
Show All 34 Lines | void EvtDecayAmp::makeDecay( EvtParticle* p, bool recursive ) | ||||
//original default value | //original default value | ||||
int ntimes = 10000; | int ntimes = 10000; | ||||
int more; | int more; | ||||
EvtSpinDensity rho; | EvtSpinDensity rho; | ||||
double prob, prob_max; | double prob, prob_max; | ||||
_amp2.init( p->getId(), getNDaug(), getDaugs() ); | m_amp2.init( p->getId(), getNDaug(), getDaugs() ); | ||||
do { | do { | ||||
_daugsDecayedByParentModel = false; | m_daugsDecayedByParentModel = false; | ||||
_weight = 1.0; | m_weight = 1.0; | ||||
decay( p ); | decay( p ); | ||||
rho = _amp2.getSpinDensity(); | rho = m_amp2.getSpinDensity(); | ||||
prob = p->getSpinDensityForward().normalizedProb( rho ); | prob = p->getSpinDensityForward().normalizedProb( rho ); | ||||
if ( prob < 0.0 ) { | if ( prob < 0.0 ) { | ||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | ||||
<< "Negative prob:" << p->getId().getId() << " " | << "Negative prob:" << p->getId().getId() << " " | ||||
<< p->getChannel() << endl; | << p->getChannel() << endl; | ||||
Show All 23 Lines | do { | ||||
if ( p->getParent() != 0 ) { | if ( p->getParent() != 0 ) { | ||||
EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) | EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) | ||||
<< "parent:" | << "parent:" | ||||
<< EvtPDL::name( p->getParent()->getId() ).c_str() << endl; | << EvtPDL::name( p->getParent()->getId() ).c_str() << endl; | ||||
EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) | EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) | ||||
<< "parent channel :" << p->getParent()->getChannel() | << "parent channel :" << p->getParent()->getChannel() | ||||
<< endl; | << endl; | ||||
size_t i; | |||||
EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) << "parent daughters :"; | EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) << "parent daughters :"; | ||||
for ( i = 0; i < p->getParent()->getNDaug(); i++ ) { | for ( size_t i = 0; i < p->getParent()->getNDaug(); i++ ) { | ||||
EvtGenReport( EVTGEN_DEBUG, "" ) | EvtGenReport( EVTGEN_DEBUG, "" ) | ||||
<< EvtPDL::name( p->getParent()->getDaug( i )->getId() ) | << EvtPDL::name( p->getParent()->getDaug( i )->getId() ) | ||||
.c_str() | .c_str() | ||||
<< " "; | << " "; | ||||
} | } | ||||
EvtGenReport( EVTGEN_DEBUG, "" ) << endl; | EvtGenReport( EVTGEN_DEBUG, "" ) << endl; | ||||
EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) << "daughters :"; | EvtGenReport( EVTGEN_DEBUG, "EvtGen" ) << "daughters :"; | ||||
Show All 11 Lines | do { | ||||
EvtGenReport( EVTGEN_DEBUG, "" ) | EvtGenReport( EVTGEN_DEBUG, "" ) | ||||
<< p->getDaug( i )->getP4() << " " | << p->getDaug( i )->getP4() << " " | ||||
<< p->getDaug( i )->mass(); | << p->getDaug( i )->mass(); | ||||
EvtGenReport( EVTGEN_DEBUG, "" ) << endl; | EvtGenReport( EVTGEN_DEBUG, "" ) << endl; | ||||
} | } | ||||
} | } | ||||
} | } | ||||
prob /= _weight; | prob /= m_weight; | ||||
prob_max = getProbMax( prob ); | prob_max = getProbMax( prob ); | ||||
p->setDecayProb( prob / prob_max ); | p->setDecayProb( prob / prob_max ); | ||||
more = prob < EvtRandom::Flat( prob_max ); | more = prob < EvtRandom::Flat( prob_max ); | ||||
ntimes--; | ntimes--; | ||||
Show All 22 Lines | void EvtDecayAmp::makeDecay( EvtParticle* p, bool recursive ) | ||||
} | } | ||||
EvtSpinDensity rho_list[10]; | EvtSpinDensity rho_list[10]; | ||||
rho_list[0] = p->getSpinDensityForward(); | rho_list[0] = p->getSpinDensityForward(); | ||||
EvtAmp ampcont; | EvtAmp ampcont; | ||||
if ( _amp2._pstates != 1 ) { | if ( m_amp2.m_pstates != 1 ) { | ||||
ampcont = _amp2.contract( 0, p->getSpinDensityForward() ); | ampcont = m_amp2.contract( 0, p->getSpinDensityForward() ); | ||||
} else { | } else { | ||||
ampcont = _amp2; | ampcont = m_amp2; | ||||
} | } | ||||
// it may be that the parent decay model has already | // it may be that the parent decay model has already | ||||
// done the decay - this should be rare and the | // done the decay - this should be rare and the | ||||
// model better know what it is doing.. | // model better know what it is doing.. | ||||
if ( !daugsDecayedByParentModel() ) { | if ( !daugsDecayedByParentModel() ) { | ||||
if ( recursive ) { | if ( recursive ) { | ||||
for ( size_t i = 0; i < p->getNDaug(); i++ ) { | for ( size_t i = 0; i < p->getNDaug(); i++ ) { | ||||
rho.setDim( _amp2.dstates[i] ); | rho.setDim( m_amp2.dstates[i] ); | ||||
if ( _amp2.dstates[i] == 1 ) { | if ( m_amp2.dstates[i] == 1 ) { | ||||
rho.set( 0, 0, EvtComplex( 1.0, 0.0 ) ); | rho.set( 0, 0, EvtComplex( 1.0, 0.0 ) ); | ||||
} else { | } else { | ||||
rho = ampcont.contract( _amp2._dnontrivial[i], _amp2 ); | rho = ampcont.contract( m_amp2.m_dnontrivial[i], m_amp2 ); | ||||
} | } | ||||
if ( !rho.check() ) { | if ( !rho.check() ) { | ||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | ||||
<< "-------start error-------" << endl; | << "-------start error-------" << endl; | ||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | ||||
<< "forward rho failed Check:" | << "forward rho failed Check:" | ||||
<< EvtPDL::name( p->getId() ).c_str() << " " | << EvtPDL::name( p->getId() ).c_str() << " " | ||||
Show All 22 Lines | if ( !daugsDecayedByParentModel() ) { | ||||
<< EvtPDL::name( grandParent->getId() ).c_str() | << EvtPDL::name( grandParent->getId() ).c_str() | ||||
<< endl; | << endl; | ||||
} | } | ||||
} | } | ||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | ||||
<< " EvtSpinDensity rho: " << rho; | << " EvtSpinDensity rho: " << rho; | ||||
_amp2.dump(); | m_amp2.dump(); | ||||
for ( size_t ii = 0; ii < i + 1; ii++ ) { | for ( size_t ii = 0; ii < i + 1; ii++ ) { | ||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | ||||
<< "rho_list[" << ii << "] = " << rho_list[ii]; | << "rho_list[" << ii << "] = " << rho_list[ii]; | ||||
} | } | ||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | ||||
<< "-------Done with error-------" << endl; | << "-------Done with error-------" << endl; | ||||
} | } | ||||
p->getDaug( i )->setSpinDensityForward( rho ); | p->getDaug( i )->setSpinDensityForward( rho ); | ||||
p->getDaug( i )->decay(); | p->getDaug( i )->decay(); | ||||
rho_list[i + 1] = p->getDaug( i )->getSpinDensityBackward(); | rho_list[i + 1] = p->getDaug( i )->getSpinDensityBackward(); | ||||
if ( _amp2.dstates[i] != 1 ) { | if ( m_amp2.dstates[i] != 1 ) { | ||||
ampcont = ampcont.contract( _amp2._dnontrivial[i], | ampcont = ampcont.contract( m_amp2.m_dnontrivial[i], | ||||
rho_list[i + 1] ); | rho_list[i + 1] ); | ||||
} | } | ||||
} | } | ||||
p->setSpinDensityBackward( _amp2.getBackwardSpinDensity( rho_list ) ); | p->setSpinDensityBackward( m_amp2.getBackwardSpinDensity( rho_list ) ); | ||||
if ( !p->getSpinDensityBackward().check() ) { | if ( !p->getSpinDensityBackward().check() ) { | ||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | ||||
<< "rho_backward failed Check" << p->getId().getId() << " " | << "rho_backward failed Check" << p->getId().getId() << " " | ||||
<< p->getChannel() << endl; | << p->getChannel() << endl; | ||||
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | EvtGenReport( EVTGEN_ERROR, "EvtGen" ) | ||||
<< p->getSpinDensityBackward(); | << p->getSpinDensityBackward(); | ||||
Show All 13 Lines |