EvtPythiaEngine.cpp

Go to the documentation of this file.

/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
*                                                                      *
* This file is part of EvtGen.                                         *
*                                                                      *
* EvtGen is free software: you can redistribute it and/or modify       *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or    *
* (at your option) any later version.                                  *
*                                                                      *
* EvtGen is distributed in the hope that it will be useful,            *
* but WITHOUT ANY WARRANTY; without even the implied warranty of       *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
* GNU General Public License for more details.                         *
*                                                                      *
* You should have received a copy of the GNU General Public License    *
* along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
***********************************************************************/

#ifdef EVTGEN_PYTHIA

#include "EvtGenExternal/EvtPythiaEngine.hh"

#include "EvtGenBase/EvtDecayTable.hh"
#include "EvtGenBase/EvtExtGeneratorCommandsTable.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticleFactory.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtSpinType.hh"

#include "EvtGenExternal/EvtPythia6CommandConverter.hh"

#include "Pythia8/Event.h"
#include "Pythia8/ParticleData.h"

#include <cmath>
#include <iostream>
#include <sstream>

using std::endl;

EvtPythiaEngine::EvtPythiaEngine( std::string xmlDir, bool convertPhysCodes,
                                  bool useEvtGenRandom )
{
    // Create two Pythia generators. One will be for generic
    // Pythia decays in the decay.dec file. The other one will be to
    // only decay aliased particles, which are in general "signal"
    // decays different from those in the decay.dec file.
    // Even though it is not possible to have two different particle
    // versions in one Pythia generator, we can use two generators to
    // get the required behaviour.

    EvtGenReport( EVTGEN_INFO, "EvtGen" )
        << "Creating generic Pythia generator" << endl;
    _genericPythiaGen = std::make_unique<Pythia8::Pythia>( xmlDir );

    EvtGenReport( EVTGEN_INFO, "EvtGen" )
        << "Creating alias Pythia generator" << endl;
    _aliasPythiaGen = std::make_unique<Pythia8::Pythia>( xmlDir, false );

    _thePythiaGenerator = 0;
    _daugPDGVector.clear();
    _daugP4Vector.clear();

    _convertPhysCodes = convertPhysCodes;

    // Specify if we are going to use the random number generator (engine)
    // from EvtGen for Pythia 8.
    _useEvtGenRandom = useEvtGenRandom;

    _evtgenRandom = std::make_unique<EvtPythiaRandom>();

    _initialised = false;
}

EvtPythiaEngine::~EvtPythiaEngine()
{
    _thePythiaGenerator = nullptr;
    this->clearDaughterVectors();
    this->clearPythiaModeMap();
}

void EvtPythiaEngine::clearDaughterVectors()
{
    _daugPDGVector.clear();
    _daugP4Vector.clear();
}

void EvtPythiaEngine::clearPythiaModeMap()
{
    PythiaModeMap::iterator iter;
    for ( iter = _pythiaModeMap.begin(); iter != _pythiaModeMap.end(); ++iter ) {
        std::vector<int> modeVector = iter->second;
        modeVector.clear();
    }

    _pythiaModeMap.clear();
}

void EvtPythiaEngine::initialise()
{
    if ( _initialised ) {
        return;
    }

    this->clearPythiaModeMap();

    this->updateParticleLists();

    // 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..
    _genericPythiaGen->readString( "ProcessLevel:all = off" );
    _aliasPythiaGen->readString( "ProcessLevel:all = off" );

    // Turn off Pythia warnings, e.g. changes to particle properties
    _genericPythiaGen->readString( "Print:quiet = on" );
    _aliasPythiaGen->readString( "Print:quiet = on" );

    // Apply any other physics (or special particle) requirements/cuts etc..
    this->updatePhysicsParameters();

    // Set the random number generator
    if ( _useEvtGenRandom == true ) {
        _genericPythiaGen->setRndmEnginePtr( _evtgenRandom.get() );
        _aliasPythiaGen->setRndmEnginePtr( _evtgenRandom.get() );
    }

    _genericPythiaGen->init();
    _aliasPythiaGen->init();

    _initialised = true;
}

