Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:20:32

0001 /** \class HLTJetSortedVBFFilter
0002  *
0003  * See header file for documentation
0004  *
0005 
0006 
0007  *
0008  *  \author Jacopo Bernardini
0009  *
0010  */
0011 
0012 #include <vector>
0013 #include <string>
0014 
0015 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "DataFormats/Math/interface/deltaPhi.h"
0020 #include "DataFormats/Math/interface/deltaR.h"
0021 #include "DataFormats/Common/interface/RefToBase.h"
0022 #include "HLTrigger/JetMET/interface/HLTJetSortedVBFFilter.h"
0023 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0024 
0025 using namespace std;
0026 //
0027 // constructors and destructor//
0028 //
0029 template <typename T>
0030 HLTJetSortedVBFFilter<T>::HLTJetSortedVBFFilter(const edm::ParameterSet& iConfig)
0031     : HLTFilter(iConfig),
0032       inputJets_(iConfig.getParameter<edm::InputTag>("inputJets")),
0033       inputJetTags_(iConfig.getParameter<edm::InputTag>("inputJetTags")),
0034       mqq_(iConfig.getParameter<double>("Mqq")),
0035       detaqq_(iConfig.getParameter<double>("Detaqq")),
0036       detabb_(iConfig.getParameter<double>("Detabb")),
0037       dphibb_(iConfig.getParameter<double>("Dphibb")),
0038       ptsqq_(iConfig.getParameter<double>("Ptsumqq")),
0039       ptsbb_(iConfig.getParameter<double>("Ptsumbb")),
0040       seta_(iConfig.getParameter<double>("Etaq1Etaq2")),
0041       njets_(iConfig.getParameter<int>("njets")),
0042       value_(iConfig.getParameter<std::string>("value")),
0043       triggerType_(iConfig.getParameter<int>("triggerType")) {
0044   m_theJetsToken = consumes<std::vector<T>>(inputJets_);
0045   m_theJetTagsToken = consumes<reco::JetTagCollection>(inputJetTags_);
0046   if (njets_ < 4) {
0047     edm::LogWarning("LowNJets") << "njets=" << njets_ << " it must be >=4. Forced njets=4.";
0048     njets_ = 4;
0049   }
0050 }
0051 
0052 template <typename T>
0053 HLTJetSortedVBFFilter<T>::~HLTJetSortedVBFFilter() = default;
0054 
0055 template <typename T>
0056 void HLTJetSortedVBFFilter<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0057   edm::ParameterSetDescription desc;
0058   makeHLTFilterDescription(desc);
0059   desc.add<edm::InputTag>("inputJets", edm::InputTag("hltJetCollection"));
0060   desc.add<edm::InputTag>("inputJetTags", edm::InputTag(""));
0061   desc.add<double>("Mqq", 200);
0062   desc.add<double>("Detaqq", 2.5);
0063   desc.add<double>("Detabb", 10.);
0064   desc.add<double>("Dphibb", 10.);
0065   desc.add<double>("Ptsumqq", 0.);
0066   desc.add<double>("Ptsumbb", 0.);
0067   desc.add<double>("Etaq1Etaq2", 40.);
0068   desc.add<std::string>("value", "second");
0069   desc.add<int>("triggerType", trigger::TriggerJet);
0070   desc.add<int>("njets", 4);
0071   descriptions.add(defaultModuleLabel<HLTJetSortedVBFFilter<T>>(), desc);
0072 }
0073 
0074 template <typename T>
0075 float HLTJetSortedVBFFilter<T>::findCSV(const typename std::vector<T>::const_iterator& jet,
0076                                         const reco::JetTagCollection& jetTags) {
0077   float minDr = 0.1;
0078   float tmpCSV = -20;
0079   for (auto jetb = jetTags.begin(); (jetb != jetTags.end()); ++jetb) {
0080     float tmpDr = reco::deltaR(*jet, *(jetb->first));
0081 
0082     if (tmpDr < minDr) {
0083       minDr = tmpDr;
0084       tmpCSV = jetb->second;
0085     }
0086   }
0087   return tmpCSV;
0088 }
0089 //
0090 // member functions
0091 //
0092 
0093 // ------------ method called to produce the data  ------------
0094 template <typename T>
0095 bool HLTJetSortedVBFFilter<T>::hltFilter(edm::Event& event,
0096                                          const edm::EventSetup& setup,
0097                                          trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0098   using namespace std;
0099   using namespace edm;
0100   using namespace reco;
0101   using namespace trigger;
0102 
0103   typedef vector<T> TCollection;
0104   typedef Ref<TCollection> TRef;
0105 
0106   bool accept(false);
0107 
0108   Handle<TCollection> jets;
0109   event.getByToken(m_theJetsToken, jets);
0110   Handle<JetTagCollection> jetTags;
0111 
0112   if (saveTags())
0113     filterproduct.addCollectionTag(inputJets_);
0114   if (jets->size() < 4)
0115     return false;
0116 
0117   const unsigned int nMax(njets_ < jets->size() ? njets_ : jets->size());
0118   vector<Jpair> sorted(nMax);
0119   vector<TRef> jetRefs(nMax);
0120   unsigned int nJetRefs = 0;
0121 
0122   unsigned int nJet = 0;
0123   double value(0.0);
0124 
0125   Particle::LorentzVector b1, b2, q1, q2;
0126   if (inputJetTags_.encode().empty()) {
0127     for (typename TCollection::const_iterator jet = jets->begin(); (jet != jets->end() && nJet < nMax); ++jet) {
0128       if (value_ == "Pt") {
0129         value = jet->pt();
0130       } else if (value_ == "Eta") {
0131         value = jet->eta();
0132       } else if (value_ == "Phi") {
0133         value = jet->phi();
0134       } else {
0135         value = 0.0;
0136       }
0137       sorted[nJet] = make_pair(value, nJet);
0138       ++nJet;
0139     }
0140     sort(sorted.begin(), sorted.end(), comparator);
0141     for (unsigned int i = 0; i < nMax; ++i) {
0142       jetRefs[i] = TRef(jets, sorted[i].second);
0143     }
0144     nJetRefs = nMax;
0145     q1 = jetRefs[3]->p4();
0146     b1 = jetRefs[2]->p4();
0147     b2 = jetRefs[1]->p4();
0148     q2 = jetRefs[0]->p4();
0149   } else if (value_ == "1BTagAndEta") {
0150     event.getByToken(m_theJetTagsToken, jetTags);
0151     vector<Jpair> sorted;
0152     unsigned int b1_idx = -1;
0153     float csv_max = -999;
0154     for (typename TCollection::const_iterator jet = jets->begin(); (jet != jets->end() && nJet < nMax);
0155          ++jet) {  //fill "sorted" and get the most b-tagged jet with higher CSV (b1)
0156       value = findCSV(jet, *jetTags);
0157       if (value > csv_max) {
0158         csv_max = value;
0159         b1_idx = nJet;
0160       }
0161       sorted.push_back(make_pair(jet->eta(), nJet));
0162       nJet++;
0163       //        cout << "jetPt=" << jet->pt() << "\tjetEta=" << jet->eta() << "\tjetCSV=" << value << endl;
0164     }
0165     if (b1_idx >= sorted.size())
0166       edm::LogError("OutOfRange") << "b1 index out of range.";
0167     sorted.erase(sorted.begin() + b1_idx);           //remove the most b-tagged jet from "sorted"
0168     sort(sorted.begin(), sorted.end(), comparator);  //sort "sorted" by eta
0169 
0170     unsigned int q1_idx = sorted.front().second;  //take the backward jet (q1)
0171     unsigned int q2_idx = sorted.back().second;   //take the forward jet (q2)
0172 
0173     unsigned int i = 0;
0174     while ((i == q1_idx) || (i == q2_idx) || (i == b1_idx))
0175       i++;  //take jet with highest pT but q1,q2,b1 (q2)
0176     unsigned int b2_idx = i;
0177 
0178     if (q1_idx < jets->size())
0179       q1 = jets->at(q1_idx).p4();
0180     else
0181       edm::LogWarning("Something wrong with q1");
0182     if (q2_idx < jets->size())
0183       q2 = jets->at(q2_idx).p4();
0184     else
0185       edm::LogWarning("Something wrong with q2");
0186     if (b1_idx < jets->size())
0187       b1 = jets->at(b1_idx).p4();
0188     else
0189       edm::LogWarning("Something wrong with b1");
0190     if (b2_idx < jets->size())
0191       b2 = jets->at(b2_idx).p4();
0192     else
0193       edm::LogWarning("Something wrong with b2");
0194 
0195     jetRefs[0] = TRef(jets, b1_idx);
0196     jetRefs[1] = TRef(jets, b2_idx);
0197     jetRefs[2] = TRef(jets, q1_idx);
0198     jetRefs[3] = TRef(jets, q2_idx);
0199     nJetRefs = 4;
0200 
0201     //      cout<<"\tPathB: b1="<<b1.pt()<<" b2="<<b2.pt()<<" q1="<<q1.pt()<<" q2="<<q2.pt()<<endl;
0202   } else if (value_ == "2BTagAndPt") {
0203     event.getByToken(m_theJetTagsToken, jetTags);
0204     vector<Jpair> sorted;
0205 
0206     unsigned int b1_idx = -1;
0207     unsigned int b2_idx = -1;
0208     float csv1 = -999;
0209     float csv2 = -999;
0210     for (typename TCollection::const_iterator jet = jets->begin(); (jet != jets->end() && nJet < nMax);
0211          ++jet) {  //fill "sorted" and get the two most b-tagged jets (b1,b2)
0212       value = findCSV(jet, *jetTags);
0213       if (value > csv1) {
0214         csv2 = csv1;
0215         b2_idx = b1_idx;
0216         csv1 = value;
0217         b1_idx = nJet;
0218       } else if (value > csv2) {
0219         csv2 = value;
0220         b2_idx = nJet;
0221       }
0222       sorted.push_back(make_pair(jet->eta(), nJet));
0223       nJet++;
0224       //        cout << "jetPt=" << jet->pt() << "\tjetEta=" << jet->eta() << "\tjetCSV=" << value << endl;
0225     }
0226     sorted.erase(sorted.begin() + b1_idx);  //remove b1 and b2 from sorted
0227     sorted.erase(sorted.begin() + (b1_idx > b2_idx ? b2_idx : b2_idx - 1));
0228 
0229     unsigned int q1_idx = sorted.at(0).second;  //get q1 and q2 as the jets with highest pT, but b1 and b2.
0230     unsigned int q2_idx = sorted.at(1).second;
0231 
0232     if (q1_idx < jets->size())
0233       q1 = jets->at(q1_idx).p4();
0234     else
0235       edm::LogWarning("Something wrong with q1");
0236     if (q2_idx < jets->size())
0237       q2 = jets->at(q2_idx).p4();
0238     else
0239       edm::LogWarning("Something wrong with q2");
0240     if (b1_idx < jets->size())
0241       b1 = jets->at(b1_idx).p4();
0242     else
0243       edm::LogWarning("Something wrong with b1");
0244     if (b2_idx < jets->size())
0245       b2 = jets->at(b2_idx).p4();
0246     else
0247       edm::LogWarning("Something wrong with b2");
0248 
0249     jetRefs[0] = TRef(jets, b1_idx);
0250     jetRefs[1] = TRef(jets, b2_idx);
0251     jetRefs[2] = TRef(jets, q1_idx);
0252     jetRefs[3] = TRef(jets, q2_idx);
0253     nJetRefs = 4;
0254 
0255     //      cout<<"\tPathA: b1="<<b1.pt()<<" b2="<<b2.pt()<<" q1="<<q1.pt()<<" q2="<<q2.pt()<<endl;
0256   } else {
0257     event.getByToken(m_theJetTagsToken, jetTags);
0258     for (typename TCollection::const_iterator jet = jets->begin(); (jet != jets->end() && nJet < nMax); ++jet) {
0259       if (value_ == "second") {
0260         value = findCSV(jet, *jetTags);
0261       } else {
0262         value = 0.0;
0263       }
0264       sorted[nJet] = make_pair(value, nJet);
0265       ++nJet;
0266     }
0267     sort(sorted.begin(), sorted.end(), comparator);
0268     for (unsigned int i = 0; i < nMax; ++i) {
0269       jetRefs[i] = TRef(jets, sorted[i].second);
0270     }
0271     nJetRefs = nMax;
0272     b1 = jetRefs[3]->p4();
0273     b2 = jetRefs[2]->p4();
0274     q1 = jetRefs[1]->p4();
0275     q2 = jetRefs[0]->p4();
0276   }
0277 
0278   double mqq_bs = (q1 + q2).M();
0279   double deltaetaqq = std::abs(q1.Eta() - q2.Eta());
0280   double deltaetabb = std::abs(b1.Eta() - b2.Eta());
0281   double deltaphibb = std::abs(reco::deltaPhi(b1.Phi(), b2.Phi()));
0282   double ptsqq_bs = (q1 + q2).Pt();
0283   double ptsbb_bs = (b1 + b2).Pt();
0284   double signeta = q1.Eta() * q2.Eta();
0285 
0286   if ((mqq_bs > mqq_) && (deltaetaqq > detaqq_) && (deltaetabb < detabb_) && (deltaphibb < dphibb_) &&
0287       (ptsqq_bs > ptsqq_) && (ptsbb_bs > ptsbb_) && (signeta < seta_)) {
0288     accept = true;
0289     for (unsigned int i = 0; i < nJetRefs; ++i) {
0290       filterproduct.addObject(triggerType_, jetRefs[i]);
0291     }
0292   }
0293 
0294   return accept;
0295 }