NanoAOD/plugins/JetConstituentTableProducer.cc

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/stream/EDProducer.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/Utilities/interface/StreamID.h"
0009
0010 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0011 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0012 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0013 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0014 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
0015
0016 #include "DataFormats/PatCandidates/interface/Jet.h"
0017 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0018
0019 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0020
0021 #include "RecoBTag/FeatureTools/interface/TrackInfoBuilder.h"
0022 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0023
0024 #include "DataFormats/BTauReco/interface/TrackIPTagInfo.h"
0025 #include "DataFormats/BTauReco/interface/SecondaryVertexTagInfo.h"
0026 #include "RecoBTag/FeatureTools/interface/deep_helpers.h"
0027 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0028 using namespace btagbtvdeep;
0029
0030 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0031 #include "DataFormats/NanoAOD/interface/FlatTable.h"
0032
0033 template <typename T>
0034 class JetConstituentTableProducer : public edm::stream::EDProducer<> {
0035 public:
0036   explicit JetConstituentTableProducer(const edm::ParameterSet &);
0037   ~JetConstituentTableProducer() override;
0038
0039   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0040
0041 private:
0042   void produce(edm::Event &, const edm::EventSetup &) override;
0043
0044   typedef reco::VertexCollection VertexCollection;
0045   //=====
0046   typedef reco::VertexCompositePtrCandidateCollection SVCollection;
0047
0048   //const std::string name_;
0049   const std::string name_;
0050   const std::string nameSV_;
0051   const std::string idx_name_;
0052   const std::string idx_nameSV_;
0053   const bool readBtag_;
0054   const double jet_radius_;
0055
0056   edm::EDGetTokenT<edm::View<T>> jet_token_;
0057   edm::EDGetTokenT<VertexCollection> vtx_token_;
0058   edm::EDGetTokenT<reco::CandidateView> cand_token_;
0059   edm::EDGetTokenT<SVCollection> sv_token_;
0060
0061   edm::Handle<VertexCollection> vtxs_;
0062   edm::Handle<reco::CandidateView> cands_;
0063   edm::Handle<SVCollection> svs_;
0064   edm::ESHandle<TransientTrackBuilder> track_builder_;
0065   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> track_builder_token_;
0066
0067   const reco::Vertex *pv_ = nullptr;
0068 };
0069
0070 //
0071 // constructors and destructor
0072 //
0073 template <typename T>
0074 JetConstituentTableProducer<T>::JetConstituentTableProducer(const edm::ParameterSet &iConfig)
0075     : name_(iConfig.getParameter<std::string>("name")),
0076       nameSV_(iConfig.getParameter<std::string>("nameSV")),
0077       idx_name_(iConfig.getParameter<std::string>("idx_name")),
0078       idx_nameSV_(iConfig.getParameter<std::string>("idx_nameSV")),
0079       readBtag_(iConfig.getParameter<bool>("readBtag")),
0080       jet_radius_(iConfig.getParameter<double>("jet_radius")),
0081       jet_token_(consumes<edm::View<T>>(iConfig.getParameter<edm::InputTag>("jets"))),
0082       vtx_token_(consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0083       cand_token_(consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("candidates"))),
0084       sv_token_(consumes<SVCollection>(iConfig.getParameter<edm::InputTag>("secondary_vertices"))),
0085       track_builder_token_(
0086           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))) {
0087   produces<nanoaod::FlatTable>(name_);
0088   produces<nanoaod::FlatTable>(nameSV_);
0089   produces<std::vector<reco::CandidatePtr>>();
0090 }
0091
0092 template <typename T>
0093 JetConstituentTableProducer<T>::~JetConstituentTableProducer() {}
0094
0095 template <typename T>
0096 void JetConstituentTableProducer<T>::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0097   // elements in all these collections must have the same order!
0098   auto outCands = std::make_unique<std::vector<reco::CandidatePtr>>();
0099   auto outSVs = std::make_unique<std::vector<const reco::VertexCompositePtrCandidate *>>();
0100   std::vector<int> jetIdx_pf, jetIdx_sv, pfcandIdx, svIdx;
0101   // PF Cands
0102   std::vector<float> btagEtaRel, btagPtRatio, btagPParRatio, btagSip3dVal, btagSip3dSig, btagJetDistVal,
0103       btagDecayLenVal, cand_pt, cand_dzFromPV, cand_dxyFromPV, cand_dzErrFromPV, cand_dxyErrFromPV;
0104   // Secondary vertices
0105   std::vector<float> sv_mass, sv_pt, sv_ntracks, sv_chi2, sv_normchi2, sv_dxy, sv_dxysig, sv_d3d, sv_d3dsig,
0106       sv_costhetasvpv;
0107   std::vector<float> sv_ptrel, sv_phirel, sv_deltaR, sv_enratio;
0108
0109   auto jets = iEvent.getHandle(jet_token_);
0110   iEvent.getByToken(vtx_token_, vtxs_);
0111   iEvent.getByToken(cand_token_, cands_);
0112   iEvent.getByToken(sv_token_, svs_);
0113
0114   if (readBtag_) {
0115     track_builder_ = iSetup.getHandle(track_builder_token_);
0116   }
0117
0118   for (unsigned i_jet = 0; i_jet < jets->size(); ++i_jet) {
0119     const auto &jet = jets->at(i_jet);
0120     math::XYZVector jet_dir = jet.momentum().Unit();
0121     GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz());
0122     VertexDistance3D vdist;
0123
0124     pv_ = &vtxs_->at(0);
0125
0126     //////////////////////
0127     // Secondary Vertices
0128     std::vector<const reco::VertexCompositePtrCandidate *> jetSVs;
0129     std::vector<const reco::VertexCompositePtrCandidate *> allSVs;
0130     for (const auto &sv : *svs_) {
0131       // Factor in cuts in NanoAOD for indexing
0132       Measurement1D dl = vdist.distance(
0133           vtxs_->front(), VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0134       if (dl.significance() > 3) {
0135         allSVs.push_back(&sv);
0136       }
0137       if (reco::deltaR2(sv, jet) < jet_radius_ * jet_radius_) {
0138         jetSVs.push_back(&sv);
0139       }
0140     }
0141     // sort by dxy significance
0142     std::sort(jetSVs.begin(),
0143               jetSVs.end(),
0144               [&](const reco::VertexCompositePtrCandidate *sva, const reco::VertexCompositePtrCandidate *svb) {
0145                 return sv_vertex_comparator(*sva, *svb, *pv_);
0146               });
0147
0148     for (const auto &sv : jetSVs) {
0149       // auto svPtrs = svs_->ptrs();
0150       auto svInNewList = std::find(allSVs.begin(), allSVs.end(), sv);
0151       if (svInNewList == allSVs.end()) {
0152         // continue;
0153         svIdx.push_back(-1);
0154       } else {
0155         svIdx.push_back(svInNewList - allSVs.begin());
0156       }
0157       outSVs->push_back(sv);
0158       jetIdx_sv.push_back(i_jet);
0159       if (readBtag_ && !vtxs_->empty()) {
0160         // Jet independent
0161         sv_mass.push_back(sv->mass());
0162         sv_pt.push_back(sv->pt());
0163
0164         sv_ntracks.push_back(sv->numberOfDaughters());
0165         sv_chi2.push_back(sv->vertexChi2());
0166         sv_normchi2.push_back(catch_infs_and_bound(sv->vertexChi2() / sv->vertexNdof(), 1000, -1000, 1000));
0167         const auto &dxy_meas = vertexDxy(*sv, *pv_);
0168         sv_dxy.push_back(dxy_meas.value());
0169         sv_dxysig.push_back(catch_infs_and_bound(dxy_meas.value() / dxy_meas.error(), 0, -1, 800));
0170         const auto &d3d_meas = vertexD3d(*sv, *pv_);
0171         sv_d3d.push_back(d3d_meas.value());
0172         sv_d3dsig.push_back(catch_infs_and_bound(d3d_meas.value() / d3d_meas.error(), 0, -1, 800));
0173         sv_costhetasvpv.push_back(vertexDdotP(*sv, *pv_));
0174         // Jet related
0175         sv_ptrel.push_back(sv->pt() / jet.pt());
0176         sv_phirel.push_back(reco::deltaPhi(*sv, jet));
0177         sv_deltaR.push_back(catch_infs_and_bound(std::fabs(reco::deltaR(*sv, jet_dir)) - 0.5, 0, -2, 0));
0178         sv_enratio.push_back(sv->energy() / jet.energy());
0179       }
0180     }
0181
0182     // PF Cands
0183     std::vector<reco::CandidatePtr> const &daughters = jet.daughterPtrVector();
0184
0185     for (const auto &cand : daughters) {
0186       auto candPtrs = cands_->ptrs();
0187       auto candInNewList = std::find(candPtrs.begin(), candPtrs.end(), cand);
0188       if (candInNewList == candPtrs.end()) {
0189         //std::cout << "Cannot find candidate : " << cand.id() << ", " << cand.key() << ", pt = " << cand->pt() << std::endl;
0190         continue;
0191       }
0192       outCands->push_back(cand);
0193       jetIdx_pf.push_back(i_jet);
0194       pfcandIdx.push_back(candInNewList - candPtrs.begin());
0195       cand_pt.push_back(cand->pt());
0196       auto const *packedCand = dynamic_cast<pat::PackedCandidate const *>(cand.get());
0197       if (packedCand && packedCand->hasTrackDetails()) {
0198         const reco::Track *track_ptr = &(packedCand->pseudoTrack());
0199         cand_dzFromPV.push_back(track_ptr->dz(pv_->position()));
0200         cand_dxyFromPV.push_back(track_ptr->dxy(pv_->position()));
0201         cand_dzErrFromPV.push_back(std::hypot(track_ptr->dzError(), pv_->zError()));
0202         cand_dxyErrFromPV.push_back(track_ptr->dxyError(pv_->position(), pv_->covariance()));
0203       } else {
0204         cand_dzFromPV.push_back(-1);
0205         cand_dxyFromPV.push_back(-1);
0206         cand_dzErrFromPV.push_back(-1);
0207         cand_dxyErrFromPV.push_back(-1);
0208       }
0209
0210       if (readBtag_ && !vtxs_->empty()) {
0211         if (cand.isNull())
0212           continue;
0213         auto const *packedCand = dynamic_cast<pat::PackedCandidate const *>(cand.get());
0214         if (packedCand == nullptr)
0215           continue;
0216         if (packedCand && packedCand->hasTrackDetails()) {
0217           btagbtvdeep::TrackInfoBuilder trkinfo(track_builder_);
0218           trkinfo.buildTrackInfo(&(*packedCand), jet_dir, jet_ref_track_dir, vtxs_->at(0));
0219           btagEtaRel.push_back(trkinfo.getTrackEtaRel());
0220           btagPtRatio.push_back(trkinfo.getTrackPtRatio());
0221           btagPParRatio.push_back(trkinfo.getTrackPParRatio());
0222           btagSip3dVal.push_back(trkinfo.getTrackSip3dVal());
0223           btagSip3dSig.push_back(trkinfo.getTrackSip3dSig());
0224           btagJetDistVal.push_back(trkinfo.getTrackJetDistVal());
0225           // decay length
0226           const reco::Track *track_ptr = packedCand->bestTrack();
0227           reco::TransientTrack transient_track = track_builder_->build(track_ptr);
0228           double decayLength = -1;
0229           TrajectoryStateOnSurface closest = IPTools::closestApproachToJet(
0230               transient_track.impactPointState(), *pv_, jet_ref_track_dir, transient_track.field());
0231           if (closest.isValid())
0232             decayLength = (closest.globalPosition() - RecoVertex::convertPos(pv_->position())).mag();
0233           else
0234             decayLength = -1;
0235           btagDecayLenVal.push_back(decayLength);
0236         } else {
0237           btagEtaRel.push_back(0);
0238           btagPtRatio.push_back(0);
0239           btagPParRatio.push_back(0);
0240           btagSip3dVal.push_back(0);
0241           btagSip3dSig.push_back(0);
0242           btagJetDistVal.push_back(0);
0243           btagDecayLenVal.push_back(0);
0244         }
0245       }
0246     }  // end jet loop
0247   }
0248
0249   auto candTable = std::make_unique<nanoaod::FlatTable>(outCands->size(), name_, false);
0250   // We fill from here only stuff that cannot be created with the SimpleFlatTableProducer
0251   candTable->addColumn<int>(idx_name_, pfcandIdx, "Index in the candidate list");
0252   candTable->addColumn<int>("jetIdx", jetIdx_pf, "Index of the parent jet");
0253   if (readBtag_) {
0254     candTable->addColumn<float>("pt", cand_pt, "pt", 10);  // to check matchind down the line
0255     candTable->addColumn<float>("dzFromPV", cand_dzFromPV, "dzFromPV", 10);
0256     candTable->addColumn<float>("dxyFromPV", cand_dxyFromPV, "dxyFromPV", 10);
0257     candTable->addColumn<float>("dzErrFromPV", cand_dzErrFromPV, "dzErrFromPV", 10);
0258     candTable->addColumn<float>("dxyErrFromPV", cand_dxyErrFromPV, "dxyErrFromPV", 10);
0259     candTable->addColumn<float>("btagEtaRel", btagEtaRel, "btagEtaRel", 10);
0260     candTable->addColumn<float>("btagPtRatio", btagPtRatio, "btagPtRatio", 10);
0261     candTable->addColumn<float>("btagPParRatio", btagPParRatio, "btagPParRatio", 10);
0262     candTable->addColumn<float>("btagSip3dVal", btagSip3dVal, "btagSip3dVal", 10);
0263     candTable->addColumn<float>("btagSip3dSig", btagSip3dSig, "btagSip3dSig", 10);
0264     candTable->addColumn<float>("btagJetDistVal", btagJetDistVal, "btagJetDistVal", 10);
0265     candTable->addColumn<float>("btagDecayLenVal", btagDecayLenVal, "btagDecayLenVal", 10);
0266   }
0267   iEvent.put(std::move(candTable), name_);
0268
0269   // SV table
0270   auto svTable = std::make_unique<nanoaod::FlatTable>(outSVs->size(), nameSV_, false);
0271   // We fill from here only stuff that cannot be created with the SimpleFlatTnameableProducer
0272   svTable->addColumn<int>("jetIdx", jetIdx_sv, "Index of the parent jet");
0273   svTable->addColumn<int>(idx_nameSV_, svIdx, "Index in the SV list");
0274   if (readBtag_) {
0275     svTable->addColumn<float>("mass", sv_mass, "SV mass", 10);
0276     svTable->addColumn<float>("pt", sv_pt, "SV pt", 10);
0277     svTable->addColumn<float>("ntracks", sv_ntracks, "Number of tracks associated to SV", 10);
0278     svTable->addColumn<float>("chi2", sv_chi2, "chi2", 10);
0279     svTable->addColumn<float>("normchi2", sv_normchi2, "chi2/ndof", 10);
0280     svTable->addColumn<float>("dxy", sv_dxy, "", 10);
0281     svTable->addColumn<float>("dxysig", sv_dxysig, "", 10);
0282     svTable->addColumn<float>("d3d", sv_d3d, "", 10);
0283     svTable->addColumn<float>("d3dsig", sv_d3dsig, "", 10);
0284     svTable->addColumn<float>("costhetasvpv", sv_costhetasvpv, "", 10);
0285     // Jet related
0286     svTable->addColumn<float>("phirel", sv_phirel, "DeltaPhi(sv, jet)", 10);
0287     svTable->addColumn<float>("ptrel", sv_ptrel, "pT relative to parent jet", 10);
0288     svTable->addColumn<float>("deltaR", sv_deltaR, "dR from parent jet", 10);
0289     svTable->addColumn<float>("enration", sv_enratio, "energy relative to parent jet", 10);
0290   }
0291   iEvent.put(std::move(svTable), nameSV_);
0292
0293   iEvent.put(std::move(outCands));
0294 }
0295
0296 template <typename T>
0297 void JetConstituentTableProducer<T>::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0298   edm::ParameterSetDescription desc;
0299   desc.add<std::string>("name", "JetPFCands");
0300   desc.add<std::string>("nameSV", "JetSV");
0301   desc.add<std::string>("idx_name", "candIdx");
0302   desc.add<std::string>("idx_nameSV", "svIdx");
0303   desc.add<double>("jet_radius", true);
0304   desc.add<bool>("readBtag", true);
0305   desc.add<edm::InputTag>("jets", edm::InputTag("slimmedJetsAK8"));
0306   desc.add<edm::InputTag>("vertices", edm::InputTag("offlineSlimmedPrimaryVertices"));
0307   desc.add<edm::InputTag>("candidates", edm::InputTag("packedPFCandidates"));
0308   desc.add<edm::InputTag>("secondary_vertices", edm::InputTag("slimmedSecondaryVertices"));
0309   descriptions.addWithDefaultLabel(desc);
0310 }
0311
0312 typedef JetConstituentTableProducer<pat::Jet> PatJetConstituentTableProducer;
0313 typedef JetConstituentTableProducer<reco::GenJet> GenJetConstituentTableProducer;
0314
0315 DEFINE_FWK_MODULE(PatJetConstituentTableProducer);
0316 DEFINE_FWK_MODULE(GenJetConstituentTableProducer);