bool EvtPythiaEngine::doDecay( EvtParticle* theParticle )
{
    // Store the mother particle within a Pythia8 Event object.
    // Then do the hadron level decays.
    // 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
    // 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
    // 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.
    // 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,
    // 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
    // the EvtGen decay channel and the Pythia decay channel may be different.

    if ( _initialised == false ) {
        this->initialise();
    }

    if ( theParticle == 0 ) {
        EvtGenReport( EVTGEN_INFO, "EvtGen" )
            << "Error in EvtPythiaEngine::doDecay. The mother particle is null. Not doing any Pythia decay."
            << endl;
        return false;
    }

    // Delete EvtParticle daughters if they already exist
    if ( theParticle->getNDaug() != 0 ) {
        bool keepChannel( false );
        theParticle->deleteDaughters( keepChannel );
    }

    EvtId particleId = theParticle->getId();
    int isAlias = particleId.isAlias();

    // Choose the generator depending if we have an aliased (parent) particle or not
    _thePythiaGenerator = ( isAlias == 1 ? _aliasPythiaGen.get()
                                         : _genericPythiaGen.get() );

    // Need to use the reference to the Pythia8::Event object,
    // otherwise it will just return a new empty, default event object.
    Pythia8::Event& theEvent = _thePythiaGenerator->event;
    theEvent.reset();

    // Initialise the event to be the particle rest frame
    int PDGCode = EvtPDL::getStdHep( particleId );

    int status( 1 );
    int colour( 0 ), anticolour( 0 );
    double px( 0.0 ), py( 0.0 ), pz( 0.0 );
    double m0 = theParticle->mass();
    double E = m0;

    theEvent.append( PDGCode, status, colour, anticolour, px, py, pz, E, m0 );

    // Generate the Pythia event
    int iTrial( 0 );
    bool generatedEvent( false );
    for ( iTrial = 0; iTrial < 10; iTrial++ ) {
        generatedEvent = _thePythiaGenerator->next();
        if ( generatedEvent ) {
            break;
        }
    }

    bool success( false );

    if ( generatedEvent ) {
        // Store the daughters for this particle from the Pythia decay tree.
        // This is a recursive function that will continue looping through
        // all available daughters until the first set of non-quark and non-gluon
        // particles are encountered in the Pythia Event structure.

        // First, clear up the internal vectors storing the daughter
        // EvtId types and 4-momenta.
        this->clearDaughterVectors();

        // Now store the daughter info. Since this is a recursive function
        // to loop through the full Pythia decay tree, we do not want to create
        // EvtParticles here but in the next step.
        this->storeDaughterInfo( theParticle, 1 );

        // Now create the EvtParticle daughters of the (parent) particle.
        // We need to use the EvtParticle::makeDaughters function
        // owing to the way EvtParticle stores parent-daughter information.
        this->createDaughterEvtParticles( theParticle );

        //theParticle->printTree();
        //theEvent.list(true, true);

        success = true;
    }

    return success;
}

