Page MenuHomeHEPForge
Differential D58 Diff 240 src/EvtGenExternal/EvtPythiaEngine.cpp

Changeset View

Standalone View

View Options

src/EvtGenExternal/EvtPythiaEngine.cpp

Show First 20 Lines • Show All 47 Lines • ▼ Show 20 LinesEvtPythiaEngine::EvtPythiaEngine( std::string xmlDir, bool convertPhysCodes,
// only decay aliased particles, which are in general "signal" // only decay aliased particles, which are in general "signal"
// decays different from those in the decay.dec file. // decays different from those in the decay.dec file.
// Even though it is not possible to have two different particle // Even though it is not possible to have two different particle
// versions in one Pythia generator, we can use two generators to // versions in one Pythia generator, we can use two generators to
// get the required behaviour. // get the required behaviour.
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Creating generic Pythia generator" << endl; << "Creating generic Pythia generator" << endl;
_genericPythiaGen = std::make_unique<Pythia8::Pythia>( xmlDir ); m_genericPythiaGen = std::make_unique<Pythia8::Pythia>( xmlDir );
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Creating alias Pythia generator" << endl; << "Creating alias Pythia generator" << endl;
_aliasPythiaGen = std::make_unique<Pythia8::Pythia>( xmlDir, false ); m_aliasPythiaGen = std::make_unique<Pythia8::Pythia>( xmlDir, false );
_thePythiaGenerator = 0; m_thePythiaGenerator = 0;
_daugPDGVector.clear(); m_daugPDGVector.clear();
_daugP4Vector.clear(); m_daugP4Vector.clear();
_convertPhysCodes = convertPhysCodes; m_convertPhysCodes = convertPhysCodes;
// Specify if we are going to use the random number generator (engine) // Specify if we are going to use the random number generator (engine)
// from EvtGen for Pythia 8. // from EvtGen for Pythia 8.
_useEvtGenRandom = useEvtGenRandom; m_useEvtGenRandom = useEvtGenRandom;
_evtgenRandom = std::make_unique<EvtPythiaRandom>(); m_evtgenRandom = std::make_unique<EvtPythiaRandom>();
_initialised = false; m_initialised = false;
} }
EvtPythiaEngine::~EvtPythiaEngine() EvtPythiaEngine::~EvtPythiaEngine()
{ {
_thePythiaGenerator = nullptr; m_thePythiaGenerator = nullptr;
this->clearDaughterVectors(); this->clearDaughterVectors();
this->clearPythiaModeMap(); this->clearPythiaModeMap();
} }
void EvtPythiaEngine::clearDaughterVectors() void EvtPythiaEngine::clearDaughterVectors()
{ {
_daugPDGVector.clear(); m_daugPDGVector.clear();
_daugP4Vector.clear(); m_daugP4Vector.clear();
} }
void EvtPythiaEngine::clearPythiaModeMap() void EvtPythiaEngine::clearPythiaModeMap()
{ {
PythiaModeMap::iterator iter; PythiaModeMap::iterator iter;
for ( iter = _pythiaModeMap.begin(); iter != _pythiaModeMap.end(); ++iter ) { for ( iter = m_pythiaModeMap.begin(); iter != m_pythiaModeMap.end(); ++iter ) {
std::vector<int> modeVector = iter->second; std::vector<int> modeVector = iter->second;
modeVector.clear(); modeVector.clear();
} }
_pythiaModeMap.clear(); m_pythiaModeMap.clear();
} }
void EvtPythiaEngine::initialise() void EvtPythiaEngine::initialise()
{ {
if ( _initialised ) { if ( m_initialised ) {
return; return;
} }
this->clearPythiaModeMap(); this->clearPythiaModeMap();
this->updateParticleLists(); this->updateParticleLists();
// Hadron-level processes only (hadronized, string fragmentation and secondary decays). // Hadron-level processes only (hadronized, string fragmentation and secondary decays).
// We do not want to generate the full pp or e+e- event structure etc.. // We do not want to generate the full pp or e+e- event structure etc..
_genericPythiaGen->readString( "ProcessLevel:all = off" ); m_genericPythiaGen->readString( "ProcessLevel:all = off" );
_aliasPythiaGen->readString( "ProcessLevel:all = off" ); m_aliasPythiaGen->readString( "ProcessLevel:all = off" );
// Turn off Pythia warnings, e.g. changes to particle properties // Turn off Pythia warnings, e.g. changes to particle properties
_genericPythiaGen->readString( "Print:quiet = on" ); m_genericPythiaGen->readString( "Print:quiet = on" );
_aliasPythiaGen->readString( "Print:quiet = on" ); m_aliasPythiaGen->readString( "Print:quiet = on" );
// Apply any other physics (or special particle) requirements/cuts etc.. // Apply any other physics (or special particle) requirements/cuts etc..
this->updatePhysicsParameters(); this->updatePhysicsParameters();
// Set the random number generator // Set the random number generator
if ( _useEvtGenRandom == true ) { if ( m_useEvtGenRandom == true ) {
_genericPythiaGen->setRndmEnginePtr( _evtgenRandom.get() ); m_genericPythiaGen->setRndmEnginePtr( m_evtgenRandom.get() );
_aliasPythiaGen->setRndmEnginePtr( _evtgenRandom.get() ); m_aliasPythiaGen->setRndmEnginePtr( m_evtgenRandom.get() );
} }
_genericPythiaGen->init(); m_genericPythiaGen->init();
_aliasPythiaGen->init(); m_aliasPythiaGen->init();
_initialised = true; m_initialised = true;
} }
bool EvtPythiaEngine::doDecay( EvtParticle* theParticle ) bool EvtPythiaEngine::doDecay( EvtParticle* theParticle )
{ {
// Store the mother particle within a Pythia8 Event object. // Store the mother particle within a Pythia8 Event object.
// Then do the hadron level decays. // Then do the hadron level decays.
// The EvtParticle must be a colour singlet (meson/baryon/lepton), i.e. not a gluon or quark // The EvtParticle must be a colour singlet (meson/baryon/lepton), i.e. not a gluon or quark
// We delete any daughters the particle may have, since we are asking Pythia // We delete any daughters the particle may have, since we are asking Pythia
// to generate the decay anew. Also note that _any_ Pythia decay allowed for the particle // to generate the decay anew. Also note that _any_ Pythia decay allowed for the particle
// will be generated and not the specific Pythia decay mode that EvtGen has already // will be generated and not the specific Pythia decay mode that EvtGen has already
// specified. This is necessary since we only want to initialise the Pythia decay table // specified. This is necessary since we only want to initialise the Pythia decay table
// once; all Pythia branching fractions for a given mother particle are renormalised to sum to 1.0. // once; all Pythia branching fractions for a given mother particle are renormalised to sum to 1.0.
// In EvtGen decay.dec files, it may be the case that Pythia decays are only used // In EvtGen decay.dec files, it may be the case that Pythia decays are only used
// for some of the particle decays (i.e. Pythia BF sum < 1.0). As we loop over many events, // for some of the particle decays (i.e. Pythia BF sum < 1.0). As we loop over many events,
// the total frequency for each Pythia decay mode will normalise correctly to what // the total frequency for each Pythia decay mode will normalise correctly to what
// we wanted via the specifications made to the decay.dec file, even though event-by-event // we wanted via the specifications made to the decay.dec file, even though event-by-event
// the EvtGen decay channel and the Pythia decay channel may be different. // the EvtGen decay channel and the Pythia decay channel may be different.
if ( _initialised == false ) { if ( m_initialised == false ) {
this->initialise(); this->initialise();
} }
if ( theParticle == 0 ) { if ( theParticle == 0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Error in EvtPythiaEngine::doDecay. The mother particle is null. Not doing any Pythia decay." << "Error in EvtPythiaEngine::doDecay. The mother particle is null. Not doing any Pythia decay."
<< endl; << endl;
return false; return false;
} }
// Delete EvtParticle daughters if they already exist // Delete EvtParticle daughters if they already exist
if ( theParticle->getNDaug() != 0 ) { if ( theParticle->getNDaug() != 0 ) {
bool keepChannel( false ); bool keepChannel( false );
theParticle->deleteDaughters( keepChannel ); theParticle->deleteDaughters( keepChannel );
} }
EvtId particleId = theParticle->getId(); EvtId particleId = theParticle->getId();
int isAlias = particleId.isAlias(); int isAlias = particleId.isAlias();
// Choose the generator depending if we have an aliased (parent) particle or not // Choose the generator depending if we have an aliased (parent) particle or not
_thePythiaGenerator = ( isAlias == 1 ? _aliasPythiaGen.get() m_thePythiaGenerator = ( isAlias == 1 ? m_aliasPythiaGen.get()
: _genericPythiaGen.get() ); : m_genericPythiaGen.get() );
// Need to use the reference to the Pythia8::Event object, // Need to use the reference to the Pythia8::Event object,
// otherwise it will just return a new empty, default event object. // otherwise it will just return a new empty, default event object.
Pythia8::Event& theEvent = _thePythiaGenerator->event; Pythia8::Event& theEvent = m_thePythiaGenerator->event;
theEvent.reset(); theEvent.reset();
// Initialise the event to be the particle rest frame // Initialise the event to be the particle rest frame
int PDGCode = EvtPDL::getStdHep( particleId ); int PDGCode = EvtPDL::getStdHep( particleId );
int status( 1 ); int status( 1 );
int colour( 0 ), anticolour( 0 ); int colour( 0 ), anticolour( 0 );
double px( 0.0 ), py( 0.0 ), pz( 0.0 ); double px( 0.0 ), py( 0.0 ), pz( 0.0 );
double m0 = theParticle->mass(); double m0 = theParticle->mass();
double E = m0; double E = m0;
theEvent.append( PDGCode, status, colour, anticolour, px, py, pz, E, m0 ); theEvent.append( PDGCode, status, colour, anticolour, px, py, pz, E, m0 );
// Generate the Pythia event // Generate the Pythia event
int iTrial( 0 ); int iTrial( 0 );
bool generatedEvent( false ); bool generatedEvent( false );
for ( iTrial = 0; iTrial < 10; iTrial++ ) { for ( iTrial = 0; iTrial < 10; iTrial++ ) {
generatedEvent = _thePythiaGenerator->next(); generatedEvent = m_thePythiaGenerator->next();
if ( generatedEvent ) { if ( generatedEvent ) {
break; break;
} }
} }
bool success( false ); bool success( false );
if ( generatedEvent ) { if ( generatedEvent ) {
Show All 22 Lines if ( generatedEvent ) {
success = true; success = true;
} }
return success; return success;
} }
void EvtPythiaEngine::storeDaughterInfo( EvtParticle* theParticle, int startInt ) void EvtPythiaEngine::storeDaughterInfo( EvtParticle* theParticle, int startInt )
{ {
Pythia8::Event& theEvent = _thePythiaGenerator->event; Pythia8::Event& theEvent = m_thePythiaGenerator->event;
std::vector<int> daugList = theEvent.daughterList( startInt ); std::vector<int> daugList = theEvent.daughterList( startInt );
std::vector<int>::iterator daugIter; std::vector<int>::iterator daugIter;
for ( daugIter = daugList.begin(); daugIter != daugList.end(); ++daugIter ) { for ( daugIter = daugList.begin(); daugIter != daugList.end(); ++daugIter ) {
int daugInt = *daugIter; int daugInt = *daugIter;
// Ask if the daughter is a quark or gluon. If so, recursively search again. // Ask if the daughter is a quark or gluon. If so, recursively search again.
Show All 17 Lines for ( daugIter = daugList.begin(); daugIter != daugList.end(); ++daugIter ) {
double px = daugParticle.px(); double px = daugParticle.px();
double py = daugParticle.py(); double py = daugParticle.py();
double pz = daugParticle.pz(); double pz = daugParticle.pz();
double E = daugParticle.e(); double E = daugParticle.e();
EvtVector4R daughterP4( E, px, py, pz ); EvtVector4R daughterP4( E, px, py, pz );
// Now store the EvtId and 4-momentum in the internal vectors // Now store the EvtId and 4-momentum in the internal vectors
_daugPDGVector.push_back( daugPDGInt ); m_daugPDGVector.push_back( daugPDGInt );
_daugP4Vector.push_back( daughterP4 ); m_daugP4Vector.push_back( daughterP4 );
// Set the status flag for the Pythia particle to let us know // Set the status flag for the Pythia particle to let us know
// that we have already considered it to avoid double counting. // that we have already considered it to avoid double counting.
daugParticle.status( 1000 ); daugParticle.status( 1000 );
} // Status code != 1000 } // Status code != 1000
} }
} }
} }
void EvtPythiaEngine::createDaughterEvtParticles( EvtParticle* theParent ) void EvtPythiaEngine::createDaughterEvtParticles( EvtParticle* theParent )
{ {
if ( theParent == 0 ) { if ( theParent == 0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Error in EvtPythiaEngine::createDaughterEvtParticles. The parent is null" << "Error in EvtPythiaEngine::createDaughterEvtParticles. The parent is null"
<< endl; << endl;
return; return;
} }
// Get the list of Pythia decay modes defined for this particle id alias. // Get the list of Pythia decay modes defined for this particle id alias.
// It would be easier to just use the decay channel number that Pythia chose to use // It would be easier to just use the decay channel number that Pythia chose to use
// for the particle decay, but this is not accessible from the Pythia interface at present. // for the particle decay, but this is not accessible from the Pythia interface at present.
int nDaughters = _daugPDGVector.size(); int nDaughters = m_daugPDGVector.size();
std::vector<EvtId> daugAliasIdVect( 0 ); std::vector<EvtId> daugAliasIdVect( 0 );
EvtId particleId = theParent->getId(); EvtId particleId = theParent->getId();
// Check to see if we have an anti-particle. If we do, charge conjugate the particle id to get the // Check to see if we have an anti-particle. If we do, charge conjugate the particle id to get the
// Pythia "alias" we can compare with the defined (particle) Pythia modes. // Pythia "alias" we can compare with the defined (particle) Pythia modes.
int PDGId = EvtPDL::getStdHep( particleId ); int PDGId = EvtPDL::getStdHep( particleId );
int aliasInt = particleId.getAlias(); int aliasInt = particleId.getAlias();
int pythiaAliasInt( aliasInt ); int pythiaAliasInt( aliasInt );
if ( PDGId < 0 ) { if ( PDGId < 0 ) {
// We have an anti-particle. // We have an anti-particle.
EvtId conjPartId = EvtPDL::chargeConj( particleId ); EvtId conjPartId = EvtPDL::chargeConj( particleId );
pythiaAliasInt = conjPartId.getAlias(); pythiaAliasInt = conjPartId.getAlias();
} }
std::vector<int> pythiaModes = _pythiaModeMap[pythiaAliasInt]; std::vector<int> pythiaModes = m_pythiaModeMap[pythiaAliasInt];
// Loop over all available Pythia decay modes and find the channel that matches // Loop over all available Pythia decay modes and find the channel that matches
// the daughter ids. Set each daughter id to also use the alias integer. // the daughter ids. Set each daughter id to also use the alias integer.
// This will then convert the Pythia generated channel to the EvtGen alias defined one. // This will then convert the Pythia generated channel to the EvtGen alias defined one.
std::vector<int>::iterator modeIter; std::vector<int>::iterator modeIter;
bool gotMode( false ); bool gotMode( false );
Show All 14 Lines for ( modeIter = pythiaModes.begin(); modeIter != pythiaModes.end();
// We need to make sure that the number of daughters match // We need to make sure that the number of daughters match
if ( nDaughters == nModeDaug ) { if ( nDaughters == nModeDaug ) {
int iModeDaug( 0 ); int iModeDaug( 0 );
for ( iModeDaug = 0; iModeDaug < nModeDaug; iModeDaug++ ) { for ( iModeDaug = 0; iModeDaug < nModeDaug; iModeDaug++ ) {
EvtId daugId = decayModel->getDaug( iModeDaug ); EvtId daugId = decayModel->getDaug( iModeDaug );
int daugPDGId = EvtPDL::getStdHep( daugId ); int daugPDGId = EvtPDL::getStdHep( daugId );
// Pythia has used the right PDG codes for this decay mode, even for conjugate modes // Pythia has used the right PDG codes for this decay mode, even for conjugate modes
int pythiaPDGId = _daugPDGVector[iModeDaug]; int pythiaPDGId = m_daugPDGVector[iModeDaug];
if ( daugPDGId == pythiaPDGId ) { if ( daugPDGId == pythiaPDGId ) {
daugAliasIdVect.push_back( daugId ); daugAliasIdVect.push_back( daugId );
} }
} // Loop over EvtGen mode daughters } // Loop over EvtGen mode daughters
int daugAliasSize = daugAliasIdVect.size(); int daugAliasSize = daugAliasIdVect.size();
Show All 12 Linesvoid EvtPythiaEngine::createDaughterEvtParticles( EvtParticle* theParent )
} // Loop over available Pythia modes } // Loop over available Pythia modes
if ( gotMode == false ) { if ( gotMode == false ) {
// We did not find a match for the daughter aliases. Just use the normal PDG codes // We did not find a match for the daughter aliases. Just use the normal PDG codes
// from the Pythia decay result // from the Pythia decay result
int iPyDaug( 0 ); int iPyDaug( 0 );
for ( iPyDaug = 0; iPyDaug < nDaughters; iPyDaug++ ) { for ( iPyDaug = 0; iPyDaug < nDaughters; iPyDaug++ ) {
int daugPDGCode = _daugPDGVector[iPyDaug]; int daugPDGCode = m_daugPDGVector[iPyDaug];
EvtId daugPyId = EvtPDL::evtIdFromStdHep( daugPDGCode ); EvtId daugPyId = EvtPDL::evtIdFromStdHep( daugPDGCode );
daugAliasIdVect.push_back( daugPyId ); daugAliasIdVect.push_back( daugPyId );
} }
} }
// Make the EvtParticle daughters (with correct alias id's). Their 4-momenta are uninitialised. // Make the EvtParticle daughters (with correct alias id's). Their 4-momenta are uninitialised.
theParent->makeDaughters( nDaughters, daugAliasIdVect ); theParent->makeDaughters( nDaughters, daugAliasIdVect );
// Now set the 4-momenta of the daughters. // Now set the 4-momenta of the daughters.
int iDaug( 0 ); int iDaug( 0 );
// Can use an iterator here, but we already had to use the vector size... // Can use an iterator here, but we already had to use the vector size...
for ( iDaug = 0; iDaug < nDaughters; iDaug++ ) { for ( iDaug = 0; iDaug < nDaughters; iDaug++ ) {
EvtParticle* theDaughter = theParent->getDaug( iDaug ); EvtParticle* theDaughter = theParent->getDaug( iDaug );
// Set the correct 4-momentum for each daughter particle. // Set the correct 4-momentum for each daughter particle.
if ( theDaughter != 0 ) { if ( theDaughter != 0 ) {
EvtId theDaugId = daugAliasIdVect[iDaug]; EvtId theDaugId = daugAliasIdVect[iDaug];
const EvtVector4R theDaugP4 = _daugP4Vector[iDaug]; const EvtVector4R theDaugP4 = m_daugP4Vector[iDaug];
theDaughter->init( theDaugId, theDaugP4 ); theDaughter->init( theDaugId, theDaugP4 );
} }
} }
} }
void EvtPythiaEngine::updateParticleLists() void EvtPythiaEngine::updateParticleLists()
{ {
// Use the EvtGen decay table (decay/user.dec) to update Pythia particle // Use the EvtGen decay table (decay/user.dec) to update Pythia particle
// decay modes. Also, make sure the generic and alias Pythia generators // decay modes. Also, make sure the generic and alias Pythia generators
// use the same particle data entries as defined by EvtGen (evt.pdl). // use the same particle data entries as defined by EvtGen (evt.pdl).
// Loop over all entries in the EvtPDL particle data table. // Loop over all entries in the EvtPDL particle data table.
// Aliases are added at the end with id numbers equal to the // Aliases are added at the end with id numbers equal to the
// original particle, but with alias integer = lastPDLEntry+1 etc.. // original particle, but with alias integer = lastPDLEntry+1 etc..
int iPDL; int iPDL;
int nPDL = EvtPDL::entries(); int nPDL = EvtPDL::entries();
// Reset the _addedPDGCodes map that keeps track // Reset the m_addedPDGCodes map that keeps track
// of any new particles added to the Pythia input data stream // of any new particles added to the Pythia input data stream
_addedPDGCodes.clear(); m_addedPDGCodes.clear();
for ( iPDL = 0; iPDL < nPDL; iPDL++ ) { for ( iPDL = 0; iPDL < nPDL; iPDL++ ) {
EvtId particleId = EvtPDL::getEntry( iPDL ); EvtId particleId = EvtPDL::getEntry( iPDL );
int aliasInt = particleId.getAlias(); int aliasInt = particleId.getAlias();
// Pythia and EvtGen should use the same PDG codes for particles // Pythia and EvtGen should use the same PDG codes for particles
int PDGCode = EvtPDL::getStdHep( particleId ); int PDGCode = EvtPDL::getStdHep( particleId );
Show All 9 Lines for ( iPDL = 0; iPDL < nPDL; iPDL++ ) {
// the particle properties table. We can extract and change individual particle // the particle properties table. We can extract and change individual particle
// entries using the particleDataEntryPtr() function within ParticleData. // entries using the particleDataEntryPtr() function within ParticleData.
// However, we must be careful when accessing the particleData table. We must do // However, we must be careful when accessing the particleData table. We must do
// this directly, since assigning it to another Pythia8::ParticleData object via copying // this directly, since assigning it to another Pythia8::ParticleData object via copying
// or assignment will give it a different memory address and it will no longer refer to // or assignment will give it a different memory address and it will no longer refer to
// the original particleData information from the generator pointer. // the original particleData information from the generator pointer.
Pythia8::ParticleDataEntry* entry_generic = Pythia8::ParticleDataEntry* entry_generic =
_genericPythiaGen->particleData.particleDataEntryPtr( PDGCode ); m_genericPythiaGen->particleData.particleDataEntryPtr( PDGCode );
Pythia8::ParticleDataEntry* entry_alias = Pythia8::ParticleDataEntry* entry_alias =
_aliasPythiaGen->particleData.particleDataEntryPtr( PDGCode ); m_aliasPythiaGen->particleData.particleDataEntryPtr( PDGCode );
// Check that the PDG code is not zero/null and exclude other // Check that the PDG code is not zero/null and exclude other
// special cases, e.g. those reserved for internal generator use // special cases, e.g. those reserved for internal generator use
if ( entry_generic != 0 && this->validPDGCode( PDGCode ) ) { if ( entry_generic != 0 && this->validPDGCode( PDGCode ) ) {
entry_generic->setM0( mass ); entry_generic->setM0( mass );
entry_generic->setMWidth( width ); entry_generic->setMWidth( width );
entry_generic->setTau0( lifetime ); entry_generic->setTau0( lifetime );
Show All 22 Lines for ( iPDL = 0; iPDL < nPDL; iPDL++ ) {
bool hasPythiaDecays = EvtDecayTable::getInstance()->hasPythia( aliasInt ); bool hasPythiaDecays = EvtDecayTable::getInstance()->hasPythia( aliasInt );
if ( hasPythiaDecays ) { if ( hasPythiaDecays ) {
int isAlias = particleId.isAlias(); int isAlias = particleId.isAlias();
// Decide what generator to use depending on whether we have // Decide what generator to use depending on whether we have
// an aliased particle or not // an aliased particle or not
_thePythiaGenerator = ( isAlias == 1 ? _aliasPythiaGen.get() m_thePythiaGenerator = ( isAlias == 1 ? m_aliasPythiaGen.get()
: _genericPythiaGen.get() ); : m_genericPythiaGen.get() );
// Find the Pythia particle name given the standard PDG code integer // Find the Pythia particle name given the standard PDG code integer
std::string dataName = _thePythiaGenerator->particleData.name( std::string dataName = m_thePythiaGenerator->particleData.name(
PDGCode ); PDGCode );
bool alreadyStored = ( _addedPDGCodes.find( abs( PDGCode ) ) != bool alreadyStored = ( m_addedPDGCodes.find( abs( PDGCode ) ) !=
_addedPDGCodes.end() ); m_addedPDGCodes.end() );
if ( dataName == " " && !alreadyStored ) { if ( dataName == " " && !alreadyStored ) {
// Particle and its antiparticle do not exist in the Pythia database. // Particle and its antiparticle do not exist in the Pythia database.
// Create a new particle, then create the new decay modes. // Create a new particle, then create the new decay modes.
this->createPythiaParticle( particleId, PDGCode ); this->createPythiaParticle( particleId, PDGCode );
} }
// For the particle, create the Pythia decay modes. // For the particle, create the Pythia decay modes.
▲ Show 20 Lines • Show All 108 Lines • ▼ Show 20 Lines for ( iMode = 0; iMode < nModes; iMode++ ) {
int iDaug( 0 ); int iDaug( 0 );
for ( iDaug = 0; iDaug < nDaug; iDaug++ ) { for ( iDaug = 0; iDaug < nDaug; iDaug++ ) {
EvtId daugId = decayModel->getDaug( iDaug ); EvtId daugId = decayModel->getDaug( iDaug );
int daugPDG = EvtPDL::getStdHep( daugId ); int daugPDG = EvtPDL::getStdHep( daugId );
oss << " " << daugPDG; oss << " " << daugPDG;
} // Daughter list } // Daughter list
_thePythiaGenerator->readString( oss.str() ); m_thePythiaGenerator->readString( oss.str() );
} // is Pythia } // is Pythia
} else { } else {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Warning in EvtPythiaEngine. Trying to redefine Pythia table for " << "Warning in EvtPythiaEngine. Trying to redefine Pythia table for "
<< EvtPDL::name( particleId ) << EvtPDL::name( particleId )
<< " for a decay channel that has no daughters." << " for a decay channel that has no daughters."
<< " Keeping Pythia default (if available)." << endl; << " Keeping Pythia default (if available)." << endl;
} }
} else { } else {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Error in EvtPythiaEngine. decayModel is null for particle " << "Error in EvtPythiaEngine. decayModel is null for particle "
<< EvtPDL::name( particleId ) << " mode number " << iMode << EvtPDL::name( particleId ) << " mode number " << iMode
<< endl; << endl;
} }
} // Loop over modes } // Loop over modes
_pythiaModeMap[aliasInt] = pythiaModes; m_pythiaModeMap[aliasInt] = pythiaModes;
// Now, renormalise the decay branching fractions to sum to 1.0 // Now, renormalise the decay branching fractions to sum to 1.0
std::ostringstream rescaleStr; std::ostringstream rescaleStr;
rescaleStr.setf( std::ios::scientific ); rescaleStr.setf( std::ios::scientific );
rescaleStr << PDGCode << ":rescaleBR = 1.0"; rescaleStr << PDGCode << ":rescaleBR = 1.0";
_thePythiaGenerator->readString( rescaleStr.str() ); m_thePythiaGenerator->readString( rescaleStr.str() );
} }
int EvtPythiaEngine::getModeInt( EvtDecayBase* decayModel ) int EvtPythiaEngine::getModeInt( EvtDecayBase* decayModel )
{ {
int tmpModeInt( 0 ), modeInt( 0 ); int tmpModeInt( 0 ), modeInt( 0 );
if ( decayModel != 0 ) { if ( decayModel != 0 ) {
int nVars = decayModel->getNArg(); int nVars = decayModel->getNArg();
// Just read the first integer, which specifies the Pythia decay model. // Just read the first integer, which specifies the Pythia decay model.
// Ignore any other values. // Ignore any other values.
if ( nVars > 0 ) { if ( nVars > 0 ) {
tmpModeInt = static_cast<int>( decayModel->getArg( 0 ) ); tmpModeInt = static_cast<int>( decayModel->getArg( 0 ) );
} }
} }
if ( _convertPhysCodes ) { if ( m_convertPhysCodes ) {
// Extra code to convert the old Pythia decay model integer MDME(ICC,2) to the new one. // Extra code to convert the old Pythia decay model integer MDME(ICC,2) to the new one.
// This should be removed eventually after updating decay.dec files to use // This should be removed eventually after updating decay.dec files to use
// the new convention. // the new convention.
if ( tmpModeInt == 0 ) { if ( tmpModeInt == 0 ) {
modeInt = 0; // phase-space modeInt = 0; // phase-space
} else if ( tmpModeInt == 1 ) { } else if ( tmpModeInt == 1 ) {
modeInt = 1; // omega or phi -> 3pi modeInt = 1; // omega or phi -> 3pi
▲ Show 20 Lines • Show All 94 Lines • ▼ Show 20 Linesvoid EvtPythiaEngine::createPythiaParticle( EvtId& particleId, int PDGCode )
std::ostringstream oss; std::ostringstream oss;
oss.setf( std::ios::scientific ); oss.setf( std::ios::scientific );
int absPDGCode = abs( PDGCode ); int absPDGCode = abs( PDGCode );
oss << absPDGCode << ":new = " << aliasName << " " << antiName << " " oss << absPDGCode << ":new = " << aliasName << " " << antiName << " "
<< spin << " " << charge << " " << colour << " " << m0 << " " << mWidth << spin << " " << charge << " " << colour << " " << m0 << " " << mWidth
<< " " << mMin << " " << mMax << " " << tau0; << " " << mMin << " " << mMax << " " << tau0;
// Pass this information to Pythia // Pass this information to Pythia
_thePythiaGenerator->readString( oss.str() ); m_thePythiaGenerator->readString( oss.str() );
// Also store the absolute value of the PDG entry // Also store the absolute value of the PDG entry
// to keep track of which new particles have been added, // to keep track of which new particles have been added,
// which also automatically includes the anti-particle. // which also automatically includes the anti-particle.
// We need to avoid creating new anti-particles when // We need to avoid creating new anti-particles when
// they already exist when the particle was added. // they already exist when the particle was added.
_addedPDGCodes[absPDGCode] = 1; m_addedPDGCodes[absPDGCode] = 1;
} }
void EvtPythiaEngine::updatePhysicsParameters() void EvtPythiaEngine::updatePhysicsParameters()
{ {
// Update any more Pythia physics (or special particle) requirements/cuts etc.. // Update any more Pythia physics (or special particle) requirements/cuts etc..
// This should be used if any of the Pythia 6 parameters like JetSetPar MSTJ(i) = x // This should be used if any of the Pythia 6 parameters like JetSetPar MSTJ(i) = x
// are needed. Such commands will need to be implemented using the new interface // are needed. Such commands will need to be implemented using the new interface
// pythiaGenerator->readString(cmd); Here cmd is a string telling Pythia 8 // pythiaGenerator->readString(cmd); Here cmd is a string telling Pythia 8
// what physics parameters to change. This will need to be done for the generic and // what physics parameters to change. This will need to be done for the generic and
// alias generator pointers, as appropriate. // alias generator pointers, as appropriate.
// Set the multiplicity level for hadronic weak decays // Set the multiplicity level for hadronic weak decays
std::string multiWeakCut( "ParticleDecays:multIncreaseWeak = 2.0" ); std::string multiWeakCut( "ParticleDecays:multIncreaseWeak = 2.0" );
_genericPythiaGen->readString( multiWeakCut ); m_genericPythiaGen->readString( multiWeakCut );
_aliasPythiaGen->readString( multiWeakCut ); m_aliasPythiaGen->readString( multiWeakCut );
// Set the multiplicity level for all other decays // Set the multiplicity level for all other decays
std::string multiCut( "ParticleDecays:multIncrease = 4.5" ); std::string multiCut( "ParticleDecays:multIncrease = 4.5" );
_genericPythiaGen->readString( multiCut ); m_genericPythiaGen->readString( multiCut );
_aliasPythiaGen->readString( multiCut ); m_aliasPythiaGen->readString( multiCut );
//Now read in any custom configuration entered in the XML //Now read in any custom configuration entered in the XML
GeneratorCommands commands = GeneratorCommands commands =
EvtExtGeneratorCommandsTable::getInstance()->getCommands( "PYTHIA" ); EvtExtGeneratorCommandsTable::getInstance()->getCommands( "PYTHIA" );
GeneratorCommands::iterator it = commands.begin(); GeneratorCommands::iterator it = commands.begin();
for ( ; it != commands.end(); it++ ) { for ( ; it != commands.end(); it++ ) {
Command command = *it; Command command = *it;
Show All 24 Lines for ( ; it != commands.end(); it++ ) {
if ( generator == "GENERIC" || generator == "Generic" || if ( generator == "GENERIC" || generator == "Generic" ||
generator == "generic" || generator == "BOTH" || generator == "generic" || generator == "BOTH" ||
generator == "Both" || generator == "both" ) { generator == "Both" || generator == "both" ) {
std::vector<std::string>::iterator it2 = commandStrings.begin(); std::vector<std::string>::iterator it2 = commandStrings.begin();
for ( ; it2 != commandStrings.end(); it2++ ) { for ( ; it2 != commandStrings.end(); it2++ ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Configuring generic Pythia generator: " << ( *it2 ) << "Configuring generic Pythia generator: " << ( *it2 )
<< endl; << endl;
_genericPythiaGen->readString( *it2 ); m_genericPythiaGen->readString( *it2 );
} }
} }
if ( generator == "ALIAS" || generator == "Alias" || if ( generator == "ALIAS" || generator == "Alias" ||
generator == "alias" || generator == "BOTH" || generator == "alias" || generator == "BOTH" ||
generator == "Both" || generator == "both" ) { generator == "Both" || generator == "both" ) {
std::vector<std::string>::iterator it2 = commandStrings.begin(); std::vector<std::string>::iterator it2 = commandStrings.begin();
for ( ; it2 != commandStrings.end(); it2++ ) { for ( ; it2 != commandStrings.end(); it2++ ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" ) EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Configuring alias Pythia generator: " << ( *it2 ) << "Configuring alias Pythia generator: " << ( *it2 )
<< endl; << endl;
_aliasPythiaGen->readString( *it2 ); m_aliasPythiaGen->readString( *it2 );
} }
} }
} }
} }
#endif #endif