TopJetCombination/plugins/TtSemiLepJetCombMaxSumPtWMass.cc

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/global/EDProducer.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004
0005 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
0006 #include "DataFormats/PatCandidates/interface/Jet.h"
0007
0008 class TtSemiLepJetCombMaxSumPtWMass : public edm::global::EDProducer<> {
0009 public:
0010   explicit TtSemiLepJetCombMaxSumPtWMass(const edm::ParameterSet&);
0011
0012 private:
0013   void produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& setup) const override;
0014
0015   bool isValid(const int& idx, const std::vector<pat::Jet>& jets) const {
0016     return (0 <= idx && idx < (int)jets.size());
0017   };
0018
0019   edm::EDGetTokenT<std::vector<pat::Jet>> jetsToken_;
0020   edm::EDGetTokenT<edm::View<reco::RecoCandidate>> lepsToken_;
0021   int maxNJets_;
0022   double wMass_;
0023   bool useBTagging_;
0024   std::string bTagAlgorithm_;
0025   double minBDiscBJets_;
0026   double maxBDiscLightJets_;
0027 };
0028
0029 TtSemiLepJetCombMaxSumPtWMass::TtSemiLepJetCombMaxSumPtWMass(const edm::ParameterSet& cfg)
0030     : jetsToken_(consumes<std::vector<pat::Jet>>(cfg.getParameter<edm::InputTag>("jets"))),
0031       lepsToken_(consumes<edm::View<reco::RecoCandidate>>(cfg.getParameter<edm::InputTag>("leps"))),
0032       maxNJets_(cfg.getParameter<int>("maxNJets")),
0033       wMass_(cfg.getParameter<double>("wMass")),
0034       useBTagging_(cfg.getParameter<bool>("useBTagging")),
0035       bTagAlgorithm_(cfg.getParameter<std::string>("bTagAlgorithm")),
0036       minBDiscBJets_(cfg.getParameter<double>("minBDiscBJets")),
0037       maxBDiscLightJets_(cfg.getParameter<double>("maxBDiscLightJets")) {
0038   if (maxNJets_ < 4 && maxNJets_ != -1)
0039     throw cms::Exception("WrongConfig") << "Parameter maxNJets can not be set to " << maxNJets_ << ". \n"
0040                                         << "It has to be larger than 4 or can be set to -1 to take all jets.";
0041
0042   produces<std::vector<std::vector<int>>>();
0043   produces<int>("NumberOfConsideredJets");
0044 }
0045
0046 void TtSemiLepJetCombMaxSumPtWMass::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& setup) const {
0047   auto pOut = std::make_unique<std::vector<std::vector<int>>>();
0048   auto pJetsConsidered = std::make_unique<int>(0);
0049
0050   std::vector<int> match;
0051   for (unsigned int i = 0; i < 4; ++i)
0052     match.push_back(-1);
0053
0054   // get jets
0055   const auto& jets = evt.get(jetsToken_);
0056
0057   // get leptons
0058   const auto& leps = evt.get(lepsToken_);
0059
0060   // skip events without lepton candidate or less than 4 jets or no MET
0061   if (leps.empty() || jets.size() < 4) {
0062     pOut->push_back(match);
0063     evt.put(std::move(pOut));
0064     *pJetsConsidered = jets.size();
0065     evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0066     return;
0067   }
0068
0069   unsigned maxNJets = maxNJets_;
0070   if (maxNJets_ == -1 || (int)jets.size() < maxNJets_)
0071     maxNJets = jets.size();
0072   *pJetsConsidered = maxNJets;
0073   evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0074
0075   std::vector<bool> isBJet;
0076   std::vector<bool> isLJet;
0077   int cntBJets = 0;
0078   if (useBTagging_) {
0079     for (unsigned int idx = 0; idx < maxNJets; ++idx) {
0080       isBJet.push_back((jets[idx].bDiscriminator(bTagAlgorithm_) > minBDiscBJets_));
0081       isLJet.push_back((jets[idx].bDiscriminator(bTagAlgorithm_) < maxBDiscLightJets_));
0082       if (jets[idx].bDiscriminator(bTagAlgorithm_) > minBDiscBJets_)
0083         cntBJets++;
0084     }
0085   }
0086
0087   // -----------------------------------------------------
0088   // associate those jets with maximum pt of the vectorial
0089   // sum to the hadronic decay chain
0090   // -----------------------------------------------------
0091   double maxPt = -1.;
0092   std::vector<int> maxPtIndices;
0093   maxPtIndices.push_back(-1);
0094   maxPtIndices.push_back(-1);
0095   maxPtIndices.push_back(-1);
0096   for (unsigned idx = 0; idx < maxNJets; ++idx) {
0097     if (useBTagging_ && (!isLJet[idx] || (cntBJets <= 2 && isBJet[idx])))
0098       continue;
0099     for (unsigned jdx = (idx + 1); jdx < maxNJets; ++jdx) {
0100       if (jdx == idx || (useBTagging_ && (!isLJet[jdx] || (cntBJets <= 2 && isBJet[jdx]) ||
0101                                           (cntBJets == 3 && isBJet[idx] && isBJet[jdx]))))
0102         continue;
0103       for (unsigned kdx = 0; kdx < maxNJets; ++kdx) {
0104         if (kdx == idx || kdx == jdx || (useBTagging_ && !isBJet[kdx]))
0105           continue;
0106         reco::Particle::LorentzVector sum = jets[idx].p4() + jets[jdx].p4() + jets[kdx].p4();
0107         if (maxPt < 0. || maxPt < sum.pt()) {
0108           maxPt = sum.pt();
0109           maxPtIndices.clear();
0110           maxPtIndices.push_back(idx);
0111           maxPtIndices.push_back(jdx);
0112           maxPtIndices.push_back(kdx);
0113         }
0114       }
0115     }
0116   }
0117
0118   // -----------------------------------------------------
0119   // associate those jets that get closest to the W mass
0120   // with their invariant mass to the W boson
0121   // -----------------------------------------------------
0122   double wDist = -1.;
0123   std::vector<int> closestToWMassIndices;
0124   closestToWMassIndices.push_back(-1);
0125   closestToWMassIndices.push_back(-1);
0126   if (isValid(maxPtIndices[0], jets) && isValid(maxPtIndices[1], jets) && isValid(maxPtIndices[2], jets)) {
0127     for (unsigned idx = 0; idx < maxPtIndices.size(); ++idx) {
0128       for (unsigned jdx = 0; jdx < maxPtIndices.size(); ++jdx) {
0129         if (jdx == idx || maxPtIndices[idx] > maxPtIndices[jdx] ||
0130             (useBTagging_ &&
0131              (!isLJet[maxPtIndices[idx]] || !isLJet[maxPtIndices[jdx]] ||
0132               (cntBJets <= 2 && isBJet[maxPtIndices[idx]]) || (cntBJets <= 2 && isBJet[maxPtIndices[jdx]]) ||
0133               (cntBJets == 3 && isBJet[maxPtIndices[idx]] && isBJet[maxPtIndices[jdx]]))))
0134           continue;
0135         reco::Particle::LorentzVector sum = jets[maxPtIndices[idx]].p4() + jets[maxPtIndices[jdx]].p4();
0136         if (wDist < 0. || wDist > fabs(sum.mass() - wMass_)) {
0137           wDist = fabs(sum.mass() - wMass_);
0138           closestToWMassIndices.clear();
0139           closestToWMassIndices.push_back(maxPtIndices[idx]);
0140           closestToWMassIndices.push_back(maxPtIndices[jdx]);
0141         }
0142       }
0143     }
0144   }
0145   int hadB = -1;
0146   for (unsigned idx = 0; idx < maxPtIndices.size(); ++idx) {
0147     // if this idx is not yet contained in the list of W mass candidates...
0148     if (std::find(closestToWMassIndices.begin(), closestToWMassIndices.end(), maxPtIndices[idx]) ==
0149         closestToWMassIndices.end()) {
0150       hadB = maxPtIndices[idx];
0151       break;  // there should be no other cadidates!
0152     }
0153   }
0154
0155   // -----------------------------------------------------
0156   // associate the remaining jet with maximum pt of the
0157   // vectorial sum with the leading lepton with the
0158   // leptonic decay chain
0159   // -----------------------------------------------------
0160   maxPt = -1.;
0161   int lepB = -1;
0162   for (unsigned idx = 0; idx < maxNJets; ++idx) {
0163     if (useBTagging_ && !isBJet[idx])
0164       continue;
0165     // make sure it's not used up already from the hadronic decay chain
0166     if (std::find(maxPtIndices.begin(), maxPtIndices.end(), idx) == maxPtIndices.end()) {
0167       reco::Particle::LorentzVector sum = jets[idx].p4() + leps[0].p4();
0168       if (maxPt < 0. || maxPt < sum.pt()) {
0169         maxPt = sum.pt();
0170         lepB = idx;
0171       }
0172     }
0173   }
0174
0175   match[TtSemiLepEvtPartons::LightQ] = closestToWMassIndices[0];
0176   match[TtSemiLepEvtPartons::LightQBar] = closestToWMassIndices[1];
0177   match[TtSemiLepEvtPartons::HadB] = hadB;
0178   match[TtSemiLepEvtPartons::LepB] = lepB;
0179
0180   pOut->push_back(match);
0181   evt.put(std::move(pOut));
0182 }
0183
0184 #include "FWCore/Framework/interface/MakerMacros.h"
0185 DEFINE_FWK_MODULE(TtSemiLepJetCombMaxSumPtWMass);