Page MenuHomeHEPForge
Differential D58 Diff 259 src/EvtGenBase/EvtDecayBase.cpp

Changeset View

Standalone View

View Options

src/EvtGenBase/EvtDecayBase.cpp

Show All 38 Linesvoid EvtDecayBase::checkQ()
int i; int i;
int q = 0; int q = 0;
int qpar; int qpar;
//If there are no daughters (jetset etc) then we do not //If there are no daughters (jetset etc) then we do not
//want to do this test. Why? Because sometimes the parent //want to do this test. Why? Because sometimes the parent
//will have a nonzero charge. //will have a nonzero charge.
if ( _ndaug != 0 ) { if ( m_ndaug != 0 ) {
for ( i = 0; i < _ndaug; i++ ) { for ( i = 0; i < m_ndaug; i++ ) {
q += EvtPDL::chg3( _daug[i] ); q += EvtPDL::chg3( m_daug[i] );
} }
qpar = EvtPDL::chg3( _parent ); qpar = EvtPDL::chg3( m_parent );
if ( q != qpar ) { if ( q != qpar ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< _modelname.c_str() << " generator expected " << m_modelname.c_str() << " generator expected "
<< " charge to be conserved, found:" << endl; << " charge to be conserved, found:" << endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Parent charge of " << ( qpar / 3 ) << endl; << "Parent charge of " << ( qpar / 3 ) << endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Sum of daughter charge of " << ( q / 3 ) << endl; << "Sum of daughter charge of " << ( q / 3 ) << endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "The parent is " << EvtPDL::name( _parent ).c_str() << endl; << "The parent is " << EvtPDL::name( m_parent ).c_str() << endl;
for ( i = 0; i < _ndaug; i++ ) { for ( i = 0; i < m_ndaug; i++ ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Daughter " << EvtPDL::name( _daug[i] ).c_str() << endl; << "Daughter " << EvtPDL::name( m_daug[i] ).c_str() << endl;
} }
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Will terminate execution!" << endl; << "Will terminate execution!" << endl;
::abort(); ::abort();
} }
} }
} }
Show All 17 Lines if ( defaultprobmax && ntimes_prob <= 500 ) {
probmax *= 1.2; probmax *= 1.2;
} }
return 1000000.0 * prob; return 1000000.0 * prob;
} }
if ( prob > probmax * 1.0001 ) { if ( prob > probmax * 1.0001 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "prob > probmax:(" << prob << ">" << probmax << ")"; << "prob > probmax:(" << prob << ">" << probmax << ")";
EvtGenReport( EVTGEN_INFO, "" ) << "(" << _modelname.c_str() << ") "; EvtGenReport( EVTGEN_INFO, "" ) << "(" << m_modelname.c_str() << ") ";
EvtGenReport( EVTGEN_INFO, "" ) EvtGenReport( EVTGEN_INFO, "" )
<< EvtPDL::name( _parent ).c_str() << " -> "; << EvtPDL::name( m_parent ).c_str() << " -> ";
for ( i = 0; i < _ndaug; i++ ) { for ( i = 0; i < m_ndaug; i++ ) {
EvtGenReport( EVTGEN_INFO, "" ) EvtGenReport( EVTGEN_INFO, "" )
<< EvtPDL::name( _daug[i] ).c_str() << " "; << EvtPDL::name( m_daug[i] ).c_str() << " ";
} }
EvtGenReport( EVTGEN_INFO, "" ) << endl; EvtGenReport( EVTGEN_INFO, "" ) << endl;
if ( defaultprobmax ) if ( defaultprobmax )
probmax = prob; probmax = prob;
} }
ntimes_prob += 1; ntimes_prob += 1;
return probmax; return probmax;
} //getProbMax } //getProbMax
double EvtDecayBase::resetProbMax( double prob ) double EvtDecayBase::resetProbMax( double prob )
{ {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Reseting prob max\n"; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Reseting prob max\n";
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "prob > probmax:(" << prob << ">" << probmax << ")"; << "prob > probmax:(" << prob << ">" << probmax << ")";
EvtGenReport( EVTGEN_INFO, "" ) << "(" << _modelname.c_str() << ")"; EvtGenReport( EVTGEN_INFO, "" ) << "(" << m_modelname.c_str() << ")";
EvtGenReport( EVTGEN_INFO, "" ) << EvtPDL::getStdHep( _parent ) << "->"; EvtGenReport( EVTGEN_INFO, "" ) << EvtPDL::getStdHep( m_parent ) << "->";
for ( int i = 0; i < _ndaug; i++ ) { for ( int i = 0; i < m_ndaug; i++ ) {
EvtGenReport( EVTGEN_INFO, "" ) << EvtPDL::getStdHep( _daug[i] ) << " "; EvtGenReport( EVTGEN_INFO, "" ) << EvtPDL::getStdHep( m_daug[i] ) << " ";
} }
EvtGenReport( EVTGEN_INFO, "" ) << endl; EvtGenReport( EVTGEN_INFO, "" ) << endl;
probmax = 0.0; probmax = 0.0;
defaultprobmax = 0; defaultprobmax = 0;
ntimes_prob = 0; ntimes_prob = 0;
return prob; return prob;
▲ Show 20 Lines • Show All 71 Lines • ▼ Show 20 Lines
} //initProbMax } //initProbMax
void EvtDecayBase::saveDecayInfo( EvtId ipar, int ndaug, EvtId* daug, int narg, void EvtDecayBase::saveDecayInfo( EvtId ipar, int ndaug, EvtId* daug, int narg,
std::vector<std::string>& args, std::vector<std::string>& args,
std::string name, double brfr ) std::string name, double brfr )
{ {
int i; int i;
_brfr = brfr; m_brfr = brfr;
_ndaug = ndaug; m_ndaug = ndaug;
_narg = narg; m_narg = narg;
_parent = ipar; m_parent = ipar;
_dsum = 0; m_dsum = 0;
if ( _ndaug > 0 ) { if ( m_ndaug > 0 ) {
_daug.resize( _ndaug ); m_daug.resize( m_ndaug );
for ( i = 0; i < _ndaug; i++ ) { for ( i = 0; i < m_ndaug; i++ ) {
_daug[i] = daug[i]; m_daug[i] = daug[i];
_dsum += daug[i].getAlias(); m_dsum += daug[i].getAlias();
} }
} else { } else {
_daug.clear(); m_daug.clear();
} }
if ( _narg > 0 ) { if ( m_narg > 0 ) {
_args.resize( _narg + 1 ); m_args.resize( m_narg + 1 );
for ( i = 0; i < _narg; i++ ) { for ( i = 0; i < m_narg; i++ ) {
_args[i] = args[i]; m_args[i] = args[i];
} }
} else { } else {
_args.clear(); m_args.clear();
} }
_modelname = name; m_modelname = name;
this->init(); this->init();
this->initProbMax(); this->initProbMax();
if ( _chkCharge ) { if ( m_chkCharge ) {
this->checkQ(); this->checkQ();
} }
if ( defaultprobmax ) { if ( defaultprobmax ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "No default probmax for "; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "No default probmax for ";
EvtGenReport( EVTGEN_INFO, "" ) << "(" << _modelname.c_str() << ") "; EvtGenReport( EVTGEN_INFO, "" ) << "(" << m_modelname.c_str() << ") ";
EvtGenReport( EVTGEN_INFO, "" ) EvtGenReport( EVTGEN_INFO, "" )
<< EvtPDL::name( _parent ).c_str() << " -> "; << EvtPDL::name( m_parent ).c_str() << " -> ";
for ( i = 0; i < _ndaug; i++ ) { for ( i = 0; i < m_ndaug; i++ ) {
EvtGenReport( EVTGEN_INFO, "" ) EvtGenReport( EVTGEN_INFO, "" )
<< EvtPDL::name( _daug[i] ).c_str() << " "; << EvtPDL::name( m_daug[i] ).c_str() << " ";
} }
EvtGenReport( EVTGEN_INFO, "" ) << endl; EvtGenReport( EVTGEN_INFO, "" ) << endl;
EvtGenReport( EVTGEN_INFO, "" ) EvtGenReport( EVTGEN_INFO, "" )
<< "This is fine for development, but must be provided for production." << "This is fine for development, but must be provided for production."
<< endl; << endl;
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Never fear though - the decay will use the \n"; << "Never fear though - the decay will use the \n";
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "500 iterations to build up a good probmax \n"; << "500 iterations to build up a good probmax \n";
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "before accepting a decay. " << endl; << "before accepting a decay. " << endl;
} }
} }
EvtDecayBase::EvtDecayBase() EvtDecayBase::EvtDecayBase()
{ {
//the default is that the user module does _not_ set //the default is that the user module does _not_ set
// any probmax. // any probmax.
defaultprobmax = 1; defaultprobmax = 1;
ntimes_prob = 0; ntimes_prob = 0;
probmax = 0.0; probmax = 0.0;
_photos = 0; m_photos = 0;
_verbose = 0; m_verbose = 0;
_summary = 0; m_summary = 0;
_parent = EvtId( -1, -1 ); m_parent = EvtId( -1, -1 );
_ndaug = 0; m_ndaug = 0;
_narg = 0; m_narg = 0;
_modelname = "**********"; m_modelname = "**********";
//Default is to check that charge is conserved //Default is to check that charge is conserved
_chkCharge = 1; m_chkCharge = 1;
//statistics collection! //statistics collection!
max_prob = 0.0; max_prob = 0.0;
sum_prob = 0.0; sum_prob = 0.0;
} }
void EvtDecayBase::printSummary() const void EvtDecayBase::printSummary() const
{ {
if ( ntimes_prob > 0 ) { if ( ntimes_prob > 0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Calls = " << ntimes_prob << "Calls = " << ntimes_prob
<< " eff: " << sum_prob / ( probmax * ntimes_prob ) << " eff: " << sum_prob / ( probmax * ntimes_prob )
<< " frac. max:" << max_prob / probmax; << " frac. max:" << max_prob / probmax;
EvtGenReport( EVTGEN_INFO, "" ) EvtGenReport( EVTGEN_INFO, "" )
<< " probmax:" << probmax << " max:" << max_prob << " : "; << " probmax:" << probmax << " max:" << max_prob << " : ";
} }
printInfo(); printInfo();
} }
void EvtDecayBase::printInfo() const void EvtDecayBase::printInfo() const
{ {
EvtGenReport( EVTGEN_INFO, "" ) << EvtPDL::name( _parent ).c_str() << " -> "; EvtGenReport( EVTGEN_INFO, "" ) << EvtPDL::name( m_parent ).c_str() << " -> ";
for ( int i = 0; i < _ndaug; i++ ) { for ( int i = 0; i < m_ndaug; i++ ) {
EvtGenReport( EVTGEN_INFO, "" ) EvtGenReport( EVTGEN_INFO, "" )
<< EvtPDL::name( _daug[i] ).c_str() << " "; << EvtPDL::name( m_daug[i] ).c_str() << " ";
} }
EvtGenReport( EVTGEN_INFO, "" ) << " (" << _modelname.c_str() << ")" << endl; EvtGenReport( EVTGEN_INFO, "" ) << " (" << m_modelname.c_str() << ")" << endl;
} }
void EvtDecayBase::setProbMax( double prbmx ) void EvtDecayBase::setProbMax( double prbmx )
{ {
defaultprobmax = 0; defaultprobmax = 0;
probmax = prbmx; probmax = prbmx;
} }
▲ Show 20 Lines • Show All 135 Lines • ▼ Show 20 Lines else {
p->setFirstOrNot(); p->setFirstOrNot();
// if only one daughter do it // if only one daughter do it
if ( ndaugs == 1 ) { if ( ndaugs == 1 ) {
masses[0] = p->mass(); masses[0] = p->mass();
return; return;
} }
//until we get a combo whose masses are less than _parent mass. //until we get a combo whose masses are less than m_parent mass.
do { do {
mass_sum = 0.0; mass_sum = 0.0;
for ( i = 0; i < ndaugs; i++ ) { for ( i = 0; i < ndaugs; i++ ) {
masses[i] = EvtPDL::getMass( daugs[i] ); masses[i] = EvtPDL::getMass( daugs[i] );
mass_sum = mass_sum + masses[i]; mass_sum = mass_sum + masses[i];
} }
Show All 33 Lines else {
} while ( mass_sum > p->mass() ); } while ( mass_sum > p->mass() );
} //else } //else
return; return;
} }
void EvtDecayBase::checkNArg( int a1, int a2, int a3, int a4 ) void EvtDecayBase::checkNArg( int a1, int a2, int a3, int a4 )
{ {
if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) { if ( m_narg != a1 && m_narg != a2 && m_narg != a3 && m_narg != a4 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< _modelname.c_str() << " generator expected " << endl; << m_modelname.c_str() << " generator expected " << endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << a1 << endl; EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << a1 << endl;
; ;
if ( a2 > -1 ) { if ( a2 > -1 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << " or " << a2 << endl; EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << " or " << a2 << endl;
} }
if ( a3 > -1 ) { if ( a3 > -1 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << " or " << a3 << endl; EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << " or " << a3 << endl;
} }
if ( a4 > -1 ) { if ( a4 > -1 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << " or " << a4 << endl; EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << " or " << a4 << endl;
} }
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< " arguments but found:" << _narg << endl; << " arguments but found:" << m_narg << endl;
printSummary(); printSummary();
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Will terminate execution!" << endl; << "Will terminate execution!" << endl;
::abort(); ::abort();
} }
} }
void EvtDecayBase::checkNDaug( int d1, int d2 ) void EvtDecayBase::checkNDaug( int d1, int d2 )
{ {
if ( _ndaug != d1 && _ndaug != d2 ) { if ( m_ndaug != d1 && m_ndaug != d2 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< _modelname.c_str() << " generator expected "; << m_modelname.c_str() << " generator expected ";
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << d1; EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << d1;
if ( d2 > -1 ) { if ( d2 > -1 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << " or " << d2; EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << " or " << d2;
} }
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< " daughters but found:" << _ndaug << endl; << " daughters but found:" << m_ndaug << endl;
printSummary(); printSummary();
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Will terminate execution!" << endl; << "Will terminate execution!" << endl;
::abort(); ::abort();
} }
} }
void EvtDecayBase::checkSpinParent( EvtSpinType::spintype sp ) void EvtDecayBase::checkSpinParent( EvtSpinType::spintype sp )
{ {
EvtSpinType::spintype parenttype = EvtPDL::getSpinType( getParentId() ); EvtSpinType::spintype parenttype = EvtPDL::getSpinType( getParentId() );
if ( parenttype != sp ) { if ( parenttype != sp ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< _modelname.c_str() << " did not get the correct parent spin\n"; << m_modelname.c_str() << " did not get the correct parent spin\n";
printSummary(); printSummary();
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Will terminate execution!" << endl; << "Will terminate execution!" << endl;
::abort(); ::abort();
} }
} }
void EvtDecayBase::checkSpinDaughter( int d1, EvtSpinType::spintype sp ) void EvtDecayBase::checkSpinDaughter( int d1, EvtSpinType::spintype sp )
{ {
EvtSpinType::spintype parenttype = EvtPDL::getSpinType( getDaug( d1 ) ); EvtSpinType::spintype parenttype = EvtPDL::getSpinType( getDaug( d1 ) );
if ( parenttype != sp ) { if ( parenttype != sp ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< _modelname.c_str() << m_modelname.c_str()
<< " did not get the correct daughter spin d=" << d1 << endl; << " did not get the correct daughter spin d=" << d1 << endl;
printSummary(); printSummary();
EvtGenReport( EVTGEN_ERROR, "EvtGen" ) EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Will terminate execution!" << endl; << "Will terminate execution!" << endl;
::abort(); ::abort();
} }
} }
double* EvtDecayBase::getArgs() double* EvtDecayBase::getArgs()
{ {
if ( !_argsD.empty() ) if ( !m_argsD.empty() )
return _argsD.data(); return m_argsD.data();
//The user has asked for a list of doubles - the arguments //The user has asked for a list of doubles - the arguments
//better all be doubles... //better all be doubles...
if ( _narg == 0 ) if ( m_narg == 0 )
return _argsD.data(); return m_argsD.data();
_argsD.resize( _narg ); m_argsD.resize( m_narg );
for ( int i = 0; i < _narg; i++ ) { for ( int i = 0; i < m_narg; i++ ) {
char* tc; char* tc;
_argsD[i] = strtod( _args[i].c_str(), &tc ); m_argsD[i] = strtod( m_args[i].c_str(), &tc );
} }
return _argsD.data(); return m_argsD.data();
} }
double EvtDecayBase::getArg( unsigned int j ) double EvtDecayBase::getArg( unsigned int j )
{ {
// Verify string // Verify string
if ( getParentId().getId() == 25 ) { if ( getParentId().getId() == 25 ) {
int i = 0; int i = 0;
++i; ++i;
} }
const char* str = _args[j].c_str(); const char* str = m_args[j].c_str();
int i = 0; int i = 0;
while ( str[i] != 0 ) { while ( str[i] != 0 ) {
if ( isalpha( str[i] ) && str[i] != 'e' ) { if ( isalpha( str[i] ) && str[i] != 'e' ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "String " << str << " is not a number" << endl; << "String " << str << " is not a number" << endl;
assert( 0 ); assert( 0 );
} }
i++; i++;
} }
char** tc = 0; char** tc = 0;
double result = strtod( _args[j].c_str(), tc ); double result = strtod( m_args[j].c_str(), tc );
if ( _storedArgs.size() < j + 1 ) { // then store the argument's value if ( m_storedArgs.size() < j + 1 ) { // then store the argument's value
_storedArgs.push_back( result ); m_storedArgs.push_back( result );
} }
return result; return result;
} }
bool EvtDecayBase::matchingDecay( const EvtDecayBase& other ) const bool EvtDecayBase::matchingDecay( const EvtDecayBase& other ) const
{ {
if ( _ndaug != other._ndaug ) if ( m_ndaug != other.m_ndaug )
return false; return false;
if ( _parent != other._parent ) if ( m_parent != other.m_parent )
return false; return false;
std::vector<int> useDs; std::vector<int> useDs;
for ( int i = 0; i < _ndaug; i++ ) for ( int i = 0; i < m_ndaug; i++ )
useDs.push_back( 0 ); useDs.push_back( 0 );
for ( int i = 0; i < _ndaug; i++ ) { for ( int i = 0; i < m_ndaug; i++ ) {
bool foundIt = false; bool foundIt = false;
for ( int j = 0; j < _ndaug; j++ ) { for ( int j = 0; j < m_ndaug; j++ ) {
if ( useDs[j] == 1 ) if ( useDs[j] == 1 )
continue; continue;
if ( _daug[i] == other._daug[j] && if ( m_daug[i] == other.m_daug[j] &&
_daug[i].getAlias() == other._daug[j].getAlias() ) { m_daug[i].getAlias() == other.m_daug[j].getAlias() ) {
foundIt = true; foundIt = true;
useDs[j] = 1; useDs[j] = 1;
break; break;
} }
} }
if ( foundIt == false ) if ( foundIt == false )
return false; return false;
} }
for ( int i = 0; i < _ndaug; i++ ) for ( int i = 0; i < m_ndaug; i++ )
if ( useDs[i] == 0 ) if ( useDs[i] == 0 )
return false; return false;
return true; return true;
} }