diff --git a/FastSimulation/ParticleDecay/plugins/TestPythiaDecays.cc b/FastSimulation/ParticleDecay/plugins/TestPythiaDecays.cc index e3849f3541294..9859714579d58 100644 --- a/FastSimulation/ParticleDecay/plugins/TestPythiaDecays.cc +++ b/FastSimulation/ParticleDecay/plugins/TestPythiaDecays.cc @@ -76,7 +76,7 @@ class TestPythiaDecays : public edm::stream::EDAnalyzer <> { std::map h_br; std::map h_br_ref; - std::map > knownDecayModes; + std::map > knownDecayModes; Pythia8::Pythia * pythia; std::string outputFile; @@ -182,7 +182,7 @@ TestPythiaDecays::TestPythiaDecays(const edm::ParameterSet& iConfig) h_br[pid]->SetCanExtend(TH1::kAllAxes); h_br_ref[pid] = (TH1D*)(h_br[pid]->Clone(strstr.str().c_str())); h_br_ref[pid]->SetTitle(h_br_ref[pid]->GetName()); - knownDecayModes[pid] = vector(); + knownDecayModes[pid] = std::vector(); for(int d = 0;dsizeChannels();++d){ Pythia8::DecayChannel & channel = pd->channel(d); std::vector prod; @@ -334,7 +334,7 @@ TestPythiaDecays::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetu size_t parentIndex = it->first; const SimTrack & parent = simtracks->at(parentIndex); int pid = abs(parent.type()); - vector & childIndices = it->second; + std::vector & childIndices = it->second; if(childIndices.size() == 0) continue; diff --git a/GeneratorInterface/Pythia8Interface/BuildFile.xml b/GeneratorInterface/Pythia8Interface/BuildFile.xml index 6d3d01f253105..048d63026c1a3 100644 --- a/GeneratorInterface/Pythia8Interface/BuildFile.xml +++ b/GeneratorInterface/Pythia8Interface/BuildFile.xml @@ -6,6 +6,7 @@ + diff --git a/GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h b/GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h index c676439e9dbef..697c2ec694d61 100644 --- a/GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h +++ b/GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h @@ -12,6 +12,8 @@ #include #include +class EvtGenDecays; + namespace CLHEP { class HepRandomEngine; } @@ -51,6 +53,13 @@ namespace gen { unsigned int maxEventsToPrint; HepMC::IO_AsciiParticles* ascii_io; + // EvtGen plugin + // + bool useEvtGen; + EvtGenDecays* evtgenDecays; + std::string evtgenDecFile; + std::string evtgenPdlFile; + private: P8RndmEngine p8RndmEngine_; diff --git a/GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.cc b/GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.cc index 483a1e62faad4..8a374d1842749 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.cc @@ -1,9 +1,13 @@ +#include "Pythia8/Pythia.h" + +using namespace Pythia8; + #include "GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h" #include "FWCore/ServiceRegistry/interface/Service.h" -void EmissionVetoHook1::fatalEmissionVeto(string message) { +void EmissionVetoHook1::fatalEmissionVeto(std::string message) { throw edm::Exception(edm::errors::Configuration,"Pythia8Interface") - << "EmissionVeto: " << message << endl; + << "EmissionVeto: " << message << std::endl; } //-------------------------------------------------------------------------- @@ -56,14 +60,14 @@ double EmissionVetoHook1::pTpythia(const Pythia8::Event &e, int RadAfterBranch, // Can get negative pT for massive splittings if (pTnow < 0.) { - cout << "Warning: pTpythia was negative" << endl; + std::cout << "Warning: pTpythia was negative" << std::endl; return -1.; } #ifdef DBGOUTPUT - cout << "pTpythia: rad = " << RadAfterBranch << ", emt = " + std::cout << "pTpythia: rad = " << RadAfterBranch << ", emt = " << EmtAfterBranch << ", rec = " << RecAfterBranch - << ", pTnow = " << sqrt(pTnow) << endl; + << ", pTnow = " << sqrt(pTnow) << std::endl; #endif // Return pT @@ -97,13 +101,13 @@ double EmissionVetoHook1::pTpowheg(const Pythia8::Event &e, int i, int j, bool F // Check result if (pTnow < 0.) { - cout << "Warning: pTpowheg was negative" << endl; + std::cout << "Warning: pTpowheg was negative" << std::endl; return -1.; } #ifdef DBGOUTPUT - cout << "pTpowheg: i = " << i << ", j = " << j - << ", pTnow = " << pTnow << endl; + std::cout << "pTpowheg: i = " << i << ", j = " << j + << ", pTnow = " << pTnow << std::endl; #endif return pTnow; @@ -163,7 +167,7 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i // ISR - only done once as just kinematical pT if (!FSR) { pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR); - if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow); + if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow); // FSR - try all outgoing partons from system before branching // as i. Note that for the hard system, there is no @@ -182,7 +186,7 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0) ? false : FSR); if (pTnow > 0.) pTemt = (pTemt < 0) - ? pTnow : min(pTemt, pTnow); + ? pTnow : std::min(pTemt, pTnow); } // for (iMem) } // if (!FSR) @@ -193,9 +197,9 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i // ISR - other incoming as recoiler if (!FSR) { pTnow = pTpythia(e, iInA, jNow, iInB, FSR); - if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow); + if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow); pTnow = pTpythia(e, iInB, jNow, iInA, FSR); - if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow); + if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow); // FSR - try all final-state coloured partons as radiator // after emission (k). @@ -208,10 +212,10 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i // Try two incoming. pTnow = pTpythia(e, kNow, jNow, iInA, FSR); if (pTnow > 0.) pTemt = (pTemt < 0) - ? pTnow : min(pTemt, pTnow); + ? pTnow : std::min(pTemt, pTnow); pTnow = pTpythia(e, kNow, jNow, iInB, FSR); if (pTnow > 0.) pTemt = (pTemt < 0) - ? pTnow : min(pTemt, pTnow); + ? pTnow : std::min(pTemt, pTnow); // Try all other outgoing. for (int rNow = 0; rNow < e.size(); rNow++) { @@ -219,7 +223,7 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i !e[rNow].isFinal() || e[rNow].colType() == 0) continue; pTnow = pTpythia(e, kNow, jNow, rNow, FSR); if (pTnow > 0.) pTemt = (pTemt < 0) - ? pTnow : min(pTemt, pTnow); + ? pTnow : std::min(pTemt, pTnow); } // for (rNow) } // for (kNow) @@ -230,9 +234,9 @@ double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, i } // for (xSR) #ifdef DBGOUTPUT - cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k + std::cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k << ", r = " << r << ", xSR = " << xSRin - << ", pTemt = " << pTemt << endl; + << ", pTemt = " << pTemt << std::endl; #endif return pTemt; @@ -274,7 +278,7 @@ bool EmissionVetoHook1::doVetoMPIStep(int nMPI, const Pythia8::Event &e) { first = ip; break; } - if(first < 0) fatalEmissionVeto(string("signal particles not found")); + if(first < 0) fatalEmissionVeto(std::string("signal particles not found")); for(int ip = first; ip < e.size(); ip++) { myid = e[ip].id(); if(abs(myid) < 6 || abs(myid) == 21) continue; @@ -285,7 +289,7 @@ bool EmissionVetoHook1::doVetoMPIStep(int nMPI, const Pythia8::Event &e) { // Extra check that we have the correct final state if (count != nFinal && count != nFinal + 1) - fatalEmissionVeto(string("Wrong number of final state particles in event")); + fatalEmissionVeto(std::string("Wrong number of final state particles in event")); // Flag if POWHEG radiation present and index bool isEmt = (count == nFinal) ? false : true; @@ -314,9 +318,9 @@ bool EmissionVetoHook1::doVetoMPIStep(int nMPI, const Pythia8::Event &e) { } if(Verbosity) - cout << "doVetoMPIStep: QFac = " << infoPtr->QFac() + std::cout << "doVetoMPIStep: QFac = " << infoPtr->QFac() << ", QRen = " << infoPtr->QRen() - << ", pThard = " << pThard << endl << endl; + << ", pThard = " << pThard << std::endl << std::endl; // Initialise other variables accepted = false; @@ -347,12 +351,12 @@ bool EmissionVetoHook1::doVetoISREmission(int, const Pythia8::Event &e, int iSys } if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) { e.list(); - fatalEmissionVeto(string("Couldn't find Pythia ISR emission")); + fatalEmissionVeto(std::string("Couldn't find Pythia ISR emission")); } // pTemtMode == 0: pT of emitted w.r.t. radiator - // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing) - // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing) + // pTemtMode == 1: std::min(pT of emitted w.r.t. all incoming/outgoing) + // pTemtMode == 2: std::min(pT of all outgoing w.r.t. all incoming/outgoing) int xSR = (pTemtMode == 0) ? 0 : -1; int i = (pTemtMode == 0) ? iRadAft : -1; int j = (pTemtMode != 2) ? iEmt : -1; @@ -361,7 +365,7 @@ bool EmissionVetoHook1::doVetoISREmission(int, const Pythia8::Event &e, int iSys double pTemt = pTcalc(e, i, j, k, r, xSR); #ifdef DBGOUTPUT - cout << "doVetoISREmission: pTemt = " << pTemt << endl << endl; + std::cout << "doVetoISREmission: pTemt = " << pTemt << std::endl << std::endl; #endif // Veto if pTemt > pThard @@ -396,13 +400,13 @@ bool EmissionVetoHook1::doVetoFSREmission(int, const Pythia8::Event &e, int iSys if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) || e[iEmt].status() != 51 || e[iRadAft].status() != 51) { e.list(); - fatalEmissionVeto(string("Couldn't find Pythia FSR emission")); + fatalEmissionVeto(std::string("Couldn't find Pythia FSR emission")); } // Behaviour based on pTemtMode: // 0 - pT of emitted w.r.t. radiator before - // 1 - min(pT of emitted w.r.t. all incoming/outgoing) - // 2 - min(pT of all outgoing w.r.t. all incoming/outgoing) + // 1 - std::min(pT of emitted w.r.t. all incoming/outgoing) + // 2 - std::min(pT of all outgoing w.r.t. all incoming/outgoing) int xSR = (pTemtMode == 0) ? 1 : -1; int i = (pTemtMode == 0) ? iRadBef : -1; int k = (pTemtMode == 0) ? iRadAft : -1; @@ -421,11 +425,11 @@ bool EmissionVetoHook1::doVetoFSREmission(int, const Pythia8::Event &e, int iSys for (int jLoop = 0; jLoop < 2; jLoop++) { if (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR); - else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR)); + else if (jLoop == 1) pTemt = std::min(pTemt, pTcalc(e, i, j, k, r, xSR)); // For emittedMode == 3, have tried iRadAft, now try iEmt if (emittedMode != 3) break; - if (k != -1) swap(j, k); else j = iEmt; + if (k != -1) std::swap(j, k); else j = iEmt; } // If pTemtMode is 2, then try all final-state partons as emitted @@ -435,7 +439,7 @@ bool EmissionVetoHook1::doVetoFSREmission(int, const Pythia8::Event &e, int iSys } #ifdef DBGOUTPUT - cout << "doVetoFSREmission: pTemt = " << pTemt << endl << endl; + std::cout << "doVetoFSREmission: pTemt = " << pTemt << std::endl << std::endl; #endif // Veto if pTemt > pThard @@ -459,8 +463,8 @@ bool EmissionVetoHook1::doVetoMPIEmission(int, const Pythia8::Event &e) { if (MPIvetoOn) { if (e[e.size() - 1].pT() > pTMPI) { #ifdef DBGOUTPUT - cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT() - << ", pTMPI = " << pTMPI << endl << endl; + std::cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT() + << ", pTMPI = " << pTMPI << std::endl << std::endl; #endif return true; } diff --git a/GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h b/GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h index 8de630824ca59..9e7bd579599f8 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h +++ b/GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h @@ -1,5 +1,3 @@ -#include "Pythia8/Pythia.h" - class EmissionVetoHook1 : public Pythia8::UserHooks { public: @@ -15,8 +13,8 @@ class EmissionVetoHook1 : public Pythia8::UserHooks { MPIvetoOn(MPIvetoOnIn), nISRveto(0), nFSRveto(0), Verbosity(VerbosityIn) {} ~EmissionVetoHook1() { - cout << "Number of ISR vetoed = " << nISRveto << endl; - cout << "Number of FSR vetoed = " << nFSRveto << endl; + std::cout << "Number of ISR vetoed = " << nISRveto << std::endl; + std::cout << "Number of FSR vetoed = " << nFSRveto << std::endl; } //-------------------------------------------------------------------------- @@ -34,7 +32,7 @@ class EmissionVetoHook1 : public Pythia8::UserHooks { bool canVetoMPIEmission() { return MPIvetoOn; } bool doVetoMPIEmission(int, const Pythia8::Event &e); - void fatalEmissionVeto(string message); + void fatalEmissionVeto(std::string message); double pTpythia(const Pythia8::Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR); diff --git a/GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.cc b/GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.cc index 8f3089b6e7d6d..912c70be82f40 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.cc @@ -42,10 +42,13 @@ bool LHAupLesHouches::setInit() infoPtr->setHeader("slha",slhaheader); } - //work around missing initialization inside pythia8 - infoPtr->eventAttributes = new std::map; - - + //will be used to work around missing initialization inside pythia8 + if(!fEvAttributes) { + fEvAttributes = new std::map; + } else { + fEvAttributes->clear(); + } + return true; } @@ -91,8 +94,13 @@ bool LHAupLesHouches::setEvent(int inProcId) hepeup.SPINUP[i],scalein); } - infoPtr->eventAttributes->clear(); - + if(!infoPtr->eventAttributes) { + fEvAttributes->clear(); + infoPtr->eventAttributes = fEvAttributes; + } else { + infoPtr->eventAttributes->clear(); + } + //fill parton multiplicities if available int npLO = event->npLO(); int npNLO = event->npNLO(); diff --git a/GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.h b/GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.h index ef3cc3c71b578..8d6eb5f2e0961 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.h +++ b/GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.h @@ -21,7 +21,7 @@ class LHAupLesHouches : public Pythia8::LHAup { public: - LHAupLesHouches() : setScalesFromLHEF_(false) {;} + LHAupLesHouches() : setScalesFromLHEF_(false),fEvAttributes(0) {;} //void loadRunInfo(const boost::shared_ptr &runInfo) void loadRunInfo(lhef::LHERunInfo* runInfo) @@ -33,6 +33,8 @@ class LHAupLesHouches : public Pythia8::LHAup { void setScalesFromLHEF(bool b) { setScalesFromLHEF_ = b; } + ~LHAupLesHouches() {if(fEvAttributes) delete fEvAttributes;} + private: bool setInit(); @@ -46,4 +48,5 @@ class LHAupLesHouches : public Pythia8::LHAup { // Flag to set particle production scales or not. bool setScalesFromLHEF_; + std::map * fEvAttributes; }; diff --git a/GeneratorInterface/Pythia8Interface/plugins/Py8EGun.cc b/GeneratorInterface/Pythia8Interface/plugins/Py8EGun.cc index 85e89b42f57d7..531bcbb7346a1 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/Py8EGun.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/Py8EGun.cc @@ -4,6 +4,10 @@ #include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h" +// EvtGen plugin +// +#include "Pythia8Plugins/EvtGen.h" + namespace gen { class Py8EGun : public Py8GunBase { @@ -80,11 +84,11 @@ bool Py8EGun::generatePartonsAndHadronize() { if ( (fMasterGen->particleData).isParticle( -particleID ) ) { - (fMasterGen->event).append( -particleID, 1, 0, 0, px, py, pz, ee, mass ); + (fMasterGen->event).append( -particleID, 1, 0, 0, -px, -py, -pz, ee, mass ); } else { - (fMasterGen->event).append( particleID, 1, 0, 0, px, py, pz, ee, mass ); + (fMasterGen->event).append( particleID, 1, 0, 0, -px, -py, -pz, ee, mass ); } } @@ -92,6 +96,8 @@ bool Py8EGun::generatePartonsAndHadronize() if ( !fMasterGen->next() ) return false; + if (evtgenDecays) evtgenDecays->decay(); + event().reset(new HepMC::GenEvent); return toHepMC.fill_next_event( fMasterGen->event, event().get() ); diff --git a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc index 5fa91d5d5e526..a338632434b06 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc @@ -11,6 +11,8 @@ #include "Pythia8/Pythia.h" #include "Pythia8Plugins/HepMC2.h" +using namespace Pythia8; + #include "GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h" #include "GeneratorInterface/Pythia8Interface/plugins/ReweightUserHooks.h" @@ -26,6 +28,10 @@ #include "Pythia8Plugins/PowhegHooks.h" #include "GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h" +// EvtGen plugin +// +#include "Pythia8Plugins/EvtGen.h" + #include "FWCore/Concurrency/interface/SharedResourceNames.h" #include "FWCore/ServiceRegistry/interface/Service.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -50,7 +56,7 @@ namespace CLHEP { } using namespace gen; -using namespace Pythia8; + class Pythia8Hadronizer : public BaseHadronizer, public Py8InterfaceBase { @@ -81,7 +87,7 @@ class Pythia8Hadronizer : public BaseHadronizer, public Py8InterfaceBase { /// Center-of-Mass energy double comEnergy; - string LHEInputFileName; + std::string LHEInputFileName; std::auto_ptr lhaUP; enum { PP, PPbar, ElectronPositron }; @@ -136,7 +142,7 @@ const std::vector Pythia8Hadronizer::p8SharedResources = { edm::Sha Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet ¶ms) : BaseHadronizer(params), Py8InterfaceBase(params), comEnergy(params.getParameter("comEnergy")), - LHEInputFileName(params.getUntrackedParameter("LHEInputFileName","")), + LHEInputFileName(params.getUntrackedParameter("LHEInputFileName","")), fInitialState(PP), fReweightUserHook(0),fReweightRapUserHook(0),fReweightPtHatRapUserHook(0), fJetMatchingHook(0),fJetMatchingPy8InternalHook(0), fMergingHook(0), @@ -182,7 +188,7 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet ¶ms) : } if( params.exists( "SLHAFileForPythia8" ) ) { - std::string slhafilenameshort = params.getParameter("SLHAFileForPythia8"); + std::string slhafilenameshort = params.getParameter("SLHAFileForPythia8"); edm::FileInPath f1( slhafilenameshort ); fMasterGen->settings.mode("SLHA:readFrom", 2); @@ -197,7 +203,7 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet ¶ms) : } } else if( params.exists( "SLHATableForPythia8" ) ) { - std::string slhatable = params.getParameter("SLHATableForPythia8"); + std::string slhatable = params.getParameter("SLHATableForPythia8"); char tempslhaname[] = "pythia8SLHAtableXXXXXX"; int fd = mkstemp(tempslhaname); @@ -374,6 +380,13 @@ bool Pythia8Hadronizer::initializeForInternalPartons() edm::LogInfo("Pythia8Interface") << "Initializing Decayer"; status1 = fDecayer->init(); + if (useEvtGen) { + edm::LogInfo("Pythia8Interface") << "Creating and initializing pythia8 EvtGen plugin"; + + evtgenDecays = new EvtGenDecays(fMasterGen.get(), evtgenDecFile.c_str(), evtgenPdlFile.c_str()); + evtgenDecays->readDecayFile("evtgen_userfile.dec"); + } + return (status&&status1); } @@ -427,7 +440,7 @@ bool Pythia8Hadronizer::initializeForExternalPartons() } - if(LHEInputFileName != string()) { + if(LHEInputFileName != std::string()) { edm::LogInfo("Pythia8Interface") << "Initialize direct pythia8 reading from LHE file " << LHEInputFileName; @@ -467,6 +480,14 @@ bool Pythia8Hadronizer::initializeForExternalPartons() edm::LogInfo("Pythia8Interface") << "Initializing Decayer"; status1 = fDecayer->init(); + if (useEvtGen) { + edm::LogInfo("Pythia8Interface") << "Creating and initializing pythia8 EvtGen plugin"; + + std::string evtgenpath(getenv("EVTGENDATA")); + evtgenDecays = new EvtGenDecays(fMasterGen.get(), evtgenDecFile.c_str(), evtgenPdlFile.c_str()); + evtgenDecays->readDecayFile("evtgen_userfile.dec"); + } + return (status&&status1); } @@ -495,6 +516,8 @@ bool Pythia8Hadronizer::generatePartonsAndHadronize() if (!fMasterGen->next()) return false; + if (evtgenDecays) evtgenDecays->decay(); + event().reset(new HepMC::GenEvent); return toHepMC.fill_next_event( *(fMasterGen.get()), event().get()); @@ -506,7 +529,7 @@ bool Pythia8Hadronizer::hadronize() DJR.resize(0); nME = -1; nMEFiltered = -1; - if(LHEInputFileName == string()) lhaUP->loadEvent(lheEvent()); + if(LHEInputFileName == std::string()) lhaUP->loadEvent(lheEvent()); if ( fJetMatchingHook ) { @@ -546,6 +569,8 @@ bool Pythia8Hadronizer::hadronize() // lheEvent()->count( lhef::LHERunInfo::kAccepted, 1.0, mergeweight ); + if (evtgenDecays) evtgenDecays->decay(); + event().reset(new HepMC::GenEvent); bool py8hepmc = toHepMC.fill_next_event( *(fMasterGen.get()), event().get()); if (!py8hepmc) { @@ -563,8 +588,7 @@ bool Pythia8Hadronizer::hadronize() } return true; - - + } diff --git a/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc b/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc index e5018b5ea8d7b..17240678f817a 100644 --- a/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc +++ b/GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc @@ -1,6 +1,11 @@ #include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/Concurrency/interface/SharedResourceNames.h" +// EvtGen plugin +// +#include "Pythia8Plugins/EvtGen.h" + using namespace Pythia8; const std::vector gen::Py8GunBase::p8SharedResources = { edm::SharedResourceNames::kPythia8 }; @@ -11,27 +16,26 @@ Py8GunBase::Py8GunBase( edm::ParameterSet const& ps ) : Py8InterfaceBase(ps) { - runInfo().setFilterEfficiency( - ps.getUntrackedParameter("filterEfficiency", -1.) ); - runInfo().setExternalXSecLO( - GenRunInfoProduct::XSec(ps.getUntrackedParameter("crossSection", -1.)) ); - runInfo().setExternalXSecNLO( - GenRunInfoProduct::XSec(ps.getUntrackedParameter("crossSectionNLO", -1.)) ); + runInfo().setFilterEfficiency( + ps.getUntrackedParameter("filterEfficiency", -1.) ); + runInfo().setExternalXSecLO( + GenRunInfoProduct::XSec(ps.getUntrackedParameter("crossSection", -1.)) ); + runInfo().setExternalXSecNLO( + GenRunInfoProduct::XSec(ps.getUntrackedParameter("crossSectionNLO", -1.)) ); - // PGun specs // - edm::ParameterSet pgun_params = - ps.getParameter("PGunParameters"); + edm::ParameterSet pgun_params = + ps.getParameter("PGunParameters"); - // although there's the method ParameterSet::empty(), - // it looks like it's NOT even necessary to check if it is, - // before trying to extract parameters - if it is empty, - // the default values seem to be taken - // - fPartIDs = pgun_params.getParameter< std::vector >("ParticleID"); - fMinPhi = pgun_params.getParameter("MinPhi"); // ,-3.14159265358979323846); - fMaxPhi = pgun_params.getParameter("MaxPhi"); // , 3.14159265358979323846); + // although there's the method ParameterSet::empty(), + // it looks like it's NOT even necessary to check if it is, + // before trying to extract parameters - if it is empty, + // the default values seem to be taken + // + fPartIDs = pgun_params.getParameter< std::vector >("ParticleID"); + fMinPhi = pgun_params.getParameter("MinPhi"); // ,-3.14159265358979323846); + fMaxPhi = pgun_params.getParameter("MaxPhi"); // , 3.14159265358979323846); } @@ -40,20 +44,27 @@ Py8GunBase::Py8GunBase( edm::ParameterSet const& ps ) bool Py8GunBase::initializeForInternalPartons() { - // NO MATTER what was this setting below, override it before init - // - this is essencial for the PGun mode + // NO MATTER what was this setting below, override it before init + // - this is essencial for the PGun mode - // Key requirement: switch off ProcessLevel, and thereby also PartonLevel. - fMasterGen->readString("ProcessLevel:all = off"); - fMasterGen->readString("ProcessLevel::resonanceDecays=on"); - fMasterGen->init(); + // Key requirement: switch off ProcessLevel, and thereby also PartonLevel. + fMasterGen->readString("ProcessLevel:all = off"); + fMasterGen->readString("ProcessLevel::resonanceDecays=on"); + fMasterGen->init(); - // init decayer - fDecayer->readString("ProcessLevel:all = off"); // Same trick! - fDecayer->readString("ProcessLevel::resonanceDecays=on"); - fDecayer->init(); + // init decayer + fDecayer->readString("ProcessLevel:all = off"); // Same trick! + fDecayer->readString("ProcessLevel::resonanceDecays=on"); + fDecayer->init(); - return true; + if (useEvtGen) { + edm::LogInfo("Pythia8Interface") << "Creating and initializing pythia8 EvtGen plugin"; + + evtgenDecays = new EvtGenDecays(fMasterGen.get(), evtgenDecFile.c_str(), evtgenPdlFile.c_str()); + evtgenDecays->readDecayFile("evtgen_userfile.dec"); + } + + return true; } diff --git a/GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc b/GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc index c6c1eaeef1cfc..2c86a61a71f49 100644 --- a/GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc +++ b/GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc @@ -3,11 +3,16 @@ #include "FWCore/Utilities/interface/Exception.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" +// EvtGen plugin +// +//#include "Pythia8Plugins/EvtGen.h" + using namespace Pythia8; namespace gen { -Py8InterfaceBase::Py8InterfaceBase( edm::ParameterSet const& ps ) +Py8InterfaceBase::Py8InterfaceBase( edm::ParameterSet const& ps ) : +useEvtGen(false), evtgenDecays(0) { fMasterGen.reset(new Pythia); fDecayer.reset(new Pythia); @@ -27,6 +32,23 @@ Py8InterfaceBase::Py8InterfaceBase( edm::ParameterSet const& ps ) if(pythiaHepMCVerbosityParticles) ascii_io = new HepMC::IO_AsciiParticles("cout", std::ios::out); + + if ( ps.exists("useEvtGenPlugin") ) { + + useEvtGen = true; + + string evtgenpath(getenv("EVTGENDATA")); + evtgenDecFile = evtgenpath + string("/DECAY_2010.DEC"); + evtgenPdlFile = evtgenpath + string("/evt.pdl"); + + if ( ps.exists( "evtgenDecFile" ) ) + evtgenDecFile = ps.getParameter("evtgenDecFile"); + + if ( ps.exists( "evtgenPdlFile" ) ) + evtgenPdlFile = ps.getParameter("evtgenPdlFile"); + + } + } bool Py8InterfaceBase::readSettings( int )