diff --git a/src/EvtGenExternal/EvtPythiaEngine.cpp b/src/EvtGenExternal/EvtPythiaEngine.cpp index df38ec8..72aa46d 100644 --- a/src/EvtGenExternal/EvtPythiaEngine.cpp +++ b/src/EvtGenExternal/EvtPythiaEngine.cpp @@ -1,847 +1,853 @@ /*********************************************************************** * 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 . * ***********************************************************************/ #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 #include #include +#if PYTHIA_VERSION_INTEGER < 8304 +typedef Pythia8::ParticleDataEntry* ParticleDataEntryPtr; +#else +typedef Pythia8::ParticleDataEntryPtr ParticleDataEntryPtr; +#endif + 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( xmlDir ); EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Creating alias Pythia generator" << endl; _aliasPythiaGen = std::make_unique( 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(); _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 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 daugList = theEvent.daughterList( startInt ); std::vector::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 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 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::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 = + ParticleDataEntryPtr entry_generic = _genericPythiaGen->particleData.particleDataEntryPtr( PDGCode ); - Pythia8::ParticleDataEntry* entry_alias = + ParticleDataEntryPtr 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"<particleData.listChanged(); //EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Writing out changed alias Pythia decay list"<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 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( 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 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::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::iterator it2 = commandStrings.begin(); for ( ; it2 != commandStrings.end(); it2++ ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Configuring alias Pythia generator: " << ( *it2 ) << endl; _aliasPythiaGen->readString( *it2 ); } } } } #endif