void EvtPythiaEngine::storeDaughterInfo( EvtParticle* theParticle, int startInt )
{
    Pythia8::Event& theEvent = _thePythiaGenerator->event;

    std::vector<int> daugList = theEvent.daughterList( startInt );

    std::vector<int>::iterator daugIter;
    for ( daugIter = daugList.begin(); daugIter != daugList.end(); ++daugIter ) {
        int daugInt = *daugIter;

        // Ask if the daughter is a quark or gluon. If so, recursively search again.
        Pythia8::Particle& daugParticle = theEvent[daugInt];

        if ( daugParticle.isQuark() || daugParticle.isGluon() ) {
            // Recursively search for correct daughter type
            this->storeDaughterInfo( theParticle, daugInt );

        } else {
            // We have a daughter that is not a quark nor gluon particle.
            // Make sure we are not double counting particles, since several quarks
            // and gluons make one particle.
            // Set the status flag for any "new" particle to say that we have stored it.
            // Use status flag = 1000 (within the user allowed range for Pythia codes).

            // Check that the status flag for the particle is not equal to 1000
            int statusCode = daugParticle.status();
            if ( statusCode != 1000 ) {
                int daugPDGInt = daugParticle.id();

                double px = daugParticle.px();
                double py = daugParticle.py();
                double pz = daugParticle.pz();
                double E = daugParticle.e();
                EvtVector4R daughterP4( E, px, py, pz );

                // Now store the EvtId and 4-momentum in the internal vectors
                _daugPDGVector.push_back( daugPDGInt );
                _daugP4Vector.push_back( daughterP4 );

                // Set the status flag for the Pythia particle to let us know
                // that we have already considered it to avoid double counting.
                daugParticle.status( 1000 );

            }    // Status code != 1000
        }
    }
}

void EvtPythiaEngine::createDaughterEvtParticles( EvtParticle* theParent )
{
    if ( theParent == 0 ) {
        EvtGenReport( EVTGEN_INFO, "EvtGen" )
            << "Error in EvtPythiaEngine::createDaughterEvtParticles. The parent is null"
            << endl;
        return;
    }

    // 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
    // for the particle decay, but this is not accessible from the Pythia interface at present.

    int nDaughters = _daugPDGVector.size();
    std::vector<EvtId> daugAliasIdVect( 0 );

    EvtId particleId = theParent->getId();
    // 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.
    int PDGId = EvtPDL::getStdHep( particleId );
    int aliasInt = particleId.getAlias();
    int pythiaAliasInt( aliasInt );

    if ( PDGId < 0 ) {
        // We have an anti-particle.
        EvtId conjPartId = EvtPDL::chargeConj( particleId );
        pythiaAliasInt = conjPartId.getAlias();
    }

    std::vector<int> pythiaModes = _pythiaModeMap[pythiaAliasInt];

    // 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.
    // This will then convert the Pythia generated channel to the EvtGen alias defined one.

    std::vector<int>::iterator modeIter;
    bool gotMode( false );

    for ( modeIter = pythiaModes.begin(); modeIter != pythiaModes.end();
          ++modeIter ) {
        // Stop the loop if we have the right decay mode channel
        if ( gotMode ) {
            break;
        }

        int pythiaModeInt = *modeIter;

        EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(
            aliasInt, pythiaModeInt );

        if ( decayModel != 0 ) {
            int nModeDaug = decayModel->getNDaug();

            // We need to make sure that the number of daughters match
            if ( nDaughters == nModeDaug ) {
                int iModeDaug( 0 );
                for ( iModeDaug = 0; iModeDaug < nModeDaug; iModeDaug++ ) {
                    EvtId daugId = decayModel->getDaug( iModeDaug );
                    int daugPDGId = EvtPDL::getStdHep( daugId );
                    // Pythia has used the right PDG codes for this decay mode, even for conjugate modes
                    int pythiaPDGId = _daugPDGVector[iModeDaug];

                    if ( daugPDGId == pythiaPDGId ) {
                        daugAliasIdVect.push_back( daugId );
                    }

                }    // Loop over EvtGen mode daughters

                int daugAliasSize = daugAliasIdVect.size();
                if ( daugAliasSize == nDaughters ) {
                    // All daughter Id codes are accounted for. Set the flag to stop the loop.
                    gotMode = true;
                } else {
                    // We do not have the correct daughter ordering. Clear the id vector
                    // and try another mode.
                    daugAliasIdVect.clear();
                }

            }    // Same number of daughters

        }    // decayModel != 0

    }    // Loop over available Pythia modes

    if ( gotMode == false ) {
        // We did not find a match for the daughter aliases. Just use the normal PDG codes
        // from the Pythia decay result
        int iPyDaug( 0 );
        for ( iPyDaug = 0; iPyDaug < nDaughters; iPyDaug++ ) {
            int daugPDGCode = _daugPDGVector[iPyDaug];
            EvtId daugPyId = EvtPDL::evtIdFromStdHep( daugPDGCode );
            daugAliasIdVect.push_back( daugPyId );
        }
    }

    // Make the EvtParticle daughters (with correct alias id's). Their 4-momenta are uninitialised.
    theParent->makeDaughters( nDaughters, daugAliasIdVect );

    // Now set the 4-momenta of the daughters.
    int iDaug( 0 );
    // Can use an iterator here, but we already had to use the vector size...
    for ( iDaug = 0; iDaug < nDaughters; iDaug++ ) {
        EvtParticle* theDaughter = theParent->getDaug( iDaug );

        // Set the correct 4-momentum for each daughter particle.
        if ( theDaughter != 0 ) {
            EvtId theDaugId = daugAliasIdVect[iDaug];
            const EvtVector4R theDaugP4 = _daugP4Vector[iDaug];
            theDaughter->init( theDaugId, theDaugP4 );
        }
    }
}

