TopJetCombination/plugins/TtSemiLepJetCombWMassMaxSumPt.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 TtSemiLepJetCombWMassMaxSumPt : public edm::global::EDProducer<> {
0009 public:
0010   explicit TtSemiLepJetCombWMassMaxSumPt(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 TtSemiLepJetCombWMassMaxSumPt::TtSemiLepJetCombWMassMaxSumPt(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 TtSemiLepJetCombWMassMaxSumPt::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 that get closest to the W mass
0089   // with their invariant mass to the hadronic W boson
0090   // -----------------------------------------------------
0091   double wDist = -1.;
0092   std::vector<int> closestToWMassIndices;
0093   closestToWMassIndices.push_back(-1);
0094   closestToWMassIndices.push_back(-1);
0095   for (unsigned idx = 0; idx < maxNJets; ++idx) {
0096     if (useBTagging_ && (!isLJet[idx] || (cntBJets <= 2 && isBJet[idx])))
0097       continue;
0098     for (unsigned jdx = (idx + 1); jdx < maxNJets; ++jdx) {
0099       if (useBTagging_ &&
0100           (!isLJet[jdx] || (cntBJets <= 2 && isBJet[jdx]) || (cntBJets == 3 && isBJet[idx] && isBJet[jdx])))
0101         continue;
0102       reco::Particle::LorentzVector sum = jets[idx].p4() + jets[jdx].p4();
0103       if (wDist < 0. || wDist > fabs(sum.mass() - wMass_)) {
0104         wDist = fabs(sum.mass() - wMass_);
0105         closestToWMassIndices.clear();
0106         closestToWMassIndices.push_back(idx);
0107         closestToWMassIndices.push_back(jdx);
0108       }
0109     }
0110   }
0111
0112   // -----------------------------------------------------
0113   // associate those jets with maximum pt of the vectorial
0114   // sum to the hadronic decay chain
0115   // -----------------------------------------------------
0116   double maxPt = -1.;
0117   int hadB = -1;
0118   if (isValid(closestToWMassIndices[0], jets) && isValid(closestToWMassIndices[1], jets)) {
0119     for (unsigned idx = 0; idx < maxNJets; ++idx) {
0120       if (useBTagging_ && !isBJet[idx])
0121         continue;
0122       // make sure it's not used up already from the hadronic W
0123       if ((int)idx != closestToWMassIndices[0] && (int)idx != closestToWMassIndices[1]) {
0124         reco::Particle::LorentzVector sum =
0125             jets[closestToWMassIndices[0]].p4() + jets[closestToWMassIndices[1]].p4() + jets[idx].p4();
0126         if (maxPt < 0. || maxPt < sum.pt()) {
0127           maxPt = sum.pt();
0128           hadB = idx;
0129         }
0130       }
0131     }
0132   }
0133
0134   // -----------------------------------------------------
0135   // associate the remaining jet with maximum pt of the
0136   // vectorial sum with the leading lepton with the
0137   // leptonic b quark
0138   // -----------------------------------------------------
0139   maxPt = -1.;
0140   int lepB = -1;
0141   for (unsigned idx = 0; idx < maxNJets; ++idx) {
0142     if (useBTagging_ && !isBJet[idx])
0143       continue;
0144     // make sure it's not used up already from the hadronic decay chain
0145     if ((int)idx != closestToWMassIndices[0] && (int)idx != closestToWMassIndices[1] && (int)idx != hadB) {
0146       reco::Particle::LorentzVector sum = jets[idx].p4() + leps[0].p4();
0147       if (maxPt < 0. || maxPt < sum.pt()) {
0148         maxPt = sum.pt();
0149         lepB = idx;
0150       }
0151     }
0152   }
0153
0154   match[TtSemiLepEvtPartons::LightQ] = closestToWMassIndices[0];
0155   match[TtSemiLepEvtPartons::LightQBar] = closestToWMassIndices[1];
0156   match[TtSemiLepEvtPartons::HadB] = hadB;
0157   match[TtSemiLepEvtPartons::LepB] = lepB;
0158
0159   pOut->push_back(match);
0160   evt.put(std::move(pOut));
0161 }
0162
0163 #include "FWCore/Framework/interface/MakerMacros.h"
0164 DEFINE_FWK_MODULE(TtSemiLepJetCombWMassMaxSumPt);