TopJetCombination/plugins/TtSemiLepJetCombWMassDeltaTopMass.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/Electron.h"
0007 #include "DataFormats/PatCandidates/interface/Muon.h"
0008 #include "DataFormats/PatCandidates/interface/Jet.h"
0009 #include "DataFormats/PatCandidates/interface/MET.h"
0010
0011 #include "TopQuarkAnalysis/TopTools/interface/MEzCalculator.h"
0012
0013 class TtSemiLepJetCombWMassDeltaTopMass : public edm::global::EDProducer<> {
0014 public:
0015   explicit TtSemiLepJetCombWMassDeltaTopMass(const edm::ParameterSet&);
0016
0017 private:
0018   void produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& setup) const override;
0019
0020   bool isValid(const int& idx, const std::vector<pat::Jet>& jets) const {
0021     return (0 <= idx && idx < (int)jets.size());
0022   };
0023
0024   edm::EDGetTokenT<std::vector<pat::Jet>> jetsToken_;
0025   edm::EDGetTokenT<edm::View<reco::RecoCandidate>> lepsToken_;
0026   edm::EDGetTokenT<std::vector<pat::MET>> metsToken_;
0027   int maxNJets_;
0028   double wMass_;
0029   bool useBTagging_;
0030   std::string bTagAlgorithm_;
0031   double minBDiscBJets_;
0032   double maxBDiscLightJets_;
0033   int neutrinoSolutionType_;
0034 };
0035
0036 TtSemiLepJetCombWMassDeltaTopMass::TtSemiLepJetCombWMassDeltaTopMass(const edm::ParameterSet& cfg)
0037     : jetsToken_(consumes<std::vector<pat::Jet>>(cfg.getParameter<edm::InputTag>("jets"))),
0038       lepsToken_(consumes<edm::View<reco::RecoCandidate>>(cfg.getParameter<edm::InputTag>("leps"))),
0039       metsToken_(consumes<std::vector<pat::MET>>(cfg.getParameter<edm::InputTag>("mets"))),
0040       maxNJets_(cfg.getParameter<int>("maxNJets")),
0041       wMass_(cfg.getParameter<double>("wMass")),
0042       useBTagging_(cfg.getParameter<bool>("useBTagging")),
0043       bTagAlgorithm_(cfg.getParameter<std::string>("bTagAlgorithm")),
0044       minBDiscBJets_(cfg.getParameter<double>("minBDiscBJets")),
0045       maxBDiscLightJets_(cfg.getParameter<double>("maxBDiscLightJets")),
0046       neutrinoSolutionType_(cfg.getParameter<int>("neutrinoSolutionType")) {
0047   if (maxNJets_ < 4 && maxNJets_ != -1)
0048     throw cms::Exception("WrongConfig") << "Parameter maxNJets can not be set to " << maxNJets_ << ". \n"
0049                                         << "It has to be larger than 4 or can be set to -1 to take all jets.";
0050
0051   produces<std::vector<std::vector<int>>>();
0052   produces<int>("NumberOfConsideredJets");
0053 }
0054
0055 void TtSemiLepJetCombWMassDeltaTopMass::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& setup) const {
0056   auto pOut = std::make_unique<std::vector<std::vector<int>>>();
0057   auto pJetsConsidered = std::make_unique<int>(0);
0058
0059   std::vector<int> match;
0060   for (unsigned int i = 0; i < 4; ++i)
0061     match.push_back(-1);
0062
0063   // get jets
0064   const auto& jets = evt.get(jetsToken_);
0065
0066   // get leptons
0067   const auto& leps = evt.get(lepsToken_);
0068
0069   // get MET
0070   const auto& mets = evt.get(metsToken_);
0071
0072   // skip events without lepton candidate or less than 4 jets or no MET
0073   if (leps.empty() || jets.size() < 4 || mets.empty()) {
0074     pOut->push_back(match);
0075     evt.put(std::move(pOut));
0076     *pJetsConsidered = jets.size();
0077     evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0078     return;
0079   }
0080
0081   unsigned maxNJets = maxNJets_;
0082   if (maxNJets_ == -1 || (int)jets.size() < maxNJets_)
0083     maxNJets = jets.size();
0084   *pJetsConsidered = maxNJets;
0085   evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0086
0087   std::vector<bool> isBJet;
0088   std::vector<bool> isLJet;
0089   int cntBJets = 0;
0090   if (useBTagging_) {
0091     for (unsigned int idx = 0; idx < maxNJets; ++idx) {
0092       isBJet.push_back((jets[idx].bDiscriminator(bTagAlgorithm_) > minBDiscBJets_));
0093       isLJet.push_back((jets[idx].bDiscriminator(bTagAlgorithm_) < maxBDiscLightJets_));
0094       if (jets[idx].bDiscriminator(bTagAlgorithm_) > minBDiscBJets_)
0095         cntBJets++;
0096     }
0097   }
0098
0099   // -----------------------------------------------------
0100   // associate those jets that get closest to the W mass
0101   // with their invariant mass to the hadronic W boson
0102   // -----------------------------------------------------
0103   double wDist = -1.;
0104   std::vector<int> closestToWMassIndices;
0105   closestToWMassIndices.push_back(-1);
0106   closestToWMassIndices.push_back(-1);
0107   for (unsigned idx = 0; idx < maxNJets; ++idx) {
0108     if (useBTagging_ && (!isLJet[idx] || (cntBJets <= 2 && isBJet[idx])))
0109       continue;
0110     for (unsigned jdx = (idx + 1); jdx < maxNJets; ++jdx) {
0111       if (useBTagging_ &&
0112           (!isLJet[jdx] || (cntBJets <= 2 && isBJet[jdx]) || (cntBJets == 3 && isBJet[idx] && isBJet[jdx])))
0113         continue;
0114       reco::Particle::LorentzVector sum = jets[idx].p4() + jets[jdx].p4();
0115       if (wDist < 0. || wDist > fabs(sum.mass() - wMass_)) {
0116         wDist = fabs(sum.mass() - wMass_);
0117         closestToWMassIndices.clear();
0118         closestToWMassIndices.push_back(idx);
0119         closestToWMassIndices.push_back(jdx);
0120       }
0121     }
0122   }
0123
0124   // -----------------------------------------------------
0125   // build a leptonic W boson from the lepton and the MET
0126   // -----------------------------------------------------
0127   double neutrino_pz = 0.;
0128   if (neutrinoSolutionType_ != -1) {
0129     MEzCalculator mez;
0130     mez.SetMET(*mets.begin());
0131     if (dynamic_cast<const reco::Muon*>(&(leps.front())))
0132       mez.SetLepton(leps[0], true);
0133     else if (dynamic_cast<const reco::GsfElectron*>(&(leps.front())))
0134       mez.SetLepton(leps[0], false);
0135     else
0136       throw cms::Exception("UnimplementedFeature")
0137           << "Type of lepton given together with MET for solving neutrino kinematics is neither muon nor electron.\n";
0138     neutrino_pz = mez.Calculate(neutrinoSolutionType_);
0139   }
0140   const math::XYZTLorentzVector neutrino(mets.begin()->px(),
0141                                          mets.begin()->py(),
0142                                          neutrino_pz,
0143                                          sqrt(mets.begin()->px() * mets.begin()->px() +
0144                                               mets.begin()->py() * mets.begin()->py() + neutrino_pz * neutrino_pz));
0145   const reco::Particle::LorentzVector lepW = neutrino + leps.front().p4();
0146
0147   // -----------------------------------------------------
0148   // associate those jets to the hadronic and the leptonic
0149   // b quark that minimize the difference between
0150   // hadronic and leptonic top mass
0151   // -----------------------------------------------------
0152   double deltaTop = -1.;
0153   int hadB = -1;
0154   int lepB = -1;
0155   if (isValid(closestToWMassIndices[0], jets) && isValid(closestToWMassIndices[1], jets)) {
0156     const reco::Particle::LorentzVector hadW =
0157         jets[closestToWMassIndices[0]].p4() + jets[closestToWMassIndices[1]].p4();
0158     // find hadronic b candidate
0159     for (unsigned idx = 0; idx < maxNJets; ++idx) {
0160       if (useBTagging_ && !isBJet[idx])
0161         continue;
0162       // make sure it's not used up already from the hadronic W
0163       if ((int)idx != closestToWMassIndices[0] && (int)idx != closestToWMassIndices[1]) {
0164         reco::Particle::LorentzVector hadTop = hadW + jets[idx].p4();
0165         // find leptonic b candidate
0166         for (unsigned jdx = 0; jdx < maxNJets; ++jdx) {
0167           if (useBTagging_ && !isBJet[jdx])
0168             continue;
0169           // make sure it's not used up already from the hadronic branch
0170           if ((int)jdx != closestToWMassIndices[0] && (int)jdx != closestToWMassIndices[1] && jdx != idx) {
0171             reco::Particle::LorentzVector lepTop = lepW + jets[jdx].p4();
0172             if (deltaTop < 0. || deltaTop > fabs(hadTop.mass() - lepTop.mass())) {
0173               deltaTop = fabs(hadTop.mass() - lepTop.mass());
0174               hadB = idx;
0175               lepB = jdx;
0176             }
0177           }
0178         }
0179       }
0180     }
0181   }
0182
0183   match[TtSemiLepEvtPartons::LightQ] = closestToWMassIndices[0];
0184   match[TtSemiLepEvtPartons::LightQBar] = closestToWMassIndices[1];
0185   match[TtSemiLepEvtPartons::HadB] = hadB;
0186   match[TtSemiLepEvtPartons::LepB] = lepB;
0187
0188   pOut->push_back(match);
0189   evt.put(std::move(pOut));
0190 }
0191
0192 #include "FWCore/Framework/interface/MakerMacros.h"
0193 DEFINE_FWK_MODULE(TtSemiLepJetCombWMassDeltaTopMass);