void EvtPythiaEngine::updateParticleLists()
{
    // Use the EvtGen decay table (decay/user.dec) to update Pythia particle
    // decay modes. Also, make sure the generic and alias Pythia generators
    // use the same particle data entries as defined by EvtGen (evt.pdl).

    // Loop over all entries in the EvtPDL particle data table.
    // Aliases are added at the end with id numbers equal to the
    // original particle, but with alias integer = lastPDLEntry+1 etc..
    int iPDL;
    int nPDL = EvtPDL::entries();

    // Reset the _addedPDGCodes map that keeps track
    // of any new particles added to the Pythia input data stream
    _addedPDGCodes.clear();

    for ( iPDL = 0; iPDL < nPDL; iPDL++ ) {
        EvtId particleId = EvtPDL::getEntry( iPDL );
        int aliasInt = particleId.getAlias();

        // Pythia and EvtGen should use the same PDG codes for particles
        int PDGCode = EvtPDL::getStdHep( particleId );

        // Update pole mass, width, lifetime and mass range
        double mass = EvtPDL::getMeanMass( particleId );
        double width = EvtPDL::getWidth( particleId );
        double lifetime = EvtPDL::getctau( particleId );
        double mmin = EvtPDL::getMinMass( particleId );
        double mmax = EvtPDL::getMaxMass( particleId );

        // Particle data pointers. The generic and alias Pythia generator pointers have
        // their own Pythia8::ParticleData particleData public data members which contain
        // the particle properties table. We can extract and change individual particle
        // entries using the particleDataEntryPtr() function within ParticleData.
        // 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
        // or assignment will give it a different memory address and it will no longer refer to
        // the original particleData information from the generator pointer.

        Pythia8::ParticleDataEntry* entry_generic =
            _genericPythiaGen->particleData.particleDataEntryPtr( PDGCode );
        Pythia8::ParticleDataEntry* entry_alias =
            _aliasPythiaGen->particleData.particleDataEntryPtr( PDGCode );

        // Check that the PDG code is not zero/null and exclude other
        // special cases, e.g. those reserved for internal generator use
        if ( entry_generic != 0 && this->validPDGCode( PDGCode ) ) {
            entry_generic->setM0( mass );
            entry_generic->setMWidth( width );
            entry_generic->setTau0( lifetime );

            if ( std::fabs( width ) > 0.0 ) {
                entry_generic->setMMin( mmin );
                entry_generic->setMMax( mmax );
            }
        }

        // Check that the PDG code is not zero/null and exclude other
        // special cases, e.g. those reserved for internal generator use
        if ( entry_alias != 0 && this->validPDGCode( PDGCode ) ) {
            entry_alias->setM0( mass );
            entry_alias->setMWidth( width );
            entry_alias->setTau0( lifetime );

            if ( std::fabs( width ) > 0.0 ) {
                entry_alias->setMMin( mmin );
                entry_alias->setMMax( mmax );
            }
        }

        // Check which particles have a Pythia decay defined.
        // Get the list of all possible decays for the particle, using the alias integer.
        // If the particle is not actually an alias, aliasInt = idInt.

        bool hasPythiaDecays = EvtDecayTable::getInstance()->hasPythia( aliasInt );

        if ( hasPythiaDecays ) {
            int isAlias = particleId.isAlias();

            // Decide what generator to use depending on whether we have
            // an aliased particle or not
            _thePythiaGenerator = ( isAlias == 1 ? _aliasPythiaGen.get()
                                                 : _genericPythiaGen.get() );

            // Find the Pythia particle name given the standard PDG code integer
            std::string dataName = _thePythiaGenerator->particleData.name(
                PDGCode );
            bool alreadyStored = ( _addedPDGCodes.find( abs( PDGCode ) ) !=
                                   _addedPDGCodes.end() );

            if ( dataName == " " && !alreadyStored ) {
                // Particle and its antiparticle do not exist in the Pythia database.
                // Create a new particle, then create the new decay modes.
                this->createPythiaParticle( particleId, PDGCode );
            }

            // For the particle, create the Pythia decay modes.
            // Update Pythia data tables.
            this->updatePythiaDecayTable( particleId, aliasInt, PDGCode );

        }    // Loop over Pythia decays

    }    // Loop over EvtPDL entries

    //EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Writing out changed generic Pythia decay list"<<endl;
    //_genericPythiaGen->particleData.listChanged();

    //EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Writing out changed alias Pythia decay list"<<endl;
    //_aliasPythiaGen->particleData.listChanged();
}

bool EvtPythiaEngine::validPDGCode( int PDGCode )
{
    // Exclude certain PDG codes: void = 0 and special values = 81 to 100, which are reserved
    // for internal generator use (pseudoparticles) according to PDG guidelines. Also exclude
    // nu'_tau (nu_L) = 18, which has different masses: Pythia8 = 400 GeV, EvtGen = 0 GeV.

    bool isValid( true );

    int absPDGCode = abs( PDGCode );

    if ( absPDGCode == 0 || absPDGCode == 18 ) {
        // Void and nu_L or nu'_tau
        isValid = false;
    } else if ( absPDGCode >= 81 && absPDGCode <= 100 ) {
        // Pseudoparticles
        isValid = false;
    }

    return isValid;
}

void EvtPythiaEngine::updatePythiaDecayTable( EvtId& particleId, int aliasInt,
                                              int PDGCode )
{
    // Update the particle data table in Pythia.
    // The tables store information about the allowed decay modes
    // where the PDGId for all particles must be positive; anti-particles are stored
    // with the corresponding particle entry.
    // Since we do not want to implement CP violation here, just use the same branching
    // fractions for particle and anti-particle modes.

    int nModes = EvtDecayTable::getInstance()->getNModes( aliasInt );
    int iMode( 0 );

    bool firstMode( true );

    // Only process positive PDG codes.
    if ( PDGCode < 0 ) {
        return;
    }

    // Keep track of which decay modes are Pythia decays for each aliasInt
    std::vector<int> pythiaModes( 0 );

    // Loop over the decay modes for this particle
    for ( iMode = 0; iMode < nModes; iMode++ ) {
        EvtDecayBase* decayModel =
            EvtDecayTable::getInstance()->findDecayModel( aliasInt, iMode );

        if ( decayModel != 0 ) {
            int nDaug = decayModel->getNDaug();

            // If the decay mode has no daughters, then that means that there will be
            // no entries for any submode re-definitions for Pythia.
            // This sometimes occurs for any mode using non-standard Pythia 6 codes.
            // Do not refine the decay mode, i.e. accept the Pythia 8 default (if it exists).
            if ( nDaug > 0 ) {
                // Check to see if we have a Pythia decay mode
                std::string modelName = decayModel->getModelName();

                if ( modelName == "PYTHIA" ) {
                    // Keep track which decay mode is a Pythia one. We need this in order to
                    // reassign alias Id values for particles generated in the decay.
                    pythiaModes.push_back( iMode );

                    std::ostringstream oss;
                    oss.setf( std::ios::scientific );
                    // Write out the absolute value of the PDG code, since
                    // particles and anti-particles occupy the same part of the Pythia table.
                    oss << PDGCode;

                    if ( firstMode ) {
                        // Create a new channel
                        oss << ":oneChannel = ";
                        firstMode = false;
                    } else {
                        // Add the channel
                        oss << ":addChannel = ";
                    }

                    // Select all channels (particle and anti-particle).
                    // For CP violation, or different BFs for particle and anti-particle,
                    // use options 2 or 3 (not here).
                    int onMode( 1 );
                    oss << onMode << " ";

                    double BF = decayModel->getBranchingFraction();
                    oss << BF << " ";

                    // Need to convert the old Pythia physics mode integers with the new ones
                    // To do this, get the model argument and write a conversion method.
                    int modeInt = this->getModeInt( decayModel );
                    oss << modeInt;

                    int iDaug( 0 );
                    for ( iDaug = 0; iDaug < nDaug; iDaug++ ) {
                        EvtId daugId = decayModel->getDaug( iDaug );
                        int daugPDG = EvtPDL::getStdHep( daugId );
                        oss << " " << daugPDG;

                    }    // Daughter list

                    _thePythiaGenerator->readString( oss.str() );

                }    // is Pythia

            } else {
                EvtGenReport( EVTGEN_INFO, "EvtGen" )
                    << "Warning in EvtPythiaEngine. Trying to redefine Pythia table for "
                    << EvtPDL::name( particleId )
                    << " for a decay channel that has no daughters."
                    << " Keeping Pythia default (if available)." << endl;
            }

        } else {
            EvtGenReport( EVTGEN_INFO, "EvtGen" )
                << "Error in EvtPythiaEngine. decayModel is null for particle "
                << EvtPDL::name( particleId ) << " mode number " << iMode
                << endl;
        }

    }    // Loop over modes

    _pythiaModeMap[aliasInt] = pythiaModes;

    // Now, renormalise the decay branching fractions to sum to 1.0
    std::ostringstream rescaleStr;
    rescaleStr.setf( std::ios::scientific );
    rescaleStr << PDGCode << ":rescaleBR = 1.0";

    _thePythiaGenerator->readString( rescaleStr.str() );
}

int EvtPythiaEngine::getModeInt( EvtDecayBase* decayModel )
{
    int tmpModeInt( 0 ), modeInt( 0 );

    if ( decayModel != 0 ) {
        int nVars = decayModel->getNArg();
        // Just read the first integer, which specifies the Pythia decay model.
        // Ignore any other values.
        if ( nVars > 0 ) {
            tmpModeInt = static_cast<int>( decayModel->getArg( 0 ) );
        }
    }

    if ( _convertPhysCodes ) {
        // 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
        // the new convention.

        if ( tmpModeInt == 0 ) {
            modeInt = 0;    // phase-space
        } else if ( tmpModeInt == 1 ) {
            modeInt = 1;    // omega or phi -> 3pi
        } else if ( tmpModeInt == 2 ) {
            modeInt = 11;    // Dalitz decay
        } else if ( tmpModeInt == 3 ) {
            modeInt = 2;    // V -> PS PS
        } else if ( tmpModeInt == 4 ) {
            modeInt = 92;    // onium -> ggg or gg gamma
        } else if ( tmpModeInt == 11 ) {
            modeInt = 42;    // phase-space of hadrons from available quarks
        } else if ( tmpModeInt == 12 ) {
            modeInt = 42;    // phase-space for onia resonances
        } else if ( tmpModeInt == 13 ) {
            modeInt = 43;    // phase-space of at least 3 hadrons
        } else if ( tmpModeInt == 14 ) {
            modeInt = 44;    // phase-space of at least 4 hadrons
        } else if ( tmpModeInt == 15 ) {
            modeInt = 45;    // phase-space of at least 5 hadrons
        } else if ( tmpModeInt >= 22 && tmpModeInt <= 30 ) {
            modeInt = tmpModeInt +
                      40;    // phase space of hadrons with fixed multiplicity (modeInt - 60)
        } else if ( tmpModeInt == 31 ) {
            modeInt = 42;    // two or more quarks phase-space; one spectactor quark
        } else if ( tmpModeInt == 32 ) {
            modeInt = 91;    // qqbar or gg pair
        } else if ( tmpModeInt == 33 ) {
            modeInt = 0;    // triplet q X qbar, where X = gluon or colour singlet (superfluous, since g's are created anyway)
        } else if ( tmpModeInt == 41 ) {
            modeInt = 21;    // weak decay phase space, weighting nu_tau spectrum
        } else if ( tmpModeInt == 42 ) {
            modeInt = 22;    // weak decay V-A matrix element
        } else if ( tmpModeInt == 43 ) {
            modeInt = 22;    // weak decay V-A matrix element, quarks as jets (superfluous)
        } else if ( tmpModeInt == 44 ) {
            modeInt = 22;    // weak decay V-A matrix element, parton showers (superfluous)
        } else if ( tmpModeInt == 48 ) {
            modeInt = 23;    // weak decay V-A matrix element, at least 3 decay products
        } else if ( tmpModeInt == 50 ) {
            modeInt = 0;    // default behaviour
        } else if ( tmpModeInt == 51 ) {
            modeInt = 0;    // step threshold (channel switched off when mass daughters > mother mass
        } else if ( tmpModeInt == 52 || tmpModeInt == 53 ) {
            modeInt = 0;    // beta-factor threshold
        } else if ( tmpModeInt == 84 ) {
            modeInt = 42;    // unknown physics process - just use phase-space
        } else if ( tmpModeInt == 101 ) {
            modeInt = 0;    // continuation line
        } else if ( tmpModeInt == 102 ) {
            modeInt = 0;    // off mass shell particles.
        } else {
            EvtGenReport( EVTGEN_INFO, "EvtGen" )
                << "Pythia mode integer " << tmpModeInt
                << " is not recognised. Using phase-space model" << endl;
            modeInt = 0;    // Use phase-space for anything else
        }

    } else {
        // No need to convert the physics mode integer code
        modeInt = tmpModeInt;
    }

    return modeInt;
}

void EvtPythiaEngine::createPythiaParticle( EvtId& particleId, int PDGCode )
{
    // Use the EvtGen name, PDGId and other variables to define the new Pythia particle.
    EvtId antiPartId = EvtPDL::chargeConj( particleId );

    std::string aliasName = EvtPDL::name(
        particleId );    // If not an alias, aliasName = normal name
    std::string antiName = EvtPDL::name( antiPartId );

    EvtSpinType::spintype spinType = EvtPDL::getSpinType( particleId );
    int spin = EvtSpinType::getSpin2( spinType );

    int charge = EvtPDL::chg3( particleId );

    // Must set the correct colour type manually here, since the evt.pdl file
    // does not store this information. This is required for quarks otherwise
    // Pythia cannot generate the decay properly.
    int PDGId = EvtPDL::getStdHep( particleId );
    int colour( 0 );
    if ( PDGId == 21 ) {
        colour = 2;    // gluons
    } else if ( PDGId <= 8 && PDGId > 0 ) {
        colour = 1;    // single quarks
    }

    double m0 = EvtPDL::getMeanMass( particleId );
    double mWidth = EvtPDL::getWidth( particleId );
    double mMin = EvtPDL::getMinMass( particleId );
    double mMax = EvtPDL::getMaxMass( particleId );

    double tau0 = EvtPDL::getctau( particleId );

    std::ostringstream oss;
    oss.setf( std::ios::scientific );
    int absPDGCode = abs( PDGCode );
    oss << absPDGCode << ":new = " << aliasName << " " << antiName << " "
        << spin << " " << charge << " " << colour << " " << m0 << " " << mWidth
        << " " << mMin << " " << mMax << " " << tau0;

    // Pass this information to Pythia
    _thePythiaGenerator->readString( oss.str() );

    // Also store the absolute value of the PDG entry
    // to keep track of which new particles have been added,
    // which also automatically includes the anti-particle.
    // We need to avoid creating new anti-particles when
    // they already exist when the particle was added.
    _addedPDGCodes[absPDGCode] = 1;
}

void EvtPythiaEngine::updatePhysicsParameters()
{
    // 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
    // are needed. Such commands will need to be implemented using the new interface
    // 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
    // alias generator pointers, as appropriate.

    // Set the multiplicity level for hadronic weak decays
    std::string multiWeakCut( "ParticleDecays:multIncreaseWeak = 2.0" );
    _genericPythiaGen->readString( multiWeakCut );
    _aliasPythiaGen->readString( multiWeakCut );

    // Set the multiplicity level for all other decays
    std::string multiCut( "ParticleDecays:multIncrease = 4.5" );
    _genericPythiaGen->readString( multiCut );
    _aliasPythiaGen->readString( multiCut );

    //Now read in any custom configuration entered in the XML
    GeneratorCommands commands =
        EvtExtGeneratorCommandsTable::getInstance()->getCommands( "PYTHIA" );
    GeneratorCommands::iterator it = commands.begin();

    for ( ; it != commands.end(); it++ ) {
        Command command = *it;
        std::vector<std::string> commandStrings;

        if ( command["VERSION"] == "PYTHIA6" ) {
            EvtGenReport( EVTGEN_INFO, "EvtGen" )
                << "Converting Pythia 6 command: " << command["MODULE"] << "("
                << command["PARAM"] << ")=" << command["VALUE"] << "..." << endl;
            commandStrings = convertPythia6Command( command );
        } else if ( command["VERSION"] == "PYTHIA8" ) {
            commandStrings.push_back( command["MODULE"] + ":" + command["PARAM"] +
                                      " = " + command["VALUE"] );
        } else {
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "Pythia command received by EvtPythiaEngine has bad version:"
                << endl;
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "Received " << command["VERSION"]
                << " but expected PYTHIA6 or PYTHIA8." << endl;
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "The error is likely to be in EvtDecayTable.cpp" << endl;
            EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                << "EvtGen will now abort." << endl;
            ::abort();
        }
        std::string generator = command["GENERATOR"];
        if ( generator == "GENERIC" || generator == "Generic" ||
             generator == "generic" || generator == "BOTH" ||
             generator == "Both" || generator == "both" ) {
            std::vector<std::string>::iterator it2 = commandStrings.begin();
            for ( ; it2 != commandStrings.end(); it2++ ) {
                EvtGenReport( EVTGEN_INFO, "EvtGen" )
                    << "Configuring generic Pythia generator: " << ( *it2 )
                    << endl;
                _genericPythiaGen->readString( *it2 );
            }
        }
        if ( generator == "ALIAS" || generator == "Alias" ||
             generator == "alias" || generator == "BOTH" ||
             generator == "Both" || generator == "both" ) {
            std::vector<std::string>::iterator it2 = commandStrings.begin();
            for ( ; it2 != commandStrings.end(); it2++ ) {
                EvtGenReport( EVTGEN_INFO, "EvtGen" )
                    << "Configuring alias Pythia generator: " << ( *it2 )
                    << endl;
                _aliasPythiaGen->readString( *it2 );
            }
        }
    }
}

#endif