From 72eefcdbceb0c73ed380f5807cb2f1baf928023e Mon Sep 17 00:00:00 2001 From: Josh Bendavid Date: Mon, 11 Jan 2016 01:56:34 +0100 Subject: [PATCH 1/3] add functionality for multiple pythia8 user hooks, plus new user hook for resonance decay filtering and example --- .../interface/MultiUserHook.h | 464 ++++++++++++++++++ .../interface/ResonanceDecayFilterHook.h | 29 ++ .../plugins/Pythia8Hadronizer.cc | 59 ++- .../Pythia8Interface/src/MultiUserHook.cc | 4 + .../src/ResonanceDecayFilterHook.cc | 107 ++++ .../test/resonanceDecayFilter.py | 61 +++ 6 files changed, 714 insertions(+), 10 deletions(-) create mode 100644 GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h create mode 100644 GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h create mode 100644 GeneratorInterface/Pythia8Interface/src/MultiUserHook.cc create mode 100644 GeneratorInterface/Pythia8Interface/src/ResonanceDecayFilterHook.cc create mode 100644 GeneratorInterface/Pythia8Interface/test/resonanceDecayFilter.py diff --git a/GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h b/GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h new file mode 100644 index 0000000000000..15afa28592566 --- /dev/null +++ b/GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h @@ -0,0 +1,464 @@ +//to allow combining multiple user hooks + +class MultiUserHook : public Pythia8::UserHooks { + +public: + + // Constructor and destructor. + MultiUserHook() {} + + unsigned int nHooks() const { return hooks_.size(); } + void addHook(Pythia8::UserHooks *hook) { hooks_.push_back(hook); } + +// Initialisation after beams have been set by Pythia::init(). + bool initAfterBeams() override { + bool test = true; + for (Pythia8::UserHooks *hook : hooks_) { + hook->initPtr(infoPtr, settingsPtr, + particleDataPtr, rndmPtr, + beamAPtr, beamBPtr, + beamPomAPtr, beamPomBPtr, + coupSMPtr, partonSystemsPtr, + sigmaTotPtr); + test &= hook->initAfterBeams(); + } + + //Certain types of hooks can only be handled by one UserHook at a time. + //Check for this and return an error if needed + int nCanVetoPT = 0; + int nCanVetoStep = 0; + int nCanVetoMPIStep = 0; + int nCanSetResonanceScale = 0; + int nCanEnhance = 0; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoPT()) ++nCanVetoPT; + if (hook->canVetoStep()) ++nCanVetoStep; + if (hook->canVetoMPIStep()) ++nCanVetoMPIStep; + if (hook->canSetResonanceScale()) ++nCanSetResonanceScale; + if (hook->canEnhanceEmission() || hook->canEnhanceTrial()) ++nCanEnhance; + } + + if (nCanVetoPT>1) { + infoPtr->errorMsg("Error in MultiUserHook::initAfterBeams " + "multiple UserHooks with canVetoPT() not allowed"); + test = false; + } + if (nCanVetoPT>1) { + infoPtr->errorMsg("Error in MultiUserHook::initAfterBeams " + "multiple UserHooks with canVetoStep() not allowed"); + test = false; + } + if (nCanVetoPT>1) { + infoPtr->errorMsg("Error in MultiUserHook::initAfterBeams " + "multiple UserHooks with canVetoMPIStep() not allowed"); + test = false; + } + if (nCanVetoPT>1) { + infoPtr->errorMsg("Error in MultiUserHook::initAfterBeams " + "multiple UserHooks with canSetResonanceScale() not allowed"); + test = false; + } + if (nCanEnhance>1) { + infoPtr->errorMsg("Error in MultiUserHook::initAfterBeams " + "multiple UserHooks with canEnhanceEmission() or canEnhanceTrial() not allowed"); + test = false; + } + + return test; + } + + // Possibility to modify cross section of process. + bool canModifySigma() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canModifySigma(); + } + return test; + } + + // Multiplicative factor modifying the cross section of a hard process. + double multiplySigmaBy(const Pythia8::SigmaProcess* sigmaProcessPtr, + const Pythia8::PhaseSpace* phaseSpacePtr, bool inEvent) override { + double factor = 1.; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canModifySigma()) factor *= hook->multiplySigmaBy(sigmaProcessPtr,phaseSpacePtr,inEvent); + } + return factor; + } + + // Possibility to bias selection of events, compensated by a weight. + bool canBiasSelection() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canBiasSelection(); + } + return test; + } + + // Multiplicative factor in the phase space selection of a hard process. + double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr, + const Pythia8::PhaseSpace* phaseSpacePtr, bool inEvent) override { + double factor = 1.; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canBiasSelection()) factor *= hook->biasSelectionBy(sigmaProcessPtr,phaseSpacePtr,inEvent); + } + return factor; + }; + + // Event weight to compensate for selection weight above. + double biasedSelectionWeight() override { + double factor = 1.; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canBiasSelection()) factor *= hook->biasedSelectionWeight(); + } + return factor; + }; + + // Possibility to veto event after process-level selection. + bool canVetoProcessLevel() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoProcessLevel(); + } + return test; + } + + // Decide whether to veto current process or not, based on process record. + // Usage: doVetoProcessLevel( process). + bool doVetoProcessLevel(Pythia8::Event& event) override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoProcessLevel()) test |= hook->doVetoProcessLevel(event); + } + return test; + } + + // Possibility to veto resonance decay chain. + bool canVetoResonanceDecays() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoResonanceDecays(); + } + return test; + } + + // Decide whether to veto current resonance decay chain or not, based on + // process record. Usage: doVetoProcessLevel( process). + bool doVetoResonanceDecays(Pythia8::Event& event) override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoResonanceDecays()) test |= hook->doVetoResonanceDecays(event); + } + return test; + } + + // Possibility to veto MPI + ISR + FSR evolution and kill event, + // making decision at a fixed pT scale. Useful for MLM-style matching. + bool canVetoPT() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoPT(); + } + return test; + } + + // Transverse-momentum scale for veto test. (only one hook allowed) + double scaleVetoPT() override { + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoPT()) return hook->scaleVetoPT(); + } + return 0.; + }; + + // Decide whether to veto current event or not, based on event record. + // Usage: doVetoPT( iPos, event), where iPos = 0: no emissions so far; + // iPos = 1/2/3 joint evolution, latest step was MPI/ISR/FSR; + // iPos = 4: FSR only afterwards; iPos = 5: FSR in resonance decay. + bool doVetoPT( int iPos, const Pythia8::Event& event) override { + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoPT()) return hook->doVetoPT(iPos,event); + } + return false; + } + + // Possibility to veto MPI + ISR + FSR evolution and kill event, + // making decision after fixed number of ISR or FSR steps. + bool canVetoStep() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoStep(); + } + return test; + } + + // Up to how many ISR + FSR steps of hardest interaction should be checked. + int numberVetoStep() override { + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoStep()) return hook->numberVetoStep(); + } + return 1; + }; + + // Decide whether to veto current event or not, based on event record. + // Usage: doVetoStep( iPos, nISR, nFSR, event), where iPos as above, + // nISR and nFSR number of emissions so far for hard interaction only. + bool doVetoStep( int iPos, int nISR, int nFSR, const Pythia8::Event& event) override { + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoStep()) return hook->doVetoStep(iPos,nISR,nFSR,event); + } + return false; + } + + // Possibility to veto MPI + ISR + FSR evolution and kill event, + // making decision after fixed number of MPI steps. + bool canVetoMPIStep() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoMPIStep(); + } + return test; + } + + // Up to how many MPI steps should be checked. + int numberVetoMPIStep() override { + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoMPIStep()) return hook->numberVetoMPIStep(); + } + return 1; + }; + + // Decide whether to veto current event or not, based on event record. + // Usage: doVetoMPIStep( nMPI, event), where nMPI is number of MPI's so far. + bool doVetoMPIStep( int nMPI, const Pythia8::Event& event) override { + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoMPIStep()) return hook->doVetoMPIStep(nMPI,event); + } + return false; + } + + // Possibility to veto event after ISR + FSR + MPI in parton level, + // but before beam remnants and resonance decays. + bool canVetoPartonLevelEarly() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoPartonLevelEarly(); + } + return test; + } + + // Decide whether to veto current partons or not, based on event record. + // Usage: doVetoPartonLevelEarly( event). + bool doVetoPartonLevelEarly( const Pythia8::Event& event) override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoPartonLevelEarly()) test |= hook->doVetoPartonLevelEarly(event); + } + return test; + } + + // Retry same ProcessLevel with a new PartonLevel after a veto in + // doVetoPT, doVetoStep, doVetoMPIStep or doVetoPartonLevelEarly + // if you overload this method to return true. + bool retryPartonLevel() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->retryPartonLevel(); + } + return test; + } + + // Possibility to veto event after parton-level selection. + bool canVetoPartonLevel() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoPartonLevel(); + } + return test; + } + + // Decide whether to veto current partons or not, based on event record. + // Usage: doVetoPartonLevel( event). + bool doVetoPartonLevel( const Pythia8::Event& event) override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoPartonLevel()) test |= hook->doVetoPartonLevel(event); + } + return test; + } + + // Possibility to set initial scale in TimeShower for resonance decay. + bool canSetResonanceScale() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canSetResonanceScale(); + } + return test; + } + + // Initial scale for TimeShower evolution. + // Usage: scaleResonance( iRes, event), where iRes is location + // of decaying resonance in the event record. + // Only one UserHook setting the resonance scale is allowed + double scaleResonance( int iRes, const Pythia8::Event& event) override { + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canSetResonanceScale()) return hook->scaleResonance(iRes,event); + } + return 0.; + }; + + // Possibility to veto an emission in the ISR machinery. + bool canVetoISREmission() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoISREmission(); + } + return test; + } + + // Decide whether to veto current emission or not, based on event record. + // Usage: doVetoISREmission( sizeOld, event, iSys) where sizeOld is size + // of event record before current emission-to-be-scrutinized was added, + // and iSys is the system of the radiation (according to PartonSystems). + bool doVetoISREmission( int sizeOld, const Pythia8::Event& event, int iSys) override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoISREmission()) test |= hook->doVetoISREmission(sizeOld,event,iSys); + } + return test; + } + + // Possibility to veto an emission in the FSR machinery. + bool canVetoFSREmission() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoFSREmission(); + } + return test; + } + + // Decide whether to veto current emission or not, based on event record. + // Usage: doVetoFSREmission( sizeOld, event, iSys, inResonance) where + // sizeOld is size of event record before current emission-to-be-scrutinized + // was added, iSys is the system of the radiation (according to + // PartonSystems), and inResonance is true if the emission takes place in a + // resonance decay. + bool doVetoFSREmission( int sizeOld, const Pythia8::Event& event, int iSys, bool inResonance = false ) override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoFSREmission()) test |= hook->doVetoFSREmission(sizeOld,event,iSys,inResonance); + } + return test; + } + + // Possibility to veto an MPI. + bool canVetoMPIEmission() { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canVetoMPIEmission(); + } + return test; + } + + // Decide whether to veto an MPI based on event record. + // Usage: doVetoMPIEmission( sizeOld, event) where sizeOld + // is size of event record before the current MPI. + bool doVetoMPIEmission( int sizeOld, const Pythia8::Event &event) override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canVetoMPIEmission()) test |= hook->doVetoMPIEmission(sizeOld,event); + } + return test; + } + + // Possibility to reconnect colours from resonance decay systems. + bool canReconnectResonanceSystems() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canReconnectResonanceSystems(); + } + return test; + } + + // Do reconnect colours from resonance decay systems. + // Usage: doVetoFSREmission( oldSizeEvt, event) + // where oldSizeEvent is the event size before resonance decays. + // Should normally return true, while false means serious failure. + // Value of PartonLevel:earlyResDec determines where method is called. + bool doReconnectResonanceSystems( int oldSizeEvt, Pythia8::Event &event) override { + bool test = true; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canReconnectResonanceSystems()) test &= hook->doReconnectResonanceSystems(oldSizeEvt,event); + } + return test; + } + + // Enhance emission rates (sec. 4 in EPJC (2013) 73). + bool canEnhanceEmission() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canEnhanceEmission(); + } + return test; + } + + // Bookkeeping of weights for enhanced actual or trial emissions + // (sec. 3 in EPJC (2013) 73). + bool canEnhanceTrial() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canEnhanceTrial(); + } + return test; + } + + double enhanceFactor( std::string str) override { + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canEnhanceEmission() || hook->canEnhanceTrial()) return hook->enhanceFactor(str); + } + return 1.; + }; + + double vetoProbability( std::string str ) override { + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canEnhanceEmission() || hook->canEnhanceTrial()) return hook->vetoProbability(str); + } + return 0.; + }; + + // Can change fragmentation parameters. + bool canChangeFragPar() override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + test |= hook->canChangeFragPar(); + } + return test; + } + + // Do change fragmentation parameters. + // Input: flavPtr, zPtr, pTPtr, idEnd, m2Had, iParton. + bool doChangeFragPar( Pythia8::StringFlav* flavPtr, Pythia8::StringZ* zPtr, Pythia8::StringPT* pTPtr, int idEnd, + double m2Had, std::vector iParton) override { + bool test = true; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canChangeFragPar()) test &= hook->doChangeFragPar(flavPtr, zPtr, pTPtr, idEnd, m2Had, iParton); + } + return test; + } + + // Do a veto on a hadron just before it is added to the final state. + bool doVetoFragmentation( Pythia8::Particle part) override { + bool test = false; + for (Pythia8::UserHooks *hook : hooks_) { + if (hook->canChangeFragPar()) test |= hook->doVetoFragmentation(part); + } + return test; + } + + + +//-------------------------------------------------------------------------- + +private: + std::vector hooks_; + +}; diff --git a/GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h b/GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h new file mode 100644 index 0000000000000..62b8f0b162461 --- /dev/null +++ b/GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h @@ -0,0 +1,29 @@ +class ResonanceDecayFilterHook : public Pythia8::UserHooks { + +public: + + // Constructor and destructor. + ResonanceDecayFilterHook() {} + +//-------------------------------------------------------------------------- + + bool initAfterBeams() override; + bool canVetoResonanceDecays() override { return true; } + bool doVetoResonanceDecays(Pythia8::Event& process) override { return checkVetoResonanceDecays(process); } + bool checkVetoResonanceDecays(const Pythia8::Event& process); + +//-------------------------------------------------------------------------- + +private: + bool filter_; + bool exclusive_; + bool eMuAsEquivalent_; + bool eMuTauAsEquivalent_; + bool allNuAsEquivalent_; + std::vector mothers_; + std::vector daughters_; + + std::map requestedDaughters_; + std::map observedDaughters_; + +}; diff --git a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc index a338632434b06..efd1b55c4b07f 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc @@ -23,11 +23,16 @@ using namespace Pythia8; #include "Pythia8Plugins/JetMatching.h" #include "Pythia8Plugins/aMCatNLOHooks.h" +#include "GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h" + // Emission Veto Hooks // #include "Pythia8Plugins/PowhegHooks.h" #include "GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h" +//decay filter hook +#include "GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h" + // EvtGen plugin // #include "Pythia8Plugins/EvtGen.h" @@ -96,6 +101,9 @@ class Pythia8Hadronizer : public BaseHadronizer, public Py8InterfaceBase { double fBeam1PZ; double fBeam2PZ; + //helper class to allow multiple user hooks simultaneously + MultiUserHook *fMultiUserHook; + // Reweight user hooks // UserHooks* fReweightUserHook; @@ -113,6 +121,9 @@ class Pythia8Hadronizer : public BaseHadronizer, public Py8InterfaceBase { PowhegHooks* fEmissionVetoHook; EmissionVetoHook1* fEmissionVetoHook1; + //resonance decay filter hook + ResonanceDecayFilterHook *fResonanceDecayFilterHook; + int EV1_nFinal; bool EV1_vetoOn; int EV1_maxVetoCount; @@ -144,9 +155,9 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet ¶ms) : comEnergy(params.getParameter("comEnergy")), LHEInputFileName(params.getUntrackedParameter("LHEInputFileName","")), fInitialState(PP), - fReweightUserHook(0),fReweightRapUserHook(0),fReweightPtHatRapUserHook(0), + fMultiUserHook(new MultiUserHook), fReweightUserHook(0),fReweightRapUserHook(0),fReweightPtHatRapUserHook(0), fJetMatchingHook(0),fJetMatchingPy8InternalHook(0), fMergingHook(0), - fEmissionVetoHook(0), fEmissionVetoHook1(0), nME(-1), nMEFiltered(-1), nISRveto(0), nFSRveto(0), + fEmissionVetoHook(0), fEmissionVetoHook1(0), fResonanceDecayFilterHook(0), nME(-1), nMEFiltered(-1), nISRveto(0), nFSRveto(0), NHooks(0) { @@ -223,6 +234,15 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet ¶ms) : << std::endl; } } + + //add settings for resonance decay filter + fMasterGen->settings.addFlag("ResonanceDecayFilter:filter",false); + fMasterGen->settings.addFlag("ResonanceDecayFilter:exclusive",false); + fMasterGen->settings.addFlag("ResonanceDecayFilter:eMuAsEquivalent",false); + fMasterGen->settings.addFlag("ResonanceDecayFilter:eMuTauAsEquivalent",false); + fMasterGen->settings.addFlag("ResonanceDecayFilter:allNuAsEquivalent",false); + fMasterGen->settings.addMVec("ResonanceDecayFilter:mothers",std::vector(),false,false,0,0); + fMasterGen->settings.addMVec("ResonanceDecayFilter:daughters",std::vector(),false,false,0,0); // Reweight user hook // @@ -310,13 +330,13 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet ¶ms) : if(NHooks > 1) throw edm::Exception(edm::errors::Configuration,"Pythia8Interface") <<" Too many User Hooks. \n Please choose one from: reweightGen, reweightGenRap, reweightGenPtHatRap, jetMatching, emissionVeto1 \n"; - if(fReweightUserHook) fMasterGen->setUserHooksPtr(fReweightUserHook); - if(fReweightRapUserHook) fMasterGen->setUserHooksPtr(fReweightRapUserHook); - if(fReweightPtHatRapUserHook) fMasterGen->setUserHooksPtr(fReweightPtHatRapUserHook); - if(fJetMatchingHook) fMasterGen->setUserHooksPtr(fJetMatchingHook); + if(fReweightUserHook) fMultiUserHook->addHook(fReweightUserHook); + if(fReweightRapUserHook) fMultiUserHook->addHook(fReweightRapUserHook); + if(fReweightPtHatRapUserHook) fMultiUserHook->addHook(fReweightPtHatRapUserHook); + if(fJetMatchingHook) fMultiUserHook->addHook(fJetMatchingHook); if(fEmissionVetoHook1) { edm::LogInfo("Pythia8Interface") << "Turning on Emission Veto Hook 1 from CMSSW Pythia8Interface"; - fMasterGen->setUserHooksPtr(fEmissionVetoHook1); + fMultiUserHook->addHook(fEmissionVetoHook1); } } @@ -361,6 +381,16 @@ bool Pythia8Hadronizer::initializeForInternalPartons() throw edm::Exception(edm::errors::Configuration,"Pythia8Interface") <<" UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron \n"; } + + bool resonanceDecayFilter = fMasterGen->settings.flag("ResonanceDecayFilter:filter"); + if (resonanceDecayFilter && !fResonanceDecayFilterHook) { + fResonanceDecayFilterHook = new ResonanceDecayFilterHook; + fMultiUserHook->addHook(fResonanceDecayFilterHook); + } + + if (fMultiUserHook->nHooks()>0) { + fMasterGen->setUserHooksPtr(fMultiUserHook); + } fMasterGen->settings.parm("Beams:eCM", comEnergy); edm::LogInfo("Pythia8Interface") << "Initializing MasterGen"; @@ -407,7 +437,7 @@ bool Pythia8Hadronizer::initializeForExternalPartons() fEmissionVetoHook = new PowhegHooks(); edm::LogInfo("Pythia8Interface") << "Turning on Emission Veto Hook from pythia8 code"; - fMasterGen->setUserHooksPtr(fEmissionVetoHook); + fMultiUserHook->addHook(fEmissionVetoHook); } @@ -422,7 +452,7 @@ bool Pythia8Hadronizer::initializeForExternalPartons() if (internalMatching && !fJetMatchingPy8InternalHook) { fJetMatchingPy8InternalHook = new Pythia8::JetMatchingMadgraph; - fMasterGen->setUserHooksPtr(fJetMatchingPy8InternalHook); + fMultiUserHook->addHook(fJetMatchingPy8InternalHook); } if (internalMerging && !fMergingHook) { @@ -436,9 +466,18 @@ bool Pythia8Hadronizer::initializeForExternalPartons() 2 : 0 ); fMergingHook = new Pythia8::amcnlo_unitarised_interface(scheme); - fMasterGen->setUserHooksPtr(fMergingHook); + fMultiUserHook->addHook(fMergingHook); + } + + bool resonanceDecayFilter = fMasterGen->settings.flag("ResonanceDecayFilter:filter"); + if (resonanceDecayFilter && !fResonanceDecayFilterHook) { + fResonanceDecayFilterHook = new ResonanceDecayFilterHook; + fMultiUserHook->addHook(fResonanceDecayFilterHook); } + if (fMultiUserHook->nHooks()>0) { + fMasterGen->setUserHooksPtr(fMultiUserHook); + } if(LHEInputFileName != std::string()) { diff --git a/GeneratorInterface/Pythia8Interface/src/MultiUserHook.cc b/GeneratorInterface/Pythia8Interface/src/MultiUserHook.cc new file mode 100644 index 0000000000000..0457768a834e4 --- /dev/null +++ b/GeneratorInterface/Pythia8Interface/src/MultiUserHook.cc @@ -0,0 +1,4 @@ +#include "Pythia8/Pythia.h" +#include "GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h" + +using namespace Pythia8; diff --git a/GeneratorInterface/Pythia8Interface/src/ResonanceDecayFilterHook.cc b/GeneratorInterface/Pythia8Interface/src/ResonanceDecayFilterHook.cc new file mode 100644 index 0000000000000..e014bb61cfc84 --- /dev/null +++ b/GeneratorInterface/Pythia8Interface/src/ResonanceDecayFilterHook.cc @@ -0,0 +1,107 @@ +#include "Pythia8/Pythia.h" +#include "GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h" + +using namespace Pythia8; + + +//-------------------------------------------------------------------------- +bool ResonanceDecayFilterHook::initAfterBeams() { + filter_ = settingsPtr->flag("ResonanceDecayFilter:filter"); + exclusive_ = settingsPtr->flag("ResonanceDecayFilter:exclusive"); + eMuAsEquivalent_ = settingsPtr->flag("ResonanceDecayFilter:eMuAsEquivalent"); + eMuTauAsEquivalent_ = settingsPtr->flag("ResonanceDecayFilter:eMuTauAsEquivalent"); + allNuAsEquivalent_ = settingsPtr->flag("ResonanceDecayFilter:allNuAsEquivalent_"); + mothers_ = settingsPtr->mvec("ResonanceDecayFilter:mothers"); + daughters_ = settingsPtr->mvec("ResonanceDecayFilter:daughters"); + + + requestedDaughters_.clear(); + + for (int id : daughters_) { + int did = std::abs(id); + if ( did == 13 && (eMuAsEquivalent_ || eMuTauAsEquivalent_)) { + did = 11; + } + if ( did == 15 && eMuTauAsEquivalent_) { + did = 11; + } + if ( (did == 14 || did == 16) && allNuAsEquivalent_) { + did = 12; + } + + ++requestedDaughters_[std::abs(did)]; + } + + return true; + +} + +//-------------------------------------------------------------------------- +bool ResonanceDecayFilterHook::checkVetoResonanceDecays(const Event& process) { + + if (!filter_) return false; + + observedDaughters_.clear(); + + //count decay products + for (int i=0; i0 ? std::abs(process[p.mother1()].id()) : 0; + + //if no list of mothers is provided, then all particles + //in hard process and resonance decays are counted together + bool validMother = mothers_.size() ? false : true; + for (int id : mothers_) { + if (mid == std::abs(id)) { + validMother = true; + break; + } + } + + if (validMother) { + ++observedDaughters_[did]; + } + } + + //check if criteria is satisfied + //inclusive mode: at least as many decay products as requested + //exclusive mode: exactly as many decay products as requested + //(but additional particle types not appearing in the list of requested daughter id's are ignored) + for (const auto &reqpair : requestedDaughters_) { + int reqid = reqpair.first; + int reqcount = reqpair.second; + + int obscount = 0; + for (const auto &obspair : observedDaughters_) { + int obsid = obspair.first; + + if (obsid == reqid) { + obscount = obspair.second; + break; + } + } + + //inclusive criteria not satisfied, veto event + if (obscount < reqcount) return true; + + //exclusive criteria not satisfied, veto event + if (exclusive_ && obscount > reqcount) return true; + + } + + //all criteria satisfied, don't veto + return false; +} \ No newline at end of file diff --git a/GeneratorInterface/Pythia8Interface/test/resonanceDecayFilter.py b/GeneratorInterface/Pythia8Interface/test/resonanceDecayFilter.py new file mode 100644 index 0000000000000..80a216046be6e --- /dev/null +++ b/GeneratorInterface/Pythia8Interface/test/resonanceDecayFilter.py @@ -0,0 +1,61 @@ +#example gen fragment that takes a gridapack which produces +#ttH events with inclusive top decays and undecayed higgs at lhe level +#and selects resonance decays such that events have at least four leptons (electrons, muons, taus) + +import FWCore.ParameterSet.Config as cms + +# link to cards: +# https://github.com/cms-sw/genproductions/blob/0d4b4288fa053d9a8aef5c6e123b66bf94c3aee8/bin/Powheg/production/V2/13TeV/Higgs/ttH_inclusive_NNPDF30_13TeV_M125/ttH_inclusive_NNPDF30_13TeV_M125.input + +externalLHEProducer = cms.EDProducer("ExternalLHEProducer", + args = cms.vstring('/cvmfs/cms.cern.ch/phys_generator/gridpacks/slc6_amd64_gcc481/13TeV/powheg/V2/ttH_inclusive_NNPDF30_13TeV_M125/v1/ttH_inclusive_NNPDF30_13TeV_M125_tarball.tar.gz'), + nEvents = cms.untracked.uint32(5000), + numberOfParameters = cms.uint32(1), + outputFile = cms.string('cmsgrid_final.lhe'), + scriptName = cms.FileInPath('GeneratorInterface/LHEInterface/data/run_generic_tarball_cvmfs.sh') +) + +from Configuration.Generator.Pythia8CommonSettings_cfi import * +from Configuration.Generator.Pythia8CUEP8M1Settings_cfi import * +from Configuration.Generator.Pythia8PowhegEmissionVetoSettings_cfi import * + +generator = cms.EDFilter("Pythia8HadronizerFilter", + maxEventsToPrint = cms.untracked.int32(1), + pythiaPylistVerbosity = cms.untracked.int32(1), + filterEfficiency = cms.untracked.double(1.0), + pythiaHepMCVerbosity = cms.untracked.bool(False), + comEnergy = cms.double(13000.), + PythiaParameters = cms.PSet( + pythia8CommonSettingsBlock, + pythia8CUEP8M1SettingsBlock, + pythia8PowhegEmissionVetoSettingsBlock, + processParameters = cms.vstring( + 'POWHEG:nFinal = 3', ## Number of final state particles + ## (BEFORE THE DECAYS) in the LHE + ## other than emitted extra parton + '25:m0 = 125.0', + '25:onMode = off', + '25:onIfMatch = 15 -15', + '25:onIfMatch = 23 23', + '25:onIfMatch = 24 -24', + '24:mMin = 0.05', + '23:mMin = 0.05', + 'ResonanceDecayFilter:filter = on', + 'ResonanceDecayFilter:exclusive = off', #off: require at least the specified number of daughters, on: require exactly the specified number of daughters + 'ResonanceDecayFilter:eMuAsEquivalent = off', #on: treat electrons and muons as equivalent + 'ResonanceDecayFilter:eMuTauAsEquivalent = on', #on: treat electrons, muons , and taus as equivalent + 'ResonanceDecayFilter:allNuAsEquivalent = on', #on: treat all three neutrino flavours as equivalent + #'ResonanceDecayFilter:mothers =', #list of mothers not specified -> count all particles in hard process+resonance decays (better to avoid specifying mothers when including leptons from the lhe in counting, since intermediate resonances are not gauranteed to appear in general + 'ResonanceDecayFilter:daughters = 11,11,11,11', + ), + parameterSets = cms.vstring('pythia8CommonSettings', + 'pythia8CUEP8M1Settings', + 'pythia8PowhegEmissionVetoSettings', + 'processParameters' + ) + ) + ) + +ProductionFilterSequence = cms.Sequence(generator) + + From 3cf16ce195651d4d7f82f6e970d69cd575b27a5e Mon Sep 17 00:00:00 2001 From: Josh Bendavid Date: Thu, 21 Jan 2016 17:18:21 +0100 Subject: [PATCH 2/3] fix mistake in protection for incompatible hooks --- .../Pythia8Interface/interface/MultiUserHook.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h b/GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h index 15afa28592566..198c90d58325c 100644 --- a/GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h +++ b/GeneratorInterface/Pythia8Interface/interface/MultiUserHook.h @@ -43,17 +43,17 @@ class MultiUserHook : public Pythia8::UserHooks { "multiple UserHooks with canVetoPT() not allowed"); test = false; } - if (nCanVetoPT>1) { + if (nCanVetoStep>1) { infoPtr->errorMsg("Error in MultiUserHook::initAfterBeams " "multiple UserHooks with canVetoStep() not allowed"); test = false; } - if (nCanVetoPT>1) { + if (nCanVetoMPIStep>1) { infoPtr->errorMsg("Error in MultiUserHook::initAfterBeams " "multiple UserHooks with canVetoMPIStep() not allowed"); test = false; } - if (nCanVetoPT>1) { + if (nCanSetResonanceScale>1) { infoPtr->errorMsg("Error in MultiUserHook::initAfterBeams " "multiple UserHooks with canSetResonanceScale() not allowed"); test = false; From d680abf1aae0cc71210d4503a2d63db8b0cbf0f6 Mon Sep 17 00:00:00 2001 From: Josh Bendavid Date: Sat, 23 Jan 2016 19:03:33 +0100 Subject: [PATCH 3/3] delete multiuserhook pointer --- .../Pythia8Interface/plugins/Pythia8Hadronizer.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc index efd1b55c4b07f..84ce6c90e98fc 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc @@ -347,6 +347,10 @@ Pythia8Hadronizer::~Pythia8Hadronizer() // do we need to delete UserHooks/JetMatchingHook here ??? if(fEmissionVetoHook) {delete fEmissionVetoHook; fEmissionVetoHook=0;} if(fEmissionVetoHook1) {delete fEmissionVetoHook1; fEmissionVetoHook1=0;} + if (fMultiUserHook) { + delete fMultiUserHook; + fMultiUserHook = 0; + } //clean up temp file if (!slhafile_.empty()) {