TopJetCombination/plugins/TtSemiLepJetCombGeom.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/Math/interface/LorentzVector.h"
0007 #include "DataFormats/PatCandidates/interface/Jet.h"
0008 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0009
0010 #include "Math/VectorUtil.h"
0011
0012 class TtSemiLepJetCombGeom : public edm::global::EDProducer<> {
0013 public:
0014   explicit TtSemiLepJetCombGeom(const edm::ParameterSet&);
0015
0016 private:
0017   void produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& setup) const override;
0018
0019   bool isValid(const int idx, const std::vector<pat::Jet>& jets) const { return (0 <= idx && idx < (int)jets.size()); };
0020   double distance(const math::XYZTLorentzVector&, const math::XYZTLorentzVector&) const;
0021
0022   edm::EDGetTokenT<std::vector<pat::Jet>> jetsToken_;
0023   edm::EDGetTokenT<edm::View<reco::RecoCandidate>> lepsToken_;
0024   int maxNJets_;
0025   bool useDeltaR_;
0026   bool useBTagging_;
0027   std::string bTagAlgorithm_;
0028   double minBDiscBJets_;
0029   double maxBDiscLightJets_;
0030 };
0031
0032 TtSemiLepJetCombGeom::TtSemiLepJetCombGeom(const edm::ParameterSet& cfg)
0033     : jetsToken_(consumes<std::vector<pat::Jet>>(cfg.getParameter<edm::InputTag>("jets"))),
0034       lepsToken_(consumes<edm::View<reco::RecoCandidate>>(cfg.getParameter<edm::InputTag>("leps"))),
0035       maxNJets_(cfg.getParameter<int>("maxNJets")),
0036       useDeltaR_(cfg.getParameter<bool>("useDeltaR")),
0037       useBTagging_(cfg.getParameter<bool>("useBTagging")),
0038       bTagAlgorithm_(cfg.getParameter<std::string>("bTagAlgorithm")),
0039       minBDiscBJets_(cfg.getParameter<double>("minBDiscBJets")),
0040       maxBDiscLightJets_(cfg.getParameter<double>("maxBDiscLightJets")) {
0041   if (maxNJets_ < 4 && maxNJets_ != -1)
0042     throw cms::Exception("WrongConfig") << "Parameter maxNJets can not be set to " << maxNJets_ << ". \n"
0043                                         << "It has to be larger than 4 or can be set to -1 to take all jets.";
0044
0045   produces<std::vector<std::vector<int>>>();
0046   produces<int>("NumberOfConsideredJets");
0047 }
0048
0049 void TtSemiLepJetCombGeom::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& setup) const {
0050   auto pOut = std::make_unique<std::vector<std::vector<int>>>();
0051   auto pJetsConsidered = std::make_unique<int>(0);
0052
0053   std::vector<int> match;
0054   for (unsigned int i = 0; i < 4; ++i)
0055     match.push_back(-1);
0056
0057   // get jets
0058   const auto& jets = evt.get(jetsToken_);
0059
0060   // get leptons
0061   const auto& leps = evt.get(lepsToken_);
0062
0063   // skip events without lepton candidate or less than 4 jets
0064   if (leps.empty() || jets.size() < 4) {
0065     pOut->push_back(match);
0066     evt.put(std::move(pOut));
0067     *pJetsConsidered = jets.size();
0068     evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0069     return;
0070   }
0071
0072   unsigned maxNJets = maxNJets_;
0073   if (maxNJets_ == -1 || (int)jets.size() < maxNJets_)
0074     maxNJets = jets.size();
0075   *pJetsConsidered = maxNJets;
0076   evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0077
0078   std::vector<bool> isBJet;
0079   std::vector<bool> isLJet;
0080   int cntBJets = 0;
0081   if (useBTagging_) {
0082     for (unsigned int idx = 0; idx < maxNJets; ++idx) {
0083       isBJet.push_back((jets[idx].bDiscriminator(bTagAlgorithm_) > minBDiscBJets_));
0084       isLJet.push_back((jets[idx].bDiscriminator(bTagAlgorithm_) < maxBDiscLightJets_));
0085       if (jets[idx].bDiscriminator(bTagAlgorithm_) > minBDiscBJets_)
0086         cntBJets++;
0087     }
0088   }
0089
0090   // -----------------------------------------------------
0091   // associate those two jets to the hadronic W boson that
0092   // have the smallest distance to each other
0093   // -----------------------------------------------------
0094   double minDist = -1.;
0095   int lightQ = -1;
0096   int lightQBar = -1;
0097   for (unsigned int idx = 0; idx < maxNJets; ++idx) {
0098     if (useBTagging_ && (!isLJet[idx] || (cntBJets <= 2 && isBJet[idx])))
0099       continue;
0100     for (unsigned int jdx = (idx + 1); jdx < maxNJets; ++jdx) {
0101       if (useBTagging_ &&
0102           (!isLJet[jdx] || (cntBJets <= 2 && isBJet[jdx]) || (cntBJets == 3 && isBJet[idx] && isBJet[jdx])))
0103         continue;
0104       double dist = distance(jets[idx].p4(), jets[jdx].p4());
0105       if (minDist < 0. || dist < minDist) {
0106         minDist = dist;
0107         lightQ = idx;
0108         lightQBar = jdx;
0109       }
0110     }
0111   }
0112
0113   reco::Particle::LorentzVector wHad;
0114   if (isValid(lightQ, jets) && isValid(lightQBar, jets))
0115     wHad = jets[lightQ].p4() + jets[lightQBar].p4();
0116
0117   // -----------------------------------------------------
0118   // associate to the hadronic b quark the remaining jet
0119   // that has the smallest distance to the hadronic W
0120   // -----------------------------------------------------
0121   minDist = -1.;
0122   int hadB = -1;
0123   if (isValid(lightQ, jets) && isValid(lightQBar, jets)) {
0124     for (unsigned int idx = 0; idx < maxNJets; ++idx) {
0125       if (useBTagging_ && !isBJet[idx])
0126         continue;
0127       // make sure it's not used up already from the hadronic W
0128       if ((int)idx != lightQ && (int)idx != lightQBar) {
0129         double dist = distance(jets[idx].p4(), wHad);
0130         if (minDist < 0. || dist < minDist) {
0131           minDist = dist;
0132           hadB = idx;
0133         }
0134       }
0135     }
0136   }
0137
0138   // -----------------------------------------------------
0139   // associate to the leptonic b quark the remaining jet
0140   // that has the smallest distance to the leading lepton
0141   // -----------------------------------------------------
0142   minDist = -1.;
0143   int lepB = -1;
0144   for (unsigned int idx = 0; idx < maxNJets; ++idx) {
0145     if (useBTagging_ && !isBJet[idx])
0146       continue;
0147     // make sure it's not used up already from the hadronic decay chain
0148     if ((int)idx != lightQ && (int)idx != lightQBar && (int)idx != hadB) {
0149       double dist = distance(jets[idx].p4(), leps[0].p4());
0150       if (minDist < 0. || dist < minDist) {
0151         minDist = dist;
0152         lepB = idx;
0153       }
0154     }
0155   }
0156
0157   match[TtSemiLepEvtPartons::LightQ] = lightQ;
0158   match[TtSemiLepEvtPartons::LightQBar] = lightQBar;
0159   match[TtSemiLepEvtPartons::HadB] = hadB;
0160   match[TtSemiLepEvtPartons::LepB] = lepB;
0161
0162   pOut->push_back(match);
0163   evt.put(std::move(pOut));
0164 }
0165
0166 double TtSemiLepJetCombGeom::distance(const math::XYZTLorentzVector& v1, const math::XYZTLorentzVector& v2) const {
0167   // calculate the distance between two lorentz vectors
0168   // using DeltaR or DeltaTheta
0169   if (useDeltaR_)
0170     return ROOT::Math::VectorUtil::DeltaR(v1, v2);
0171   return fabs(v1.theta() - v2.theta());
0172 }
0173
0174 #include "FWCore/Framework/interface/MakerMacros.h"
0175 DEFINE_FWK_MODULE(TtSemiLepJetCombGeom);