JetMET/plugins/HLTL1TMatchedJetsVBFFilter.h

0001 #ifndef HLTrigger_JetMET_HLTL1TMatchedJetsVBFFilter_h
0002 #define HLTrigger_JetMET_HLTL1TMatchedJetsVBFFilter_h
0003
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "DataFormats/Common/interface/Ref.h"
0012 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0013 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0014 #include "DataFormats/Math/interface/deltaR.h"
0015 #include "HLTrigger/HLTcore/interface/HLTFilter.h"
0016
0017 #include <vector>
0018 #include <utility>
0019 #include <string>
0020
0021 template <typename T>
0022 class HLTL1TMatchedJetsVBFFilter : public HLTFilter {
0023 public:
0024   explicit HLTL1TMatchedJetsVBFFilter(const edm::ParameterSet&);
0025   ~HLTL1TMatchedJetsVBFFilter() override = default;
0026
0027   bool hltFilter(edm::Event& iEvent,
0028                  const edm::EventSetup& iSetup,
0029                  trigger::TriggerFilterObjectWithRefs& filterproduct) const override;
0030
0031   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0032
0033 private:
0034   auto logTrace() const {
0035     auto const& moduleType = moduleDescription().moduleName();
0036     auto const& moduleLabel = moduleDescription().moduleLabel();
0037     return LogTrace(moduleType) << "[" << moduleType << "] (" << moduleLabel << ") ";
0038   }
0039
0040   enum Algorithm { VBF = 0, VBFPlus2CentralJets = 1 };
0041
0042   Algorithm getAlgorithmFromString(std::string const& str) const;
0043   void fillJetIndices(std::vector<unsigned int>& outputJetIndices,
0044                       std::vector<T> const& jets,
0045                       std::vector<unsigned int> const& jetIndices,
0046                       double const pt1,
0047                       double const pt3,
0048                       double const mjj) const;
0049
0050   // input HLT jets
0051   edm::InputTag const jetTag_;
0052   edm::EDGetTokenT<std::vector<T>> const jetToken_;
0053
0054   // matching with L1T jets by delta-R distance
0055   bool const matchJetsByDeltaR_;
0056   double const maxJetDeltaR_, maxJetDeltaR2_;
0057   edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> const l1tJetRefsToken_;
0058
0059   // selection of HLT jets compatible with the VBF topology (see fillJetIndices)
0060   Algorithm const algorithm_;
0061   double const minPt1_;
0062   double const minPt2_;
0063   double const minPt3_;
0064   double const minInvMass_;
0065
0066   // min/max number of selected jets, and triggerType of output triggerRefs
0067   int const minNJets_;
0068   int const maxNJets_;
0069   int const triggerType_;
0070 };
0071
0072 template <typename T>
0073 void HLTL1TMatchedJetsVBFFilter<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0074   edm::ParameterSetDescription desc;
0075   makeHLTFilterDescription(desc);
0076   desc.add<edm::InputTag>("src", edm::InputTag("hltJets"))->setComment("InputTag of input HLT jets");
0077
0078   desc.add<bool>("matchJetsByDeltaR", true)
0079       ->setComment("Enable delta-R matching between HLT jets ('src') and L1T jets ('l1tJetRefs')");
0080   desc.add<double>("maxJetDeltaR", 0.5)
0081       ->setComment(
0082           "Maximum delta-R distance to match HLT jets ('src') and L1T jets ('l1tJetRefs'), "
0083           "used only if 'matchJetsByDeltaR == true'.");
0084   desc.add<edm::InputTag>("l1tJetRefs", edm::InputTag("hltL1sJetSeed"))
0085       ->setComment("InputTag of references to L1T jets");
0086
0087   desc.add<std::string>("algorithm", "VBF")
0088       ->setComment("Keyword to choose the algorithm used to select jets compatible with the VBF topology");
0089   desc.add<double>("minPt1", 110.)
0090       ->setComment("Minimum pT threshold 'pt1' in algorithm to select VBF jets (must respect 'pt2 <= pt3 <= pt1')");
0091   desc.add<double>("minPt2", 35.)->setComment("Minimum pT of all selected VBF jets (must respect 'pt2 <= pt3 <= pt1')");
0092   desc.add<double>("minPt3", 110.)
0093       ->setComment("Minimum pT threshold 'pt3' in algorithm to select VBF jets (must respect 'pt2 <= pt3 <= pt1')");
0094   desc.add<double>("minInvMass", 650.)->setComment("Minimum jet invariant mass");
0095
0096   desc.add<int>("minNJets", 0)->setComment("Minimum number of output jets to accept the event (ignored if negative)");
0097   desc.add<int>("maxNJets", -1)->setComment("Maximum number of output jets to accept the event (ignored if negative)");
0098
0099   desc.add<int>("triggerType", trigger::TriggerJet);
0100
0101   descriptions.setComment(
0102       "This HLTFilter selects events with HLT jets compatible with the VBF topology. "
0103       "These HLT jets can be required to be matched to L1T jets by delta-R distance.");
0104   descriptions.addWithDefaultLabel(desc);
0105 }
0106
0107 template <typename T>
0108 HLTL1TMatchedJetsVBFFilter<T>::HLTL1TMatchedJetsVBFFilter(const edm::ParameterSet& iConfig)
0109     : HLTFilter(iConfig),
0110       jetTag_(iConfig.getParameter<edm::InputTag>("src")),
0111       jetToken_(consumes(jetTag_)),
0112       matchJetsByDeltaR_(iConfig.getParameter<bool>("matchJetsByDeltaR")),
0113       maxJetDeltaR_(iConfig.getParameter<double>("maxJetDeltaR")),
0114       maxJetDeltaR2_(maxJetDeltaR_ * maxJetDeltaR_),
0115       l1tJetRefsToken_(matchJetsByDeltaR_ ? consumes<trigger::TriggerFilterObjectWithRefs>(
0116                                                 iConfig.getParameter<edm::InputTag>("l1tJetRefs"))
0117                                           : edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs>()),
0118       algorithm_(getAlgorithmFromString(iConfig.getParameter<std::string>("algorithm"))),
0119       minPt1_(iConfig.getParameter<double>("minPt1")),
0120       minPt2_(iConfig.getParameter<double>("minPt2")),
0121       minPt3_(iConfig.getParameter<double>("minPt3")),
0122       minInvMass_(iConfig.getParameter<double>("minInvMass")),
0123       minNJets_(iConfig.getParameter<int>("minNJets")),
0124       maxNJets_(iConfig.getParameter<int>("maxNJets")),
0125       triggerType_(iConfig.getParameter<int>("triggerType")) {
0126   // delta-R matching
0127   if (matchJetsByDeltaR_ and maxJetDeltaR_ < 0) {
0128     throw cms::Exception("InvalidConfigurationParameter")
0129         << "invalid value for parameter \"maxJetDeltaR\" (must be > 0): " << maxJetDeltaR_;
0130   }
0131
0132   // values of minPt{1,2,3} thresholds
0133   if (minPt2_ > minPt1_) {
0134     throw cms::Exception("InvalidConfigurationParameter")
0135         << "minPt1 (" << minPt1_ << ") must not be smaller than minPt2 (" << minPt2_ << ")";
0136   } else if (minPt3_ > minPt1_) {
0137     throw cms::Exception("InvalidConfigurationParameter")
0138         << "minPt1 (" << minPt1_ << ") must not be smaller than minPt3 (" << minPt3_ << ")";
0139   } else if (minPt2_ > minPt3_) {
0140     throw cms::Exception("InvalidConfigurationParameter")
0141         << "minPt3 (" << minPt3_ << ") must not be smaller than minPt2 (" << minPt2_ << ")";
0142   }
0143 }
0144
0145 template <typename T>
0146 typename HLTL1TMatchedJetsVBFFilter<T>::Algorithm HLTL1TMatchedJetsVBFFilter<T>::getAlgorithmFromString(
0147     std::string const& str) const {
0148   if (str == "VBF") {
0149     return Algorithm::VBF;
0150   } else if (str == "VBFPlus2CentralJets") {
0151     return Algorithm::VBFPlus2CentralJets;
0152   }
0153
0154   throw cms::Exception("HLTL1TMatchedJetsVBFFilterInvalidAlgorithmName")
0155       << "invalid value for argument of getAlgorithmFromString method: " << str
0156       << " (valid values are \"VBF\" and \"VBFPlus2CentralJets\")";
0157 }
0158
0159 template <typename T>
0160 void HLTL1TMatchedJetsVBFFilter<T>::fillJetIndices(std::vector<unsigned int>& outputJetIndices,
0161                                                    std::vector<T> const& jets,
0162                                                    std::vector<unsigned int> const& jetIndices,
0163                                                    double const pt1,
0164                                                    double const pt3,
0165                                                    double const mjj) const {
0166   // reset output jet indices
0167   outputJetIndices.clear();
0168
0169   // find dijet pair with highest Mjj
0170   int i1 = -1;
0171   int i2 = -1;
0172   double m2jj_max = -1;
0173
0174   for (unsigned int i = 0; i < jetIndices.size(); i++) {
0175     auto const& jet1 = jets[jetIndices[i]];
0176
0177     for (unsigned int j = i + 1; j < jetIndices.size(); j++) {
0178       auto const& jet2 = jets[jetIndices[j]];
0179
0180       double const m2jj_tmp = (jet1.p4() + jet2.p4()).M2();
0181
0182       if (m2jj_tmp > m2jj_max) {
0183         m2jj_max = m2jj_tmp;
0184         i1 = jetIndices[i];
0185         i2 = jetIndices[j];
0186       }
0187     }
0188   }
0189
0190   // only proceed if highest-Mjj dijet pair was found
0191   if (i1 >= 0 and i2 >= 0) {
0192     logTrace() << "dijet pair with highest invariant mass: indices={" << i1 << "," << i2 << "} mjj^2=" << m2jj_max;
0193
0194     unsigned int const j1 = i1;
0195     unsigned int const j2 = i2;
0196
0197     auto const& jet1 = jets[j1];
0198     double const m2jj = (mjj >= 0. ? mjj * mjj : -1.);
0199
0200     // algorithm: "VBF"
0201     if (algorithm_ == Algorithm::VBF) {
0202       if (m2jj_max > m2jj) {
0203         if (jet1.pt() >= pt1) {
0204           outputJetIndices = {j1, j2};
0205         } else if (jet1.pt() < pt3 and jets[jetIndices[0]].pt() > pt3) {
0206           outputJetIndices = {j1, j2, jetIndices[0]};
0207         }
0208       }
0209     } else if (algorithm_ == Algorithm::VBFPlus2CentralJets) {
0210       if (jetIndices.size() > 3 and m2jj_max > m2jj) {
0211         // indices of additional jets
0212         std::vector<unsigned int> idx_jets;
0213         idx_jets.reserve(jetIndices.size() - 2);
0214
0215         for (auto const idx : jetIndices) {
0216           if (idx == j1 or idx == j2)
0217             continue;
0218           idx_jets.emplace_back(idx);
0219         }
0220
0221         if (jet1.pt() >= pt3) {
0222           outputJetIndices.emplace_back(j1);
0223           outputJetIndices.emplace_back(j2);
0224
0225           for (auto const idx : idx_jets) {
0226             if (outputJetIndices.size() > 3)
0227               break;
0228             outputJetIndices.emplace_back(idx);
0229           }
0230         } else if (idx_jets.size() >= 3) {
0231           outputJetIndices.emplace_back(j1);
0232           outputJetIndices.emplace_back(j2);
0233           outputJetIndices.emplace_back(idx_jets[0]);
0234           outputJetIndices.emplace_back(idx_jets[1]);
0235           outputJetIndices.emplace_back(idx_jets[2]);
0236
0237           if (idx_jets.size() > 3) {
0238             outputJetIndices.emplace_back(idx_jets[3]);
0239           }
0240         }
0241       }
0242     }
0243   }
0244 }
0245
0246 template <typename T>
0247 bool HLTL1TMatchedJetsVBFFilter<T>::hltFilter(edm::Event& iEvent,
0248                                               const edm::EventSetup&,
0249                                               trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0250   if (saveTags())
0251     filterproduct.addCollectionTag(jetTag_);
0252
0253   // handle to input HLT jets
0254   auto const& jetsHandle = iEvent.getHandle(jetToken_);
0255
0256   // indices of input HLT jets satisfying delta-R matching requirements
0257   std::vector<unsigned int> jetIndices;
0258
0259   if (matchJetsByDeltaR_) {
0260     // refs to L1T jets used for delta-R matching
0261     l1t::JetVectorRef l1tJetRefs;
0262     iEvent.get(l1tJetRefsToken_).getObjects(trigger::TriggerL1Jet, l1tJetRefs);
0263
0264     jetIndices.reserve(jetsHandle->size());
0265     for (unsigned int iJet = 0; iJet < jetsHandle->size(); ++iJet) {
0266       auto const& jet = (*jetsHandle)[iJet];
0267
0268       if (jet.pt() <= minPt2_)
0269         continue;
0270
0271       for (unsigned int l1tJetIdx = 0; l1tJetIdx < l1tJetRefs.size(); ++l1tJetIdx) {
0272         if (reco::deltaR2(jet.p4(), l1tJetRefs[l1tJetIdx]->p4()) < maxJetDeltaR2_) {
0273           logTrace() << "input jet: index=" << iJet << " pt=" << jet.pt() << " eta=" << jet.eta()
0274                      << " phi=" << jet.phi() << " mass=" << jet.mass()
0275                      << " (matched to L1T jet with pt=" << l1tJetRefs[l1tJetIdx]->pt()
0276                      << " eta=" << l1tJetRefs[l1tJetIdx]->eta() << " phi=" << l1tJetRefs[l1tJetIdx]->phi()
0277                      << " mass=" << l1tJetRefs[l1tJetIdx]->mass() << ")";
0278
0279           jetIndices.emplace_back(iJet);
0280           break;
0281         }
0282       }
0283     }
0284   } else {
0285     jetIndices.reserve(jetsHandle->size());
0286     for (unsigned int iJet = 0; iJet < jetsHandle->size(); ++iJet) {
0287       auto const& jet = (*jetsHandle)[iJet];
0288
0289       if (jet.pt() <= minPt2_)
0290         continue;
0291
0292       logTrace() << "input jet: index=" << iJet << " pt=" << jet.pt() << " eta=" << jet.eta() << " phi=" << jet.phi()
0293                  << " mass=" << jet.mass();
0294
0295       jetIndices.emplace_back(iJet);
0296     }
0297   }
0298
0299   // order jet indices by jet pT (highest pT first)
0300   std::sort(jetIndices.begin(), jetIndices.end(), [&jetsHandle](unsigned int const i1, unsigned int const i2) {
0301     return (*jetsHandle)[i1].pt() > (*jetsHandle)[i2].pt();
0302   });
0303
0304   // indices of jets identified as compatible with the VBF topology
0305   std::vector<unsigned int> outputJetIndices;
0306   fillJetIndices(outputJetIndices, *jetsHandle, jetIndices, minPt1_, minPt3_, minInvMass_);
0307
0308   // add selected jets to trigger::TriggerFilterObjectWithRefs
0309   for (auto const idx : outputJetIndices) {
0310     filterproduct.addObject(triggerType_, edm::Ref<std::vector<T>>(jetsHandle, idx));
0311
0312     logTrace() << "output jet: index=" << idx << " pt=" << (*jetsHandle)[idx].pt()
0313                << " eta=" << (*jetsHandle)[idx].eta() << " phi=" << (*jetsHandle)[idx].phi()
0314                << " mass=" << (*jetsHandle)[idx].mass();
0315   }
0316
0317   // accept event only if minNJets and maxNJets requirements are satisfied
0318   auto const ret = ((minNJets_ < 0 or outputJetIndices.size() >= (unsigned int)minNJets_) and
0319                     (maxNJets_ < 0 or outputJetIndices.size() <= (unsigned int)maxNJets_));
0320
0321   logTrace() << "filter decision = " << ret << " (trigger objects = " << outputJetIndices.size() << ")";
0322
0323   return ret;
0324 }
0325
0326 #endif