HLTProducers/interface/L1TJetsMatching.h

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoTauTag/HLTProducers
0004 // Class:      L1TJetsMatching
0005 //
0006 /**\class L1TJetsMatching L1TJetsMatching.h
0007  RecoTauTag/HLTProducers/interface/L1TJetsMatching.h
0008  Description:
0009  Matching L1 to PF/Calo Jets. Used for HLT_VBF paths.
0010     *Matches PF/Calo Jets to L1 jets from the dedicated seed
0011     *Adds selection criteria to the leading/subleading jets as well as the maximum dijet mass
0012     *Separates collections of PF/Calo jets into two categories
0013
0014
0015 */
0016 //
0017 // Original Author:  Vukasin Milosevic
0018 //         Created:  Thu, 01 Jun 2017 17:23:00 GMT
0019 //
0020 //
0021
0022 #ifndef RecoTauTag_HLTProducers_L1TJetsMatching_h
0023 #define RecoTauTag_HLTProducers_L1TJetsMatching_h
0024
0025 // user include files
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/global/EDProducer.h"
0028 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/Utilities/interface/InputTag.h"
0033 #include "DataFormats/Common/interface/Handle.h"
0034 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0035 #include "DataFormats/TauReco/interface/PFTauFwd.h"
0036 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0037
0038 #include "Math/GenVector/VectorUtil.h"
0039 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0040 #include "FWCore/Utilities/interface/Exception.h"
0041 #include "DataFormats/JetReco/interface/PFJet.h"
0042
0043 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0044 #include "DataFormats/Math/interface/deltaR.h"
0045
0046 #include <vector>
0047 #include <utility>
0048 #include <tuple>
0049 #include <string>
0050
0051 template <typename T>
0052 class L1TJetsMatching : public edm::global::EDProducer<> {
0053 public:
0054   explicit L1TJetsMatching(const edm::ParameterSet&);
0055   ~L1TJetsMatching() override = default;
0056   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0057   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0058
0059 private:
0060   std::pair<std::vector<T>, std::vector<T>> categorise(
0061       const std::vector<T>& pfMatchedJets, double pt1, double pt2, double pt3, double Mjj) const;
0062   std::tuple<std::vector<T>, std::vector<T>, std::vector<T>> categoriseVBFPlus2CentralJets(
0063       const std::vector<T>& pfMatchedJets, double pt1, double pt2, double pt3, double Mjj) const;
0064   const edm::EDGetTokenT<std::vector<T>> jetSrc_;
0065   const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> jetTrigger_;
0066   const std::string matchingMode_;
0067   const double pt1Min_;
0068   const double pt2Min_;
0069   const double pt3Min_;
0070   const double mjjMin_;
0071   const double matchingR2_;
0072 };
0073 //
0074 // class declaration
0075 //
0076 template <typename T>
0077 std::pair<std::vector<T>, std::vector<T>> L1TJetsMatching<T>::categorise(
0078     const std::vector<T>& pfMatchedJets, double pt1, double pt2, double pt3, double Mjj) const {
0079   std::pair<std::vector<T>, std::vector<T>> output;
0080   unsigned int i1 = 0;
0081   unsigned int i2 = 0;
0082   double m2jj = 0;
0083   if (pfMatchedJets.size() > 1) {
0084     for (unsigned int i = 0; i < pfMatchedJets.size() - 1; i++) {
0085       const T& myJet1 = (pfMatchedJets)[i];
0086
0087       for (unsigned int j = i + 1; j < pfMatchedJets.size(); j++) {
0088         const T& myJet2 = (pfMatchedJets)[j];
0089
0090         const double m2jj_test = (myJet1.p4() + myJet2.p4()).M2();
0091
0092         if (m2jj_test > m2jj) {
0093           m2jj = m2jj_test;
0094           i1 = i;
0095           i2 = j;
0096         }
0097       }
0098     }
0099
0100     const T& myJet1 = (pfMatchedJets)[i1];
0101     const T& myJet2 = (pfMatchedJets)[i2];
0102     const double M2jj = (Mjj >= 0. ? Mjj * Mjj : -1.);
0103
0104     if ((m2jj > M2jj) && (myJet1.pt() >= pt1) && (myJet2.pt() > pt2)) {
0105       output.first.push_back(myJet1);
0106       output.first.push_back(myJet2);
0107     }
0108
0109     if ((m2jj > M2jj) && (myJet1.pt() < pt3) && (myJet1.pt() > pt2) && (myJet2.pt() > pt2)) {
0110       const T& myJetTest = (pfMatchedJets)[0];
0111       if (myJetTest.pt() > pt3) {
0112         output.second.push_back(myJet1);
0113         output.second.push_back(myJet2);
0114         output.second.push_back(myJetTest);
0115       }
0116     }
0117   }
0118
0119   return output;
0120 }
0121 template <typename T>
0122 std::tuple<std::vector<T>, std::vector<T>, std::vector<T>> L1TJetsMatching<T>::categoriseVBFPlus2CentralJets(
0123     const std::vector<T>& pfMatchedJets, double pt1, double pt2, double pt3, double Mjj) const {  //60, 30, 50, 500
0124   std::tuple<std::vector<T>, std::vector<T>, std::vector<T>> output;
0125   unsigned int i1 = 0;
0126   unsigned int i2 = 0;
0127
0128   double m2jj = 0;
0129   if (pfMatchedJets.size() > 1) {
0130     for (unsigned int i = 0; i < pfMatchedJets.size() - 1; i++) {
0131       const T& myJet1 = (pfMatchedJets)[i];
0132
0133       for (unsigned int j = i + 1; j < pfMatchedJets.size(); j++) {
0134         const T& myJet2 = (pfMatchedJets)[j];
0135
0136         const double m2jj_test = (myJet1.p4() + myJet2.p4()).M2();
0137
0138         if (m2jj_test > m2jj) {
0139           m2jj = m2jj_test;
0140           i1 = i;
0141           i2 = j;
0142         }
0143       }
0144     }
0145
0146     const T& myJet1 = (pfMatchedJets)[i1];
0147     const T& myJet2 = (pfMatchedJets)[i2];
0148     const double M2jj = (Mjj >= 0. ? Mjj * Mjj : -1.);
0149
0150     std::vector<T> vec4jets;
0151     vec4jets.reserve(4);
0152     std::vector<T> vec5jets;
0153     vec5jets.reserve(5);
0154     std::vector<T> vec6jets;
0155     vec6jets.reserve(6);
0156     if (pfMatchedJets.size() > 3) {
0157       if ((m2jj > M2jj) && (myJet1.pt() >= pt3) && (myJet2.pt() > pt2)) {
0158         vec4jets.push_back(myJet1);
0159         vec4jets.push_back(myJet2);
0160
0161         for (unsigned int i = 0; i < pfMatchedJets.size(); i++) {
0162           if (vec4jets.size() > 3)
0163             break;
0164           if (i == i1 or i == i2)
0165             continue;
0166           vec4jets.push_back(pfMatchedJets[i]);
0167         }
0168       }
0169
0170       if ((m2jj > M2jj) && (myJet1.pt() < pt1) && (myJet1.pt() < pt3) && (myJet1.pt() > pt2) &&
0171           (myJet2.pt() > pt2)) {  //60, 30, 50, 500
0172
0173         std::vector<unsigned int> idx_jets;
0174         idx_jets.reserve(pfMatchedJets.size() - 2);
0175
0176         for (unsigned int i = 0; i < pfMatchedJets.size(); i++) {
0177           if (i == i1 || i == i2)
0178             continue;
0179           if (pfMatchedJets[i].pt() > pt2) {
0180             idx_jets.push_back(i);
0181           }
0182         }
0183         if (idx_jets.size() == 3) {
0184           vec5jets.push_back(myJet1);
0185           vec5jets.push_back(myJet2);
0186           vec5jets.push_back(pfMatchedJets[idx_jets[0]]);
0187           vec5jets.push_back(pfMatchedJets[idx_jets[1]]);
0188           vec5jets.push_back(pfMatchedJets[idx_jets[2]]);
0189
0190         } else if (idx_jets.size() > 3) {
0191           vec6jets.push_back(myJet1);
0192           vec6jets.push_back(myJet2);
0193           vec6jets.push_back(pfMatchedJets[idx_jets[0]]);
0194           vec6jets.push_back(pfMatchedJets[idx_jets[1]]);
0195           vec6jets.push_back(pfMatchedJets[idx_jets[2]]);
0196           vec6jets.push_back(pfMatchedJets[idx_jets[3]]);
0197         }
0198       }
0199     }
0200
0201     output = std::make_tuple(vec4jets, vec5jets, vec6jets);
0202   }
0203
0204   return output;
0205 }
0206
0207 template <typename T>
0208 L1TJetsMatching<T>::L1TJetsMatching(const edm::ParameterSet& iConfig)
0209     : jetSrc_(consumes<std::vector<T>>(iConfig.getParameter<edm::InputTag>("JetSrc"))),
0210       jetTrigger_(consumes<trigger::TriggerFilterObjectWithRefs>(iConfig.getParameter<edm::InputTag>("L1JetTrigger"))),
0211       matchingMode_(iConfig.getParameter<std::string>("matchingMode")),
0212       pt1Min_(iConfig.getParameter<double>("pt1Min")),
0213       pt2Min_(iConfig.getParameter<double>("pt2Min")),
0214       pt3Min_(iConfig.getParameter<double>("pt3Min")),
0215       mjjMin_(iConfig.getParameter<double>("mjjMin")),
0216       matchingR2_(iConfig.getParameter<double>("matchingR") * iConfig.getParameter<double>("matchingR")) {
0217   if (matchingMode_ == "VBF") {  // Default
0218     produces<std::vector<T>>("TwoJets");
0219     produces<std::vector<T>>("ThreeJets");
0220   } else if (matchingMode_ == "VBFPlus2CentralJets") {
0221     produces<std::vector<T>>("FourJets");
0222     produces<std::vector<T>>("FiveJets");
0223     produces<std::vector<T>>("SixJets");
0224   } else {
0225     throw cms::Exception("InvalidConfiguration") << "invalid value for \"matchingMode\": " << matchingMode_
0226                                                  << " (valid values are \"VBF\" and \"VBFPlus2CentralJets\")";
0227   }
0228 }
0229
0230 template <typename T>
0231 void L1TJetsMatching<T>::produce(edm::StreamID iSId, edm::Event& iEvent, const edm::EventSetup& iES) const {
0232   unique_ptr<std::vector<T>> pfMatchedJets(new std::vector<T>);
0233
0234   // Getting HLT jets to be matched
0235   edm::Handle<std::vector<T>> pfJets;
0236   iEvent.getByToken(jetSrc_, pfJets);
0237
0238   edm::Handle<trigger::TriggerFilterObjectWithRefs> l1TriggeredJets;
0239   iEvent.getByToken(jetTrigger_, l1TriggeredJets);
0240
0241   //l1t::TauVectorRef jetCandRefVec;
0242   l1t::JetVectorRef jetCandRefVec;
0243   l1TriggeredJets->getObjects(trigger::TriggerL1Jet, jetCandRefVec);
0244
0245   math::XYZPoint a(0., 0., 0.);
0246
0247   //std::cout<<"PFsize= "<<pfJets->size()<<endl<<" L1size= "<<jetCandRefVec.size()<<std::endl;
0248   for (unsigned int iJet = 0; iJet < pfJets->size(); iJet++) {
0249     const T& myJet = (*pfJets)[iJet];
0250     for (unsigned int iL1Jet = 0; iL1Jet < jetCandRefVec.size(); iL1Jet++) {
0251       // Find the relative L2pfJets, to see if it has been reconstructed
0252       //  if ((iJet<3) && (iL1Jet==0))  std::cout<<myJet.p4().Pt()<<" ";
0253       if ((reco::deltaR2(myJet.p4(), jetCandRefVec[iL1Jet]->p4()) < matchingR2_) && (myJet.pt() > pt2Min_)) {
0254         pfMatchedJets->push_back(myJet);
0255         break;
0256       }
0257     }
0258   }
0259   // order pfMatchedJets by pT
0260   std::sort(pfMatchedJets->begin(), pfMatchedJets->end(), [](const T& j1, const T& j2) { return j1.pt() > j2.pt(); });
0261
0262   if (matchingMode_ == "VBF") {  // Default
0263     std::pair<std::vector<T>, std::vector<T>> output = categorise(*pfMatchedJets, pt1Min_, pt2Min_, pt3Min_, mjjMin_);
0264     auto output1 = std::make_unique<std::vector<T>>(output.first);
0265     auto output2 = std::make_unique<std::vector<T>>(output.second);
0266
0267     iEvent.put(std::move(output1), "TwoJets");
0268     iEvent.put(std::move(output2), "ThreeJets");
0269
0270   } else if (matchingMode_ == "VBFPlus2CentralJets") {
0271     std::tuple<std::vector<T>, std::vector<T>, std::vector<T>> output =
0272         categoriseVBFPlus2CentralJets(*pfMatchedJets, pt1Min_, pt2Min_, pt3Min_, mjjMin_);
0273     auto output1 = std::make_unique<std::vector<T>>(std::get<0>(output));
0274     auto output2 = std::make_unique<std::vector<T>>(std::get<1>(output));
0275     auto output3 = std::make_unique<std::vector<T>>(std::get<2>(output));
0276
0277     iEvent.put(std::move(output1), "FourJets");
0278     iEvent.put(std::move(output2), "FiveJets");
0279     iEvent.put(std::move(output3), "SixJets");
0280   }
0281 }
0282
0283 template <typename T>
0284 void L1TJetsMatching<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0285   edm::ParameterSetDescription desc;
0286   desc.add<edm::InputTag>("L1JetTrigger", edm::InputTag("hltL1DiJetVBF"))->setComment("Name of trigger filter");
0287   desc.add<edm::InputTag>("JetSrc", edm::InputTag("hltAK4PFJetsTightIDCorrected"))
0288       ->setComment("Input collection of PFJets");
0289   desc.add<std::string>("matchingMode", "VBF")
0290       ->setComment("Switch from Di/tri-jet (VBF) to Multi-jet (VBFPlus2CentralJets) matching");
0291   desc.add<double>("pt1Min", 110.0)->setComment("Minimal pT1 of PFJets to match");
0292   desc.add<double>("pt2Min", 35.0)->setComment("Minimal pT2 of PFJets to match");
0293   desc.add<double>("pt3Min", 110.0)->setComment("Minimum pT3 of PFJets to match");
0294   desc.add<double>("mjjMin", 650.0)->setComment("Minimal mjj of matched PFjets");
0295   desc.add<double>("matchingR", 0.5)->setComment("dR value used for matching");
0296   descriptions.setComment(
0297       "This module produces collection of PFJets matched to L1 Taus / Jets passing a HLT filter (Only p4 and vertex "
0298       "of returned PFJets are set).");
0299   descriptions.add(defaultModuleLabel<L1TJetsMatching<T>>(), desc);
0300 }
0301
0302 #endif