File indexing completed on 2021-02-14 13:32:44
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021
0022
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/EDProducer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/StreamID.h"
0031
0032 #include "DataFormats/Math/interface/deltaR.h"
0033
0034 #include "DataFormats/PatCandidates/interface/Muon.h"
0035 #include "DataFormats/PatCandidates/interface/Electron.h"
0036 #include "DataFormats/PatCandidates/interface/Jet.h"
0037
0038 #include "DataFormats/VertexReco/interface/Vertex.h"
0039 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0040
0041 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0042 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0043 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
0044
0045 #include <vector>
0046
0047
0048 #include "FWCore/Framework/interface/Frameworkfwd.h"
0049 #include "FWCore/Framework/interface/global/EDProducer.h"
0050
0051 #include "FWCore/Framework/interface/Event.h"
0052 #include "FWCore/Framework/interface/MakerMacros.h"
0053
0054 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0055 #include "FWCore/Utilities/interface/StreamID.h"
0056
0057 #include "DataFormats/PatCandidates/interface/Muon.h"
0058 #include "DataFormats/PatCandidates/interface/Jet.h"
0059 #include "DataFormats/PatCandidates/interface/Electron.h"
0060 #include "DataFormats/VertexReco/interface/Vertex.h"
0061
0062 #include "DataFormats/Common/interface/View.h"
0063 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0064 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0065
0066 #include "PhysicsTools/NanoAOD/interface/MatchingUtils.h"
0067
0068
0069
0070
0071
0072 template <typename T>
0073 class BJetEnergyRegressionVarProducer : public edm::global::EDProducer<> {
0074 public:
0075 explicit BJetEnergyRegressionVarProducer(const edm::ParameterSet& iConfig)
0076 : srcJet_(consumes<edm::View<pat::Jet>>(iConfig.getParameter<edm::InputTag>("src"))),
0077 srcVtx_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("pvsrc"))),
0078 srcSV_(consumes<edm::View<reco::VertexCompositePtrCandidate>>(iConfig.getParameter<edm::InputTag>("svsrc"))),
0079 srcGP_(consumes<std::vector<reco::GenParticle>>(iConfig.getParameter<edm::InputTag>("gpsrc"))) {
0080
0081 produces<edm::ValueMap<float>>("leptonPtRel");
0082 produces<edm::ValueMap<float>>("leptonPtRatio");
0083 produces<edm::ValueMap<float>>("leptonPtRelInv");
0084 produces<edm::ValueMap<float>>("leptonPtRelv0");
0085 produces<edm::ValueMap<float>>("leptonPtRatiov0");
0086 produces<edm::ValueMap<float>>("leptonPtRelInvv0");
0087 produces<edm::ValueMap<float>>("leptonPt");
0088 produces<edm::ValueMap<int>>("leptonPdgId");
0089 produces<edm::ValueMap<float>>("leptonDeltaR");
0090 produces<edm::ValueMap<float>>("leadTrackPt");
0091 produces<edm::ValueMap<float>>("vtxPt");
0092 produces<edm::ValueMap<float>>("vtxMass");
0093 produces<edm::ValueMap<float>>("vtx3dL");
0094 produces<edm::ValueMap<float>>("vtx3deL");
0095 produces<edm::ValueMap<int>>("vtxNtrk");
0096 produces<edm::ValueMap<float>>("ptD");
0097 produces<edm::ValueMap<float>>("genPtwNu");
0098 }
0099 ~BJetEnergyRegressionVarProducer() override{};
0100
0101 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0102
0103 private:
0104 void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0105
0106 std::tuple<float, float, float> calculatePtRatioRel(edm::Ptr<reco::Candidate> lep, edm::Ptr<pat::Jet> jet) const;
0107 std::tuple<float, float, float> calculatePtRatioRelSimple(edm::Ptr<reco::Candidate> lep,
0108 edm::Ptr<pat::Jet> jet) const;
0109
0110
0111
0112 edm::EDGetTokenT<edm::View<pat::Jet>> srcJet_;
0113 edm::EDGetTokenT<std::vector<reco::Vertex>> srcVtx_;
0114 edm::EDGetTokenT<edm::View<reco::VertexCompositePtrCandidate>> srcSV_;
0115 edm::EDGetTokenT<std::vector<reco::GenParticle>> srcGP_;
0116 };
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131 template <typename T>
0132 void BJetEnergyRegressionVarProducer<T>::produce(edm::StreamID streamID,
0133 edm::Event& iEvent,
0134 const edm::EventSetup& iSetup) const {
0135 edm::Handle<edm::View<pat::Jet>> srcJet;
0136 iEvent.getByToken(srcJet_, srcJet);
0137 edm::Handle<std::vector<reco::Vertex>> srcVtx;
0138 iEvent.getByToken(srcVtx_, srcVtx);
0139 edm::Handle<edm::View<reco::VertexCompositePtrCandidate>> srcSV;
0140 iEvent.getByToken(srcSV_, srcSV);
0141 edm::Handle<std::vector<reco::GenParticle>> srcGP;
0142 iEvent.getByToken(srcGP_, srcGP);
0143
0144 unsigned int nJet = srcJet->size();
0145
0146
0147 std::vector<float> leptonPtRel(nJet, 0);
0148 std::vector<float> leptonPtRatio(nJet, 0);
0149 std::vector<float> leptonPtRelInv(nJet, 0);
0150 std::vector<float> leptonPtRel_v0(nJet, 0);
0151 std::vector<float> leptonPtRatio_v0(nJet, 0);
0152 std::vector<float> leptonPtRelInv_v0(nJet, 0);
0153 std::vector<int> leptonPdgId(nJet, 0);
0154 std::vector<float> leptonPt(nJet, 0);
0155 std::vector<float> leptonDeltaR(nJet, 0);
0156 std::vector<float> leadTrackPt(nJet, 0);
0157 std::vector<float> vtxPt(nJet, 0);
0158 std::vector<float> vtxMass(nJet, 0);
0159 std::vector<float> vtx3dL(nJet, 0);
0160 std::vector<float> vtx3deL(nJet, 0);
0161 std::vector<int> vtxNtrk(nJet, 0);
0162 std::vector<float> ptD(nJet, 0);
0163 std::vector<float> genPtwNu(nJet, 0);
0164
0165 const auto& pv = (*srcVtx)[0];
0166 for (unsigned int ij = 0; ij < nJet; ij++) {
0167 auto jet = srcJet->ptrAt(ij);
0168
0169 if (jet->genJet() != nullptr) {
0170 auto genp4 = jet->genJet()->p4();
0171 auto gep4wNu = genp4;
0172 for (const auto& gp : *srcGP) {
0173 if ((abs(gp.pdgId()) == 12 || abs(gp.pdgId()) == 14 || abs(gp.pdgId()) == 16) && gp.status() == 1) {
0174 if (reco::deltaR(genp4, gp.p4()) < 0.4) {
0175
0176 gep4wNu = gep4wNu + gp.p4();
0177
0178 }
0179 }
0180 }
0181
0182 genPtwNu[ij] = gep4wNu.pt();
0183 }
0184
0185 float ptMax = 0;
0186 float sumWeight = 0;
0187 float sumPt = 0;
0188
0189 for (const auto& d : jet->daughterPtrVector()) {
0190 sumWeight += (d->pt()) * (d->pt());
0191 sumPt += d->pt();
0192 if (d->pt() > ptMax)
0193 ptMax = d->pt();
0194 }
0195 leadTrackPt[ij] = ptMax;
0196 ptD[ij] = (sumWeight > 0 ? sqrt(sumWeight) / sumPt : 0);
0197
0198
0199 float maxLepPt = 0;
0200 leptonPtRel[ij] = 0;
0201
0202 for (const auto& d : jet->daughterPtrVector()) {
0203 if (abs(d->pdgId()) == 11 || abs(d->pdgId()) == 13) {
0204 if (d->pt() < maxLepPt)
0205 continue;
0206 auto res = calculatePtRatioRel(d, jet);
0207 leptonPtRatio[ij] = std::get<0>(res);
0208 leptonPtRel[ij] = std::get<1>(res);
0209 leptonPtRelInv[ij] = std::get<2>(res);
0210 auto res2 = calculatePtRatioRelSimple(d, jet);
0211 leptonPtRatio_v0[ij] = std::get<0>(res2);
0212 leptonPtRel_v0[ij] = std::get<1>(res2);
0213 leptonPtRelInv_v0[ij] = std::get<2>(res2);
0214 leptonPdgId[ij] = d->pdgId();
0215 leptonDeltaR[ij] = reco::deltaR(jet->p4(), d->p4());
0216 leptonPt[ij] = d->pt();
0217 maxLepPt = d->pt();
0218 }
0219 }
0220
0221
0222 VertexDistance3D vdist;
0223 float maxFoundSignificance = 0;
0224
0225 vtxPt[ij] = 0;
0226 vtxMass[ij] = 0;
0227 vtx3dL[ij] = 0;
0228 vtx3deL[ij] = 0;
0229 vtxNtrk[ij] = 0;
0230
0231 for (const auto& sv : *srcSV) {
0232 GlobalVector flightDir(sv.vertex().x() - pv.x(), sv.vertex().y() - pv.y(), sv.vertex().z() - pv.z());
0233 GlobalVector jetDir(jet->px(), jet->py(), jet->pz());
0234 if (reco::deltaR2(flightDir, jetDir) < 0.09) {
0235 Measurement1D dl = vdist.distance(
0236 pv, VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0237 if (dl.significance() > maxFoundSignificance) {
0238 maxFoundSignificance = dl.significance();
0239 vtxPt[ij] = sv.pt();
0240 vtxMass[ij] = sv.p4().M();
0241 vtx3dL[ij] = dl.value();
0242 vtx3deL[ij] = dl.error();
0243 vtxNtrk[ij] = sv.numberOfSourceCandidatePtrs();
0244 }
0245 }
0246 }
0247 }
0248
0249 std::unique_ptr<edm::ValueMap<float>> leptonPtRelV(new edm::ValueMap<float>());
0250 edm::ValueMap<float>::Filler fillerRel(*leptonPtRelV);
0251 fillerRel.insert(srcJet, leptonPtRel.begin(), leptonPtRel.end());
0252 fillerRel.fill();
0253 iEvent.put(std::move(leptonPtRelV), "leptonPtRel");
0254
0255 std::unique_ptr<edm::ValueMap<float>> leptonPtRatioV(new edm::ValueMap<float>());
0256 edm::ValueMap<float>::Filler fillerRatio(*leptonPtRatioV);
0257 fillerRatio.insert(srcJet, leptonPtRatio.begin(), leptonPtRatio.end());
0258 fillerRatio.fill();
0259 iEvent.put(std::move(leptonPtRatioV), "leptonPtRatio");
0260
0261 std::unique_ptr<edm::ValueMap<float>> leptonPtRelInvV(new edm::ValueMap<float>());
0262 edm::ValueMap<float>::Filler fillerRelInv(*leptonPtRelInvV);
0263 fillerRelInv.insert(srcJet, leptonPtRelInv.begin(), leptonPtRelInv.end());
0264 fillerRelInv.fill();
0265 iEvent.put(std::move(leptonPtRelInvV), "leptonPtRelInv");
0266
0267 std::unique_ptr<edm::ValueMap<float>> leptonPtRelV_v0(new edm::ValueMap<float>());
0268 edm::ValueMap<float>::Filler fillerRel_v0(*leptonPtRelV_v0);
0269 fillerRel_v0.insert(srcJet, leptonPtRel_v0.begin(), leptonPtRel_v0.end());
0270 fillerRel_v0.fill();
0271 iEvent.put(std::move(leptonPtRelV_v0), "leptonPtRelv0");
0272
0273 std::unique_ptr<edm::ValueMap<float>> leptonPtRatioV_v0(new edm::ValueMap<float>());
0274 edm::ValueMap<float>::Filler fillerRatio_v0(*leptonPtRatioV_v0);
0275 fillerRatio_v0.insert(srcJet, leptonPtRatio_v0.begin(), leptonPtRatio_v0.end());
0276 fillerRatio_v0.fill();
0277 iEvent.put(std::move(leptonPtRatioV_v0), "leptonPtRatiov0");
0278
0279 std::unique_ptr<edm::ValueMap<float>> leptonPtRelInvV_v0(new edm::ValueMap<float>());
0280 edm::ValueMap<float>::Filler fillerRelInv_v0(*leptonPtRelInvV_v0);
0281 fillerRelInv_v0.insert(srcJet, leptonPtRelInv_v0.begin(), leptonPtRelInv_v0.end());
0282 fillerRelInv_v0.fill();
0283 iEvent.put(std::move(leptonPtRelInvV_v0), "leptonPtRelInvv0");
0284
0285 std::unique_ptr<edm::ValueMap<float>> leptonPtV(new edm::ValueMap<float>());
0286 edm::ValueMap<float>::Filler fillerLpt(*leptonPtV);
0287 fillerLpt.insert(srcJet, leptonPt.begin(), leptonPt.end());
0288 fillerLpt.fill();
0289 iEvent.put(std::move(leptonPtV), "leptonPt");
0290
0291 std::unique_ptr<edm::ValueMap<float>> leptonDeltaRV(new edm::ValueMap<float>());
0292 edm::ValueMap<float>::Filler fillerLdR(*leptonDeltaRV);
0293 fillerLdR.insert(srcJet, leptonDeltaR.begin(), leptonDeltaR.end());
0294 fillerLdR.fill();
0295 iEvent.put(std::move(leptonDeltaRV), "leptonDeltaR");
0296
0297 std::unique_ptr<edm::ValueMap<int>> leptonPtPdgIdV(new edm::ValueMap<int>());
0298 edm::ValueMap<int>::Filler fillerId(*leptonPtPdgIdV);
0299 fillerId.insert(srcJet, leptonPdgId.begin(), leptonPdgId.end());
0300 fillerId.fill();
0301 iEvent.put(std::move(leptonPtPdgIdV), "leptonPdgId");
0302
0303 std::unique_ptr<edm::ValueMap<float>> leadTrackPtV(new edm::ValueMap<float>());
0304 edm::ValueMap<float>::Filler fillerLT(*leadTrackPtV);
0305 fillerLT.insert(srcJet, leadTrackPt.begin(), leadTrackPt.end());
0306 fillerLT.fill();
0307 iEvent.put(std::move(leadTrackPtV), "leadTrackPt");
0308
0309 std::unique_ptr<edm::ValueMap<float>> vtxPtV(new edm::ValueMap<float>());
0310 edm::ValueMap<float>::Filler fillerVtxPt(*vtxPtV);
0311 fillerVtxPt.insert(srcJet, vtxPt.begin(), vtxPt.end());
0312 fillerVtxPt.fill();
0313 iEvent.put(std::move(vtxPtV), "vtxPt");
0314
0315 std::unique_ptr<edm::ValueMap<float>> vtxMassV(new edm::ValueMap<float>());
0316 edm::ValueMap<float>::Filler fillerVtxMass(*vtxMassV);
0317 fillerVtxMass.insert(srcJet, vtxMass.begin(), vtxMass.end());
0318 fillerVtxMass.fill();
0319 iEvent.put(std::move(vtxMassV), "vtxMass");
0320
0321 std::unique_ptr<edm::ValueMap<float>> vtx3dLV(new edm::ValueMap<float>());
0322 edm::ValueMap<float>::Filler fillerVtx3dL(*vtx3dLV);
0323 fillerVtx3dL.insert(srcJet, vtx3dL.begin(), vtx3dL.end());
0324 fillerVtx3dL.fill();
0325 iEvent.put(std::move(vtx3dLV), "vtx3dL");
0326
0327 std::unique_ptr<edm::ValueMap<float>> vtx3deLV(new edm::ValueMap<float>());
0328 edm::ValueMap<float>::Filler fillerVtx3deL(*vtx3deLV);
0329 fillerVtx3deL.insert(srcJet, vtx3deL.begin(), vtx3deL.end());
0330 fillerVtx3deL.fill();
0331 iEvent.put(std::move(vtx3deLV), "vtx3deL");
0332
0333 std::unique_ptr<edm::ValueMap<int>> vtxNtrkV(new edm::ValueMap<int>());
0334 edm::ValueMap<int>::Filler fillerVtxNT(*vtxNtrkV);
0335 fillerVtxNT.insert(srcJet, vtxNtrk.begin(), vtxNtrk.end());
0336 fillerVtxNT.fill();
0337 iEvent.put(std::move(vtxNtrkV), "vtxNtrk");
0338
0339 std::unique_ptr<edm::ValueMap<float>> ptDV(new edm::ValueMap<float>());
0340 edm::ValueMap<float>::Filler fillerPtD(*ptDV);
0341 fillerPtD.insert(srcJet, ptD.begin(), ptD.end());
0342 fillerPtD.fill();
0343 iEvent.put(std::move(ptDV), "ptD");
0344
0345 std::unique_ptr<edm::ValueMap<float>> genptV(new edm::ValueMap<float>());
0346 edm::ValueMap<float>::Filler fillergenpt(*genptV);
0347 fillergenpt.insert(srcJet, genPtwNu.begin(), genPtwNu.end());
0348 fillergenpt.fill();
0349 iEvent.put(std::move(genptV), "genPtwNu");
0350 }
0351
0352 template <typename T>
0353 std::tuple<float, float, float> BJetEnergyRegressionVarProducer<T>::calculatePtRatioRel(edm::Ptr<reco::Candidate> lep,
0354 edm::Ptr<pat::Jet> jet) const {
0355 auto rawp4 = jet->correctedP4("Uncorrected");
0356 auto lepp4 = lep->p4();
0357
0358 if ((rawp4 - lepp4).R() < 1e-4)
0359 return std::tuple<float, float, float>(1.0, 0.0, 0.0);
0360
0361 auto jetp4 = (rawp4 - lepp4 * (1.0 / jet->jecFactor("L1FastJet"))) * (jet->pt() / rawp4.pt()) + lepp4;
0362 auto ptratio = lepp4.pt() / jetp4.pt();
0363 auto ptrel = lepp4.Vect().Cross((jetp4 - lepp4).Vect().Unit()).R();
0364 auto ptrelinv = (jetp4 - lepp4).Vect().Cross((lepp4).Vect().Unit()).R();
0365
0366 return std::tuple<float, float, float>(ptratio, ptrel, ptrelinv);
0367 }
0368
0369 template <typename T>
0370 std::tuple<float, float, float> BJetEnergyRegressionVarProducer<T>::calculatePtRatioRelSimple(
0371 edm::Ptr<reco::Candidate> lep, edm::Ptr<pat::Jet> jet) const {
0372 auto lepp4 = lep->p4();
0373 auto rawp4 = jet->correctedP4("Uncorrected");
0374
0375 if ((rawp4 - lepp4).R() < 1e-4)
0376 return std::tuple<float, float, float>(1.0, 0.0, 0.0);
0377
0378 auto jetp4 = jet->p4();
0379 auto ptratio = lepp4.pt() / jetp4.pt();
0380 auto ptrel = lepp4.Vect().Cross((jetp4).Vect().Unit()).R();
0381 auto ptrelinv = jetp4.Vect().Cross((lepp4).Vect().Unit()).R();
0382
0383 return std::tuple<float, float, float>(ptratio, ptrel, ptrelinv);
0384 }
0385
0386
0387 template <typename T>
0388 void BJetEnergyRegressionVarProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0389
0390
0391 edm::ParameterSetDescription desc;
0392 desc.add<edm::InputTag>("src")->setComment("jet input collection");
0393
0394
0395 desc.add<edm::InputTag>("pvsrc")->setComment("primary vertex input collection");
0396 desc.add<edm::InputTag>("svsrc")->setComment("secondary vertex input collection");
0397 desc.add<edm::InputTag>("gpsrc")->setComment("genparticles for nu recovery");
0398 std::string modname;
0399 if (typeid(T) == typeid(pat::Jet))
0400 modname += "Jet";
0401 modname += "RegressionVarProducer";
0402 descriptions.add(modname, desc);
0403 }
0404
0405 typedef BJetEnergyRegressionVarProducer<pat::Jet> JetRegressionVarProducer;
0406
0407
0408 DEFINE_FWK_MODULE(JetRegressionVarProducer);