Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-26 01:51:27 

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 namespace {
0034 
0035   std::string toupper(std::string s) {
0036     std::transform(s.begin(), s.end(), s.begin(), ::toupper);
0037     return s;
0038   }
0039 
0040 }  // namespace
0041 
0042 using std::toupper;
0043 
0044 template <typename T>
0045 class JetConstituentTableProducer : public edm::stream::EDProducer<> {
0046 public:
0047   explicit JetConstituentTableProducer(const edm::ParameterSet &);
0048   ~JetConstituentTableProducer() override;
0049 
0050   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0051 
0052 private:
0053   void produce(edm::Event &, const edm::EventSetup &) override;
0054 
0055   typedef reco::VertexCollection VertexCollection;
0056   //=====
0057   typedef reco::VertexCompositePtrCandidateCollection SVCollection;
0058 
0059   const std::string name_;
0060   const std::string nameSV_;
0061   const std::string idx_name_;
0062   const std::string idx_nameSV_;
0063   const bool readBtag_;
0064   const double jet_radius_;
0065 
0066   edm::EDGetTokenT<edm::View<T>> jet_token_;
0067   edm::EDGetTokenT<VertexCollection> vtx_token_;
0068   edm::EDGetTokenT<reco::CandidateView> cand_token_;
0069   edm::EDGetTokenT<SVCollection> sv_token_;
0070 
0071   edm::Handle<VertexCollection> vtxs_;
0072   edm::Handle<reco::CandidateView> cands_;
0073   edm::Handle<SVCollection> svs_;
0074   edm::ESHandle<TransientTrackBuilder> track_builder_;
0075   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> track_builder_token_;
0076 
0077   const std::string sv_sort_;
0078   const std::string pf_sort_;
0079 
0080   const reco::Vertex *pv_ = nullptr;
0081 };
0082 
0083 //
0084 // constructors and destructor
0085 //
0086 template <typename T>
0087 JetConstituentTableProducer<T>::JetConstituentTableProducer(const edm::ParameterSet &iConfig)
0088     : name_(iConfig.getParameter<std::string>("name")),
0089       nameSV_(iConfig.getParameter<std::string>("nameSV")),
0090       idx_name_(iConfig.getParameter<std::string>("idx_name")),
0091       idx_nameSV_(iConfig.getParameter<std::string>("idx_nameSV")),
0092       readBtag_(iConfig.getParameter<bool>("readBtag")),
0093       jet_radius_(iConfig.getParameter<double>("jet_radius")),
0094       jet_token_(consumes<edm::View<T>>(iConfig.getParameter<edm::InputTag>("jets"))),
0095       vtx_token_(consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0096       cand_token_(consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("candidates"))),
0097       sv_token_(consumes<SVCollection>(iConfig.getParameter<edm::InputTag>("secondary_vertices"))),
0098       track_builder_token_(
0099           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
0100       sv_sort_(iConfig.getUntrackedParameter<std::string>("sv_sort")),
0101       pf_sort_(iConfig.getUntrackedParameter<std::string>("pf_sort")) {
0102   produces<nanoaod::FlatTable>(name_);
0103   produces<nanoaod::FlatTable>(nameSV_);
0104   produces<std::vector<reco::CandidatePtr>>();
0105   std::clog << "sv_sort: " << sv_sort_ << std::endl;
0106   std::clog << "pf_sort: " << pf_sort_ << std::endl;
0107 }
0108 
0109 template <typename T>
0110 JetConstituentTableProducer<T>::~JetConstituentTableProducer() {}
0111 
0112 template <typename T>
0113 void JetConstituentTableProducer<T>::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0114   // elements in all these collections must have the same order!
0115   auto outCands = std::make_unique<std::vector<reco::CandidatePtr>>();
0116   auto outSVs = std::make_unique<std::vector<const reco::VertexCompositePtrCandidate *>>();
0117   std::vector<int> jetIdx_pf, jetIdx_sv, pfcandIdx, svIdx;
0118   // PF Cands
0119   std::vector<float> btagEtaRel, btagPtRatio, btagPParRatio, btagSip3dVal, btagSip3dSig, btagJetDistVal,
0120       btagDecayLenVal, cand_pt, cand_dzFromPV, cand_dxyFromPV, cand_dzErrFromPV, cand_dxyErrFromPV;
0121   std::vector<int> cand_jetSVIdx;
0122   // Secondary vertices
0123   std::vector<float> sv_mass, sv_pt, sv_ntracks, sv_chi2, sv_normchi2, sv_dxy, sv_dxysig, sv_d3d, sv_d3dsig,
0124       sv_costhetasvpv;
0125   std::vector<float> sv_ptrel, sv_phirel, sv_deltaR, sv_enratio;
0126 
0127   auto jets = iEvent.getHandle(jet_token_);
0128   iEvent.getByToken(vtx_token_, vtxs_);
0129   iEvent.getByToken(cand_token_, cands_);
0130   iEvent.getByToken(sv_token_, svs_);
0131 
0132   if (readBtag_) {
0133     track_builder_ = iSetup.getHandle(track_builder_token_);
0134   }
0135 
0136   for (unsigned i_jet = 0; i_jet < jets->size(); ++i_jet) {
0137     const auto &jet = jets->at(i_jet);
0138     math::XYZVector jet_dir = jet.momentum().Unit();
0139     GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz());
0140     VertexDistance3D vdist;
0141 
0142     pv_ = &vtxs_->at(0);
0143 
0144     //////////////////////
0145     // Secondary Vertices
0146     std::vector<const reco::VertexCompositePtrCandidate *> jetSVs;
0147     std::vector<const reco::VertexCompositePtrCandidate *> allSVs;
0148     for (const auto &sv : *svs_) {
0149       // Factor in cuts in NanoAOD for indexing
0150       // TODO: seems fragile, better use NanoAOD vertexTable instead of slimmedSecondaryVertices as input?
0151       Measurement1D dl = vdist.distance(
0152           vtxs_->front(), VertexState(RecoVertex::convertPos(sv.position()), RecoVertex::convertError(sv.error())));
0153       if (dl.significance() > 3) {
0154         allSVs.push_back(&sv);
0155       }
0156       if (reco::deltaR2(sv, jet) < jet_radius_ * jet_radius_) {
0157         jetSVs.push_back(&sv);
0158       }
0159     }
0160     if (toupper(sv_sort_) == "IP") {
0161       // sort by dxy significance
0162       std::sort(jetSVs.begin(),
0163                 jetSVs.end(),
0164                 [this](const reco::VertexCompositePtrCandidate *sva, const reco::VertexCompositePtrCandidate *svb) {
0165                   return sv_vertex_comparator(*sva, *svb, *pv_);
0166                 });
0167     } else if (toupper(sv_sort_) == "PT") {
0168       std::sort(jetSVs.begin(),
0169                 jetSVs.end(),
0170                 [](const reco::VertexCompositePtrCandidate *sva, const reco::VertexCompositePtrCandidate *svb) {
0171                   return sva->pt() > svb->pt();
0172                 });
0173     } else if (!sv_sort_.empty()) {
0174       throw cms::Exception("Configuration")
0175           << "Unknown sorting option for secondary vertices: " << sv_sort_ << std::endl;
0176     }
0177 
0178     for (const auto &sv : jetSVs) {
0179       // auto svPtrs = svs_->ptrs();
0180       auto svInNewList = std::find(allSVs.begin(), allSVs.end(), sv);
0181       if (svInNewList == allSVs.end()) {
0182         // continue;
0183         svIdx.push_back(-1);
0184       } else {
0185         svIdx.push_back(svInNewList - allSVs.begin());
0186       }
0187       outSVs->push_back(sv);
0188       jetIdx_sv.push_back(i_jet);
0189       if (readBtag_ && !vtxs_->empty()) {
0190         // Jet independent
0191         sv_mass.push_back(sv->mass());
0192         sv_pt.push_back(sv->pt());
0193 
0194         sv_ntracks.push_back(sv->numberOfDaughters());
0195         sv_chi2.push_back(sv->vertexChi2());
0196         sv_normchi2.push_back(catch_infs_and_bound(sv->vertexChi2() / sv->vertexNdof(), 1000, -1000, 1000));
0197         const auto &dxy_meas = vertexDxy(*sv, *pv_);
0198         sv_dxy.push_back(dxy_meas.value());
0199         sv_dxysig.push_back(catch_infs_and_bound(dxy_meas.value() / dxy_meas.error(), 0, -1, 800));
0200         const auto &d3d_meas = vertexD3d(*sv, *pv_);
0201         sv_d3d.push_back(d3d_meas.value());
0202         sv_d3dsig.push_back(catch_infs_and_bound(d3d_meas.value() / d3d_meas.error(), 0, -1, 800));
0203         sv_costhetasvpv.push_back(vertexDdotP(*sv, *pv_));
0204         // Jet related
0205         sv_ptrel.push_back(sv->pt() / jet.pt());
0206         sv_phirel.push_back(reco::deltaPhi(*sv, jet));
0207         sv_deltaR.push_back(catch_infs_and_bound(std::fabs(reco::deltaR(*sv, jet_dir)) - 0.5, 0, -2, 0));
0208         sv_enratio.push_back(sv->energy() / jet.energy());
0209       }
0210     }
0211 
0212     // PF Cands
0213     std::vector<reco::CandidatePtr> const &daughters = jet.daughterPtrVector();
0214     std::vector<size_t> dauidx;
0215     dauidx.reserve(daughters.size());
0216     for (size_t i = 0; i < daughters.size(); ++i)
0217       dauidx.push_back(i);
0218     if (toupper(pf_sort_) == "PT") {
0219       std::sort(dauidx.begin(), dauidx.end(), [&daughters](size_t i, size_t j) {
0220         return daughters[i]->pt() > daughters[j]->pt();
0221       });
0222     } else if (!pf_sort_.empty()) {
0223       throw cms::Exception("Configuration")
0224           << "Unknown sorting option for particle flow candidates: " << pf_sort_ << std::endl;
0225     }
0226 
0227     for (size_t di : dauidx) {
0228       const auto &cand = daughters[di];
0229       auto candPtrs = cands_->ptrs();
0230       auto candInNewList = std::find(candPtrs.begin(), candPtrs.end(), cand);
0231       if (candInNewList == candPtrs.end()) {
0232         //std::cout << "Cannot find candidate : " << cand.id() << ", " << cand.key() << ", pt = " << cand->pt() << std::endl;
0233         continue;
0234       }
0235       outCands->push_back(cand);
0236       jetIdx_pf.push_back(i_jet);
0237       pfcandIdx.push_back(candInNewList - candPtrs.begin());
0238       cand_pt.push_back(cand->pt());
0239       auto const *packedCand = dynamic_cast<pat::PackedCandidate const *>(cand.get());
0240       if (packedCand && packedCand->hasTrackDetails()) {
0241         const reco::Track *track_ptr = &(packedCand->pseudoTrack());
0242         cand_dzFromPV.push_back(track_ptr->dz(pv_->position()));
0243         cand_dxyFromPV.push_back(track_ptr->dxy(pv_->position()));
0244         cand_dzErrFromPV.push_back(std::hypot(track_ptr->dzError(), pv_->zError()));
0245         cand_dxyErrFromPV.push_back(track_ptr->dxyError(pv_->position(), pv_->covariance()));
0246       } else {
0247         cand_dzFromPV.push_back(-1);
0248         cand_dxyFromPV.push_back(-1);
0249         cand_dzErrFromPV.push_back(-1);
0250         cand_dxyErrFromPV.push_back(-1);
0251       }
0252 
0253       if (readBtag_ && !vtxs_->empty()) {
0254         if (cand.isNull())
0255           continue;
0256         auto const *packedCand = dynamic_cast<pat::PackedCandidate const *>(cand.get());
0257         if (packedCand == nullptr)
0258           continue;
0259         if (packedCand && packedCand->hasTrackDetails()) {
0260           btagbtvdeep::TrackInfoBuilder trkinfo(track_builder_);
0261           trkinfo.buildTrackInfo(&(*packedCand), jet_dir, jet_ref_track_dir, vtxs_->at(0));
0262           btagEtaRel.push_back(trkinfo.getTrackEtaRel());
0263           btagPtRatio.push_back(trkinfo.getTrackPtRatio());
0264           btagPParRatio.push_back(trkinfo.getTrackPParRatio());
0265           btagSip3dVal.push_back(trkinfo.getTrackSip3dVal());
0266           btagSip3dSig.push_back(trkinfo.getTrackSip3dSig());
0267           btagJetDistVal.push_back(trkinfo.getTrackJetDistVal());
0268           // decay length
0269           const reco::Track *track_ptr = packedCand->bestTrack();
0270           reco::TransientTrack transient_track = track_builder_->build(track_ptr);
0271           double decayLength = -1;
0272           TrajectoryStateOnSurface closest = IPTools::closestApproachToJet(
0273               transient_track.impactPointState(), *pv_, jet_ref_track_dir, transient_track.field());
0274           if (closest.isValid())
0275             decayLength = (closest.globalPosition() - RecoVertex::convertPos(pv_->position())).mag();
0276           else
0277             decayLength = -1;
0278           btagDecayLenVal.push_back(decayLength);
0279           // Associate PF candidates to secondary vertices (SVs) by matching their tracks
0280           int jsvMatchIndex = -1;
0281           int jsvIndex = 0;
0282           for (const auto &sv : *outSVs) {
0283             for (const auto &track : sv->daughterPtrVector()) {
0284               double eps = 1e-3;
0285               double dR = deltaR(track->eta(), track->phi(), cand->eta(), cand->phi());
0286               if (dR < eps && abs(track->pt() - cand->pt()) < eps) {
0287                 jsvMatchIndex = jsvIndex;
0288                 break;
0289               }
0290             }
0291             if (jsvMatchIndex >= 0)
0292               break;
0293             jsvIndex++;
0294           }
0295           cand_jetSVIdx.push_back(jsvMatchIndex);
0296         } else {
0297           btagEtaRel.push_back(0);
0298           btagPtRatio.push_back(0);
0299           btagPParRatio.push_back(0);
0300           btagSip3dVal.push_back(0);
0301           btagSip3dSig.push_back(0);
0302           btagJetDistVal.push_back(0);
0303           btagDecayLenVal.push_back(0);
0304           cand_jetSVIdx.push_back(-1);
0305         }
0306       }
0307     }  // end jet loop
0308   }
0309 
0310   auto candTable = std::make_unique<nanoaod::FlatTable>(outCands->size(), name_, false);
0311   // We fill from here only stuff that cannot be created with the SimpleFlatTableProducer
0312   candTable->addColumn<int>(idx_name_, pfcandIdx, "Index in the candidate list");
0313   candTable->addColumn<int>("jetIdx", jetIdx_pf, "Index of the parent jet");
0314   candTable->addColumn<float>("pt", cand_pt, "pt", 10);  // to check matching down the line
0315   if (readBtag_) {
0316     candTable->addColumn<float>("dzFromPV", cand_dzFromPV, "dzFromPV", 10);
0317     candTable->addColumn<float>("dxyFromPV", cand_dxyFromPV, "dxyFromPV", 10);
0318     candTable->addColumn<float>("dzErrFromPV", cand_dzErrFromPV, "dzErrFromPV", 10);
0319     candTable->addColumn<float>("dxyErrFromPV", cand_dxyErrFromPV, "dxyErrFromPV", 10);
0320     candTable->addColumn<float>("btagEtaRel", btagEtaRel, "btagEtaRel", 10);
0321     candTable->addColumn<float>("btagPtRatio", btagPtRatio, "btagPtRatio", 10);
0322     candTable->addColumn<float>("btagPParRatio", btagPParRatio, "btagPParRatio", 10);
0323     candTable->addColumn<float>("btagSip3dVal", btagSip3dVal, "btagSip3dVal", 10);
0324     candTable->addColumn<float>("btagSip3dSig", btagSip3dSig, "btagSip3dSig", 10);
0325     candTable->addColumn<float>("btagJetDistVal", btagJetDistVal, "btagJetDistVal", 10);
0326     candTable->addColumn<float>("btagDecayLenVal", btagDecayLenVal, "btagDecayLenVal", 10);
0327     candTable->addColumn<int>("jetSVIdx", cand_jetSVIdx, "Index of the parent in the " + nameSV_ + " list");
0328   }
0329   iEvent.put(std::move(candTable), name_);
0330 
0331   // SV table
0332   auto svTable = std::make_unique<nanoaod::FlatTable>(outSVs->size(), nameSV_, false);
0333   // We fill from here only stuff that cannot be created with the SimpleFlatTnameableProducer
0334   svTable->addColumn<int>("jetIdx", jetIdx_sv, "Index of the parent jet");
0335   svTable->addColumn<int>(idx_nameSV_, svIdx, "Index in the SV list");
0336   if (readBtag_) {
0337     svTable->addColumn<float>("mass", sv_mass, "SV mass", 10);
0338     svTable->addColumn<float>("pt", sv_pt, "SV pt", 10);
0339     svTable->addColumn<float>("ntracks", sv_ntracks, "Number of tracks associated to SV", 10);
0340     svTable->addColumn<float>("chi2", sv_chi2, "chi2", 10);
0341     svTable->addColumn<float>("normchi2", sv_normchi2, "chi2/ndof", 10);
0342     svTable->addColumn<float>("dxy", sv_dxy, "", 10);
0343     svTable->addColumn<float>("dxysig", sv_dxysig, "", 10);
0344     svTable->addColumn<float>("d3d", sv_d3d, "", 10);
0345     svTable->addColumn<float>("d3dsig", sv_d3dsig, "", 10);
0346     svTable->addColumn<float>("costhetasvpv", sv_costhetasvpv, "", 10);
0347     // Jet related
0348     svTable->addColumn<float>("phirel", sv_phirel, "DeltaPhi(sv, jet)", 10);
0349     svTable->addColumn<float>("ptrel", sv_ptrel, "pT relative to parent jet", 10);
0350     svTable->addColumn<float>("deltaR", sv_deltaR, "dR from parent jet - 0.5", 10);
0351     svTable->addColumn<float>("enration", sv_enratio, "energy relative to parent jet", 10);
0352   }
0353   iEvent.put(std::move(svTable), nameSV_);
0354 
0355   iEvent.put(std::move(outCands));
0356 }
0357 
0358 template <typename T>
0359 void JetConstituentTableProducer<T>::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0360   edm::ParameterSetDescription desc;
0361   desc.add<std::string>("name", "JetPFCands");
0362   desc.add<std::string>("nameSV", "JetSV");
0363   desc.add<std::string>("idx_name", "candIdx");
0364   desc.add<std::string>("idx_nameSV", "svIdx");
0365   desc.add<double>("jet_radius", true);
0366   desc.add<bool>("readBtag", true);
0367   desc.add<edm::InputTag>("jets", edm::InputTag("slimmedJetsAK8"));
0368   desc.add<edm::InputTag>("vertices", edm::InputTag("offlineSlimmedPrimaryVertices"));
0369   desc.add<edm::InputTag>("candidates", edm::InputTag("packedPFCandidates"));
0370   desc.add<edm::InputTag>("secondary_vertices", edm::InputTag("slimmedSecondaryVertices"));
0371   desc.addUntracked<std::string>("sv_sort", "IP");
0372   desc.addUntracked<std::string>("pf_sort", "");
0373   descriptions.addWithDefaultLabel(desc);
0374 }
0375 
0376 typedef JetConstituentTableProducer<pat::Jet> PatJetConstituentTableProducer;
0377 typedef JetConstituentTableProducer<reco::GenJet> GenJetConstituentTableProducer;
0378 
0379 DEFINE_FWK_MODULE(PatJetConstituentTableProducer);
0380 DEFINE_FWK_MODULE(GenJetConstituentTableProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!