File indexing completed on 2023-03-17 11:26:11
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
0058 const auto& jets = evt.get(jetsToken_);
0059
0060
0061 const auto& leps = evt.get(lepsToken_);
0062
0063
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
0092
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
0119
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
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
0140
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
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
0168
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);