include/Pythia8Plugins/ResonanceDecayFilterHook.h

0001 // ResonanceDecayFilterHook.h is part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Stephen Mrenna, Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 // Author: Stephen Mrenna, January 2021.
0006
0007 // This class is used to filter out resonance decay products to obtain
0008 // a desired parton-level final state.  It copies an implementation
0009 // already in the CMS software originally authored by Josh Bendavid,
0010 // with some minor modifications.
0011
0012 // The following settings are available:
0013 // flag ResonanceDecayFilter:filter = false
0014 //   Activate filter if true.
0015 // flag ResonanceDecayFilter:exclusive = false
0016 //   Demand exactly the particles requested by the filter if true.
0017 //   If false, additional particles can be present in the resonance decays.
0018 // flag ResonanceDecayFilter:eMuAsEquivalent = false
0019 //   If true, treat the particle id for electron and muon as the same.
0020 // flag ResonanceDecayFilter:eMuTauAsEquivalent = false
0021 //  As above, but include tau leptons.
0022 // flag ResonanceDecayFilter:allNuAsEquivalent = false
0023 // flag ResonanceDecayFilter:udscAsEquivalent = false
0024 //   If true, treat the particle id for all up, down, strange, and charm
0025 //   quarks as the same.</flag>
0026 // flag ResonanceDecayFilter:udscbAsEquivalent = false
0027 //   As above, but include bottom quarks.
0028 // flag ResonanceDecayFilter:wzAsEquivalent = false
0029 //   If true, treat the particle id for W and Z bosons as the same.
0030 // mvec ResonanceDecayFilter:mothers = {}
0031 //   If provided, the set of mothers that will be filtered.
0032 //   Otherwise, all resonance decay products are considered.
0033 // mvec ResonanceDecayFilter:daughters = {}
0034 //   A list of particle ids for resonance decay products that
0035 //   will pass the filter.   If empty, all particles will pass the
0036 //   filter. Beware: an infinite loop is possible if the daughter never
0037 //   appears in resonance decays.
0038
0039 #ifndef Pythia8_ResonanceDecayFilterHooks_H
0040 #define Pythia8_ResonanceDecayFilterHooks_H
0041
0042 // Includes.
0043 #include "Pythia8/Pythia.h"
0044 #include "Pythia8/UserHooks.h"
0045 #include "Pythia8/Event.h"
0046
0047 namespace Pythia8 {
0048
0049 class ResonanceDecayFilterHook : public UserHooks {
0050 public:
0051
0052   // Constructor.
0053   ResonanceDecayFilterHook(Settings &settings);
0054
0055   // Override base class methods.
0056   bool canVetoResonanceDecays() override {return true;}
0057   bool doVetoResonanceDecays(Event& process) override {
0058     return checkVetoResonanceDecays(process);}
0059   bool initAfterBeams() override;
0060
0061   // Class specific.
0062   bool checkVetoResonanceDecays(const Event& process);
0063   unsigned long int returnCounter() {return counter;};
0064
0065 private:
0066
0067   // Return the particle ID category.
0068   int idCat(int id);
0069
0070   // Data members.
0071   bool filter, exclusive, eMuAsEquivalent, eMuTauAsEquivalent,
0072     allNuAsEquivalent, udscAsEquivalent, udscbAsEquivalent, wzAsEquivalent;
0073   unsigned long int counter;
0074   set<int> mothers;
0075   vector<int> daughters;
0076   unordered_map<int, int> requestedDaughters, observedDaughters;
0077
0078 };
0079
0080 //--------------------------------------------------------------------------
0081
0082 // Constructor.
0083
0084 ResonanceDecayFilterHook::ResonanceDecayFilterHook(Settings &settings) {
0085   counter = 0;
0086   settings.addFlag("ResonanceDecayFilter:filter", false);
0087   settings.addFlag("ResonanceDecayFilter:exclusive", false);
0088   settings.addFlag("ResonanceDecayFilter:eMuAsEquivalent", false);
0089   settings.addFlag("ResonanceDecayFilter:eMuTauAsEquivalent", false);
0090   settings.addFlag("ResonanceDecayFilter:allNuAsEquivalent", false);
0091   settings.addFlag("ResonanceDecayFilter:udscAsEquivalent", false);
0092   settings.addFlag("ResonanceDecayFilter:udscbAsEquivalent", false);
0093   settings.addFlag("ResonanceDecayFilter:wzAsEquivalent", false);
0094   settings.addMVec("ResonanceDecayFilter:mothers", vector<int>(), false,
0095     false, 0, 0);
0096   settings.addMVec("ResonanceDecayFilter:daughters", vector<int>(), false,
0097     false, 0, 0);
0098 }
0099
0100 //--------------------------------------------------------------------------
0101
0102 // Return a particle ID given equivalence user settings.
0103
0104 int ResonanceDecayFilterHook::idCat(int id) {
0105   id = abs(id);
0106   if (id == 13 && (eMuAsEquivalent || eMuTauAsEquivalent)) id = 11;
0107   else if (id == 15 && eMuTauAsEquivalent) id = 11;
0108   else if ((id == 14 || id == 16) && allNuAsEquivalent) id = 12;
0109   else if ((id == 2 || id == 3 || id == 4) && udscAsEquivalent) id = 1;
0110   else if ((id == 2 || id == 3 || id == 4 || id == 5) &&
0111            udscbAsEquivalent) id = 1;
0112   else if ((id == 23 || id == 24) && wzAsEquivalent) id = 23;
0113   return id;
0114 }
0115
0116 //--------------------------------------------------------------------------
0117
0118 // Intialize the user hook after the beams.
0119
0120 bool ResonanceDecayFilterHook::initAfterBeams() {
0121   filter = settingsPtr->flag("ResonanceDecayFilter:filter");
0122   exclusive = settingsPtr->flag("ResonanceDecayFilter:exclusive");
0123   eMuAsEquivalent = settingsPtr->flag("ResonanceDecayFilter:eMuAsEquivalent");
0124   eMuTauAsEquivalent = settingsPtr->
0125     flag("ResonanceDecayFilter:eMuTauAsEquivalent");
0126   allNuAsEquivalent = settingsPtr->flag
0127     ("ResonanceDecayFilter:allNuAsEquivalent");
0128   udscAsEquivalent = settingsPtr->
0129     flag("ResonanceDecayFilter:udscAsEquivalent");
0130   udscbAsEquivalent = settingsPtr->
0131     flag("ResonanceDecayFilter:udscbAsEquivalent");
0132   wzAsEquivalent = settingsPtr->flag("ResonanceDecayFilter:wzAsEquivalent");
0133   auto mothersIn = settingsPtr->mvec("ResonanceDecayFilter:mothers");
0134   mothers.clear();
0135   mothers.insert(mothersIn.begin(), mothersIn.end());
0136   daughters = settingsPtr->mvec("ResonanceDecayFilter:daughters");
0137   requestedDaughters.clear();
0138
0139   // Loop over the daughters.
0140   for (int id : daughters) ++requestedDaughters[idCat(id)];
0141   return true;
0142
0143 }
0144
0145 //--------------------------------------------------------------------------
0146
0147 // Return true of the resonance decays are vetoed.
0148
0149 bool ResonanceDecayFilterHook::checkVetoResonanceDecays(const Event &process) {
0150   if (!filter) return false;
0151
0152   // Count the number of times hook is called.
0153   counter++;
0154   observedDaughters.clear();
0155
0156   // Loop over particles and determine equivalent types.
0157   for (int i = 0; i < process.size(); ++i) {
0158     const Particle &p = process[i];
0159     int mid = p.mother1() > 0 ? abs(process[p.mother1()].id()) : 0;
0160
0161     // If no list of mothers is provided, then all particles
0162     // in hard process and resonance decays are counted together
0163     if (mothers.empty() || mothers.count(mid) || mothers.count(-mid))
0164       ++observedDaughters[idCat(p.id())];
0165   }
0166
0167   // Check if criteria is satisfied.
0168   // inclusive mode: at least as many decay products as requested
0169   // exclusive mode: exactly as many decay products as requested
0170   // (but additional particle types not appearing in the list of requested
0171   // daughter id's are ignored)
0172   for (const auto &reqDau : requestedDaughters) {
0173     int reqId = reqDau.first;
0174     int reqCount = reqDau.second;
0175     auto obsItr = observedDaughters.find(reqId);
0176     int obsCount = ( obsItr != observedDaughters.end() ) ? obsItr->second : 0;
0177
0178     // Inclusive criteria not satisfied, veto event
0179     if (obsCount < reqCount) return true;
0180
0181     // Exclusive criteria not satisfied, veto event
0182     if (exclusive && obsCount > reqCount) return true;
0183   }
0184
0185   // All criteria satisfied, don't veto
0186   return false;
0187
0188 }
0189
0190 //==========================================================================
0191
0192 } // end namespace Pythia8
0193
0194 #endif // end Pythia8_ResonanceDecayHooks_H