File indexing completed on 2023-10-25 10:05:26
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
0055 const auto& jets = evt.get(jetsToken_);
0056
0057
0058 const auto& leps = evt.get(lepsToken_);
0059
0060
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
0089
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
0120
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
0148 if (std::find(closestToWMassIndices.begin(), closestToWMassIndices.end(), maxPtIndices[idx]) ==
0149 closestToWMassIndices.end()) {
0150 hadB = maxPtIndices[idx];
0151 break;
0152 }
0153 }
0154
0155
0156
0157
0158
0159
0160 maxPt = -1.;
0161 int lepB = -1;
0162 for (unsigned idx = 0; idx < maxNJets; ++idx) {
0163 if (useBTagging_ && !isBJet[idx])
0164 continue;
0165
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);