Rivet/Tools/JetUtils.hh

0001 #ifndef RIVET_JETUTILS_HH
0002 #define RIVET_JETUTILS_HH
0003
0004 #include "Rivet/Jet.hh"
0005 #include "Rivet/Tools/ParticleBaseUtils.hh"
0006
0007 namespace Rivet {
0008
0009
0010   /// @defgroup jetutils Functions for Jets
0011   /// @{
0012
0013   /// @defgroup jetutils_conv Converting between Jets, Particles and PseudoJets
0014   /// @{
0015
0016   inline PseudoJets mkPseudoJets(const Particles& ps) {
0017     PseudoJets rtn; rtn.reserve(ps.size());
0018     for (const Particle& p : ps) rtn.push_back(p.pseudojet());
0019     return rtn;
0020   }
0021
0022   inline PseudoJets mkPseudoJets(const Jets& js) {
0023     PseudoJets rtn; rtn.reserve(js.size());
0024     for (const Jet& j : js) rtn.push_back(j.pseudojet());
0025     return rtn;
0026   }
0027
0028   inline Jets mkJets(const PseudoJets& pjs) {
0029     Jets rtn; rtn.reserve(pjs.size());
0030     for (const PseudoJet& pj : pjs) rtn.push_back(pj);
0031     return rtn;
0032   }
0033
0034   /// @}
0035
0036
0037   /// @defgroup jetutils_j2bool Jet classifier -> bool functors
0038   /// @{
0039
0040   /// std::function instantiation for functors taking a Jet and returning a bool
0041   using JetSelector = function<bool(const Jet&)>;
0042   /// std::function instantiation for functors taking two Jets and returning a bool
0043   using JetSorter = function<bool(const Jet&, const Jet&)>;
0044
0045
0046   /// Base type for Jet -> bool functors
0047   struct BoolJetFunctor {
0048     virtual bool operator()(const Jet& p) const = 0;
0049     virtual ~BoolJetFunctor() {}
0050   };
0051
0052
0053   /// Functor for and-combination of selector logic
0054   struct BoolJetAND : public BoolJetFunctor {
0055     BoolJetAND(const std::vector<JetSelector>& sels) : selectors(sels) {}
0056     BoolJetAND(const JetSelector& a, const JetSelector& b) : selectors({a,b}) {}
0057     BoolJetAND(const JetSelector& a, const JetSelector& b, const JetSelector& c) : selectors({a,b,c}) {}
0058     bool operator()(const Jet& j) const {
0059       for (const JetSelector& sel : selectors) if (!sel(j)) return false;
0060       return true;
0061     }
0062     std::vector<JetSelector> selectors;
0063   };
0064   /// Operator syntactic sugar for AND construction
0065   inline BoolJetAND operator && (const JetSelector& a, const JetSelector& b) {
0066     return BoolJetAND(a, b);
0067   }
0068
0069
0070   /// Functor for or-combination of selector logic
0071   struct BoolJetOR : public BoolJetFunctor {
0072     BoolJetOR(const std::vector<JetSelector>& sels) : selectors(sels) {}
0073     BoolJetOR(const JetSelector& a, const JetSelector& b) : selectors({a,b}) {}
0074     BoolJetOR(const JetSelector& a, const JetSelector& b, const JetSelector& c) : selectors({a,b,c}) {}
0075     bool operator()(const Jet& j) const {
0076       for (const JetSelector& sel : selectors) if (sel(j)) return true;
0077       return false;
0078     }
0079     std::vector<JetSelector> selectors;
0080   };
0081   /// Operator syntactic sugar for OR construction
0082   inline BoolJetOR operator || (const JetSelector& a, const JetSelector& b) {
0083     return BoolJetOR(a, b);
0084   }
0085
0086
0087   /// Functor for inverting selector logic
0088   struct BoolJetNOT : public BoolJetFunctor {
0089     BoolJetNOT(const JetSelector& sel) : selector(sel) {}
0090     bool operator()(const Jet& j) const { return !selector(j); }
0091     JetSelector selector;
0092   };
0093   /// Operator syntactic sugar for NOT construction
0094   inline BoolJetNOT operator ! (const JetSelector& a) {
0095     return BoolJetNOT(a);
0096   }
0097
0098
0099
0100   /// B-tagging functor, with a tag selection cut as the stored state
0101   struct HasBTag : BoolJetFunctor {
0102     HasBTag(const Cut& c=Cuts::open(), double dR=-1) : cut(c), deltaR(dR) {}
0103     // HasBTag(const std::function<bool(const Jet& j)>& f) : selector(f) {}
0104     bool operator() (const Jet& j) const { return j.bTagged(cut, deltaR); }
0105     // const std::function<bool(const Jet& j)> selector;
0106     const Cut cut;
0107     const double deltaR;
0108   };
0109   using hasBTag = HasBTag;
0110
0111   /// C-tagging functor, with a tag selection cut as the stored state
0112   struct HasCTag : BoolJetFunctor {
0113     HasCTag(const Cut& c=Cuts::open(), double dR=-1) : cut(c), deltaR(dR) {}
0114     // HasCTag(const std::function<bool(const Jet& j)>& f) : selector(f) {}
0115     bool operator() (const Jet& j) const { return j.cTagged(cut, deltaR); }
0116     // const std::function<bool(const Jet& j)> selector;
0117     const Cut cut;
0118     const double deltaR;
0119   };
0120   using hasCTag = HasCTag;
0121
0122   /// Tau-tagging functor, with a tag selection cut as the stored state
0123   struct HasTauTag : BoolJetFunctor {
0124     HasTauTag(const Cut& c=Cuts::open(), double dR=-1) : cut(c), deltaR(dR) {}
0125     // HasTauTag(const std::function<bool(const Jet& j)>& f) : selector(f) {}
0126     bool operator() (const Jet& j) const { return j.tauTagged(cut, deltaR); }
0127     // const std::function<bool(const Jet& j)> selector;
0128     const Cut cut;
0129     const double deltaR;
0130   };
0131   using hasTauTag = HasTauTag;
0132
0133   /// Anti-B/C-tagging functor, with a tag selection cut as the stored state
0134   struct HasNoTag : BoolJetFunctor {
0135     HasNoTag(const Cut& c=Cuts::open(), double dR=-1, bool quarktagsonly=false) : cut(c), deltaR(dR), qtagsonly(quarktagsonly) {}
0136     HasNoTag(const Cut& c=Cuts::open(), bool quarktagsonly=false) : HasNoTag(c, -1, quarktagsonly) {}
0137     // HasNoTag(const std::function<bool(const Jet& j)>& f) : selector(f) {}
0138     bool operator() (const Jet& j) const { return !j.bTagged(cut, deltaR) && !j.cTagged(cut, deltaR) && (qtagsonly || !j.tauTagged(cut, deltaR)); }
0139     // const std::function<bool(const Jet& j)> selector;
0140     const Cut cut;
0141     const double deltaR;
0142     bool qtagsonly;
0143   };
0144   using hasNoTag = HasNoTag;
0145
0146   /// @}
0147
0148
0149   /// @defgroup jetutils_filt Unbound functions for filtering jets
0150   /// @{
0151
0152   /// Filter a jet collection in-place to the subset that passes the supplied Cut
0153   Jets& iselect(Jets& jets, const Cut& c);
0154
0155
0156   /// Filter a jet collection in-place to the subset that passes the supplied Cut
0157   inline Jets select(const Jets& jets, const Cut& c) {
0158     Jets rtn = jets;
0159     return iselect(rtn, c);
0160   }
0161
0162
0163   /// Filter a jet collection in-place to the subset that passes the supplied Cut
0164   inline Jets select(const Jets& jets, const Cut& c, Jets& out) {
0165     out = select(jets, c);
0166     return out;
0167   }
0168
0169
0170
0171   /// Filter a jet collection in-place to the subset that fails the supplied Cut
0172   Jets& idiscard(Jets& jets, const Cut& c);
0173
0174
0175   /// Filter a jet collection in-place to the subset that fails the supplied Cut
0176   inline Jets discard(const Jets& jets, const Cut& c) {
0177     Jets rtn = jets;
0178     return idiscard(rtn, c);
0179   }
0180
0181
0182   /// Filter a jet collection in-place to the subset that fails the supplied Cut
0183   inline Jets discard(const Jets& jets, const Cut& c, Jets& out) {
0184     out = discard(jets, c);
0185     return out;
0186   }
0187
0188   /// @}
0189
0190
0191   /// @defgroup jetutils_trim Trimming operations for Rivet::Jets
0192   /// @{
0193
0194   /// @brief Trim subjets with insufficient pT fraction
0195   ///
0196   /// Removes subjet constituents of pT less than @param frac*(jetpt) in-place.
0197   ///
0198   /// @note Mainly useful for reclustered jets.
0199   ///
0200   /// @warning Keeps all tags(), which might be incorrect, depending on use case.
0201   /// TODO: Fix this behaviour if necessary.
0202   Jet& itrimJetsFrac(Jet &jet, const double frac);
0203
0204   /// @brief Trim subjets with insufficient pT fraction @param frac
0205   ///
0206   /// @note Sorts by optional @param sortFunc, default is sorting by (descending) pT
0207   ///
0208   /// @note Works on (almost) arbritrary containers of Jet
0209   ///
0210   /// TODO: Does not work on set<Jet> because of const restrictions of set. Fix if necessary.
0211   template <
0212     typename CONTAINER,
0213     typename = std::enable_if_t<
0214       is_citerable_v<CONTAINER>,
0215       Jet
0216     >
0217   >
0218   CONTAINER& itrimJetsFrac(CONTAINER &jets, const double frac, const JetSorter &sortFunc=cmpMomByPt){
0219     for ( Jet &jet : jets ) itrimJetsFrac(jet, frac);
0220     isortBy(jets, sortFunc);
0221     return jets;
0222   }
0223
0224   template <typename T, typename U>
0225   std::map<T, U>& itrimJetsFrac(std::map<T, U> &jetMap, const double frac, const JetSorter &sortFunc=cmpMomByPt){
0226     for ( auto &item : jetMap ) itrimJetsFrac(item.second, frac, sortFunc);
0227     return jetMap;
0228   }
0229
0230   /// @brief Trim subjets with insufficient pT fraction
0231   ///
0232   /// Removes subjet constituents of pT less than @param frac*(jetpt) out-of-place.
0233   ///
0234   /// @note Mainly useful for reclustered jets.
0235   ///
0236   /// @warning Keeps all tags(), which might be incorrect, depending on use case.
0237   /// TODO: Fix this behaviour if necessary.
0238   inline Jet trimJetsFrac(const Jet &jet, double frac){
0239     Jet rtn = jet;
0240     return itrimJetsFrac(rtn, frac);
0241   }
0242
0243   /// @brief Trim subjets with insufficient pT fraction
0244   ///
0245   /// @note Sorts by optional @param sortFunc, default is sorting by (descending) pT
0246   ///
0247   /// @note Works on (almost) arbritrary containers of Jet
0248   ///
0249   /// TODO: Does not work on set<Jet> because of const restrictions of set. Fix if necessary.
0250   template <
0251     typename... Args, typename CONTAINER,
0252     typename = std::enable_if_t<
0253       is_citerable_v<CONTAINER>,
0254       Jet
0255     >
0256   >
0257   CONTAINER trimJetsFrac(const CONTAINER &jets, Args&&... args){
0258     CONTAINER rtn = jets;
0259     return itrimJetsFrac(rtn, std::forward<Args>(args)...);
0260   }
0261
0262   template <typename T, typename U, typename... Args>
0263   std::map<T, U> trimJetsFrac(const std::map<T, U> &jetMap, Args&&... args){
0264     std::map<T, U> rtn = jetMap;
0265     return itrimJetsFrac(rtn, std::forward<Args>(args)...);
0266   }
0267
0268   /// @}
0269
0270   /// @defgroup jetutils_filt Operations to apply a FastJet::Filter to PseudoJet(s)
0271   /// @{
0272
0273   /// @brief Apply given FastJet::Filter @param filter to PseudoJet @param pj in-place
0274   ///
0275   /// @note Can be any kind of filtering/ grooming algorithm, e.g. trimming.
0276   PseudoJet& ifilterPseudoJets(PseudoJet &pj, const fastjet::Filter &filter);
0277
0278   /// @brief Apply given FastJet::Filter @param filter to container of PseudoJet @param pjs in-place
0279   ///
0280   /// @note Sorts by optional @param sortFunc, default is sorting by (descending) pT
0281   template <
0282     typename CONTAINER,
0283     typename = std::enable_if_t<
0284       is_citerable_v<CONTAINER>,
0285       PseudoJet
0286     >
0287   >
0288   CONTAINER& ifilterPseudoJets(CONTAINER &pjs, const fastjet::Filter &filter, const JetSorter &sortFunc=cmpMomByPt){
0289     for ( PseudoJet& pj : pjs ) ifilterPseudoJets(pj, filter);
0290     isortBy(pjs, sortFunc);
0291     return pjs;
0292   }
0293
0294   /// @brief Apply given FastJet::Filter @param filter to map @param pjMap in-place
0295   ///
0296   /// @note Sorts by optional @param sortFunc, default is sorting by (descending) pT
0297   template <typename T, typename U, typename... Args>
0298   std::map<T, U>& ifilterPseudoJets(std::map<T, U> &pjMap, Args&&... args){
0299     for ( auto &item : pjMap ) ifilterPseudoJets(item.second, std::forward<Args>(args)...);
0300     return pjMap;
0301   }
0302
0303   /// @brief Apply given FastJet::Filter @param filter to PseudoJet @param pj out-of-place
0304   inline PseudoJet filterPseudoJets(const PseudoJet &pj, const fastjet::Filter &filter){
0305     PseudoJet rtn = pj;
0306     return ifilterPseudoJets(rtn, filter);
0307   }
0308
0309   /// @brief Apply given FastJet::Filter @param filter to container of PseudoJet @param pjs out-of-place
0310   ///
0311   /// @note Sorts by optional @param sortFunc, default is sorting by (descending) pT
0312   template <
0313     typename... Args, typename CONTAINER,
0314     typename = std::enable_if_t<
0315       is_citerable_v<CONTAINER>,
0316       PseudoJet
0317     >
0318   >
0319   CONTAINER filterPseudoJets(const CONTAINER &pjs, Args&&... args){
0320     CONTAINER rtn = pjs;
0321     return ifilterPseudoJets(rtn, std::forward<Args>(args)...);
0322   }
0323
0324   /// @brief Apply given FastJet::Filter @param filter to map @param pjMap out-of-place
0325   ///
0326   /// @note Sorts by optional @param sortFunc, default is sorting by (descending) pT
0327   template <typename T, typename U, typename... Args>
0328   std::map<T, U> filterPseudoJets(const std::map<T, U> &pjMap, Args&&... args){
0329     std::map<T, U> rtn = pjMap;
0330     return ifilterPseudoJets(rtn, std::forward<Args>(args)...);
0331   }
0332   /// @}
0333
0334
0335   /// @defgroup jetutils_coll Operations on collections of Jet
0336   /// @note This can't be done on generic collections of ParticleBase -- thanks, C++ :-/
0337   /// @{
0338   namespace Kin {
0339
0340     /// @todo This shouldn't be necessary, if the sum() function SFINAE picked up be ParticleBase versions...
0341     inline double pT(const Jet& j) {
0342       return j.pT();
0343     }
0344
0345     inline double sumPt(const Jets& js) {
0346       return sum(js, Kin::pT, 0.0);
0347     }
0348
0349     inline FourMomentum sumP4(const Jets& js) {
0350       return sum(js, Kin::p4, FourMomentum());
0351     }
0352
0353     inline Vector3 sumP3(const Jets& js) {
0354       return sum(js, Kin::p3, Vector3());
0355     }
0356
0357     /// @todo Min dPhi, min dR?
0358     /// @todo Isolation routines?
0359
0360   }
0361   /// @}
0362
0363
0364   // Import Kin namespace into Rivet
0365   using namespace Kin;
0366
0367
0368   /// @}
0369
0370 }
0371
0372 #endif