Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-14 23:39:46

0001 #include <TVector3.h>
0002 
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005 
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/Utilities/interface/StreamID.h"
0011 #include "FWCore/Utilities/interface/ESGetToken.h"
0012 
0013 #include "DataFormats/Candidate/interface/Candidate.h"
0014 #include "DataFormats/PatCandidates/interface/Jet.h"
0015 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0016 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0017 
0018 #include "RecoBTag/FeatureTools/interface/TrackInfoBuilder.h"
0019 #include "RecoBTag/FeatureTools/interface/deep_helpers.h"
0020 #include "RecoBTag/FeatureTools/interface/sorting_modules.h"
0021 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0022 
0023 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0024 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0025 
0026 #include "DataFormats/BTauReco/interface/DeepBoostedJetTagInfo.h"
0027 
0028 using namespace btagbtvdeep;
0029 
0030 class DeepBoostedJetTagInfoProducer : public edm::stream::EDProducer<> {
0031 public:
0032   explicit DeepBoostedJetTagInfoProducer(const edm::ParameterSet &);
0033   ~DeepBoostedJetTagInfoProducer() override;
0034 
0035   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0036 
0037 private:
0038   typedef std::vector<reco::DeepBoostedJetTagInfo> DeepBoostedJetTagInfoCollection;
0039   typedef reco::VertexCompositePtrCandidateCollection SVCollection;
0040   typedef reco::VertexCollection VertexCollection;
0041   typedef edm::View<reco::Candidate> CandidateView;
0042 
0043   void beginStream(edm::StreamID) override {}
0044   void produce(edm::Event &, const edm::EventSetup &) override;
0045   void endStream() override {}
0046 
0047   void fillParticleFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet);
0048   void fillSVFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet);
0049   void fillParticleFeaturesHLT(DeepBoostedJetFeatures &fts, const reco::Jet &jet, const reco::VertexRefProd &PVRefProd);
0050   void fillSVFeaturesHLT(DeepBoostedJetFeatures &fts, const reco::Jet &jet);
0051 
0052   float puppiWgt(const reco::CandidatePtr &cand);
0053   bool useTrackProperties(const reco::PFCandidate *reco_cand);
0054 
0055   const double jet_radius_;
0056   const double min_jet_pt_;
0057   const double max_jet_eta_;
0058   const double min_pt_for_track_properties_;
0059   const double min_pt_for_pfcandidates_;
0060   const bool use_puppiP4_;
0061   const double min_puppi_wgt_;
0062   const bool include_neutrals_;
0063   const bool sort_by_sip2dsig_;
0064   const bool flip_ip_sign_;
0065   const double max_sip3dsig_;
0066   const bool use_hlt_features_;
0067 
0068   edm::EDGetTokenT<edm::View<reco::Jet>> jet_token_;
0069   edm::EDGetTokenT<VertexCollection> vtx_token_;
0070   edm::EDGetTokenT<SVCollection> sv_token_;
0071   edm::EDGetTokenT<CandidateView> pfcand_token_;
0072 
0073   bool use_puppi_value_map_;
0074   bool use_pvasq_value_map_;
0075   bool is_packed_pf_candidate_collection_;
0076 
0077   edm::EDGetTokenT<edm::ValueMap<float>> puppi_value_map_token_;
0078   edm::EDGetTokenT<edm::ValueMap<int>> pvasq_value_map_token_;
0079   edm::EDGetTokenT<edm::Association<VertexCollection>> pvas_token_;
0080   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> track_builder_token_;
0081 
0082   edm::Handle<VertexCollection> vtxs_;
0083   edm::Handle<SVCollection> svs_;
0084   edm::Handle<CandidateView> pfcands_;
0085   edm::ESHandle<TransientTrackBuilder> track_builder_;
0086   edm::Handle<edm::ValueMap<float>> puppi_value_map_;
0087   edm::Handle<edm::ValueMap<int>> pvasq_value_map_;
0088   edm::Handle<edm::Association<VertexCollection>> pvas_;
0089 
0090   const static std::vector<std::string> particle_features_;
0091   const static std::vector<std::string> sv_features_;
0092   const static std::vector<std::string> particle_features_hlt_;
0093   const static std::vector<std::string> sv_features_hlt_;
0094   const reco::Vertex *pv_ = nullptr;
0095   const static float min_track_pt_property_;
0096   const static int min_valid_pixel_hits_;
0097 
0098   std::map<reco::CandidatePtr::key_type, float> puppi_wgt_cache;
0099 };
0100 
0101 const std::vector<std::string> DeepBoostedJetTagInfoProducer::particle_features_{
0102     "pfcand_puppiw",        "pfcand_hcalFrac",       "pfcand_VTX_ass",      "pfcand_lostInnerHits",
0103     "pfcand_quality",       "pfcand_charge",         "pfcand_isEl",         "pfcand_isMu",
0104     "pfcand_isChargedHad",  "pfcand_isGamma",        "pfcand_isNeutralHad", "pfcand_phirel",
0105     "pfcand_etarel",        "pfcand_deltaR",         "pfcand_abseta",       "pfcand_ptrel_log",
0106     "pfcand_erel_log",      "pfcand_pt_log",         "pfcand_drminsv",      "pfcand_drsubjet1",
0107     "pfcand_drsubjet2",     "pfcand_normchi2",       "pfcand_dz",           "pfcand_dzsig",
0108     "pfcand_dxy",           "pfcand_dxysig",         "pfcand_dptdpt",       "pfcand_detadeta",
0109     "pfcand_dphidphi",      "pfcand_dxydxy",         "pfcand_dzdz",         "pfcand_dxydz",
0110     "pfcand_dphidxy",       "pfcand_dlambdadz",      "pfcand_btagEtaRel",   "pfcand_btagPtRatio",
0111     "pfcand_btagPParRatio", "pfcand_btagSip2dVal",   "pfcand_btagSip2dSig", "pfcand_btagSip3dVal",
0112     "pfcand_btagSip3dSig",  "pfcand_btagJetDistVal", "pfcand_mask",         "pfcand_pt_log_nopuppi",
0113     "pfcand_e_log_nopuppi", "pfcand_ptrel",          "pfcand_erel"};
0114 
0115 const std::vector<std::string> DeepBoostedJetTagInfoProducer::particle_features_hlt_{"jet_pfcand_pt_log",
0116                                                                                      "jet_pfcand_energy_log",
0117                                                                                      "jet_pfcand_deta",
0118                                                                                      "jet_pfcand_dphi",
0119                                                                                      "jet_pfcand_eta",
0120                                                                                      "jet_pfcand_charge",
0121                                                                                      "jet_pfcand_frompv",
0122                                                                                      "jet_pfcand_nlostinnerhits",
0123                                                                                      "jet_pfcand_track_chi2",
0124                                                                                      "jet_pfcand_track_qual",
0125                                                                                      "jet_pfcand_dz",
0126                                                                                      "jet_pfcand_dzsig",
0127                                                                                      "jet_pfcand_dxy",
0128                                                                                      "jet_pfcand_dxysig",
0129                                                                                      "jet_pfcand_etarel",
0130                                                                                      "jet_pfcand_pperp_ratio",
0131                                                                                      "jet_pfcand_ppara_ratio",
0132                                                                                      "jet_pfcand_trackjet_d3d",
0133                                                                                      "jet_pfcand_trackjet_d3dsig",
0134                                                                                      "jet_pfcand_trackjet_dist",
0135                                                                                      "jet_pfcand_nhits",
0136                                                                                      "jet_pfcand_npixhits",
0137                                                                                      "jet_pfcand_nstriphits",
0138                                                                                      "jet_pfcand_trackjet_decayL",
0139                                                                                      "jet_pfcand_puppiw",
0140                                                                                      "pfcand_mask"};
0141 
0142 const std::vector<std::string> DeepBoostedJetTagInfoProducer::sv_features_{"sv_mask",
0143                                                                            "sv_ptrel",
0144                                                                            "sv_erel",
0145                                                                            "sv_phirel",
0146                                                                            "sv_etarel",
0147                                                                            "sv_deltaR",
0148                                                                            "sv_abseta",
0149                                                                            "sv_mass",
0150                                                                            "sv_ptrel_log",
0151                                                                            "sv_erel_log",
0152                                                                            "sv_pt_log",
0153                                                                            "sv_pt",
0154                                                                            "sv_ntracks",
0155                                                                            "sv_normchi2",
0156                                                                            "sv_dxy",
0157                                                                            "sv_dxysig",
0158                                                                            "sv_d3d",
0159                                                                            "sv_d3dsig",
0160                                                                            "sv_costhetasvpv"};
0161 
0162 const std::vector<std::string> DeepBoostedJetTagInfoProducer::sv_features_hlt_{"jet_sv_pt_log",
0163                                                                                "jet_sv_mass",
0164                                                                                "jet_sv_deta",
0165                                                                                "jet_sv_dphi",
0166                                                                                "jet_sv_eta",
0167                                                                                "jet_sv_ntrack",
0168                                                                                "jet_sv_chi2",
0169                                                                                "jet_sv_dxy",
0170                                                                                "jet_sv_dxysig",
0171                                                                                "jet_sv_d3d",
0172                                                                                "jet_sv_d3dsig",
0173                                                                                "sv_mask"};
0174 
0175 const float DeepBoostedJetTagInfoProducer::min_track_pt_property_ = 0.5;
0176 const int DeepBoostedJetTagInfoProducer::min_valid_pixel_hits_ = 0;
0177 
0178 DeepBoostedJetTagInfoProducer::DeepBoostedJetTagInfoProducer(const edm::ParameterSet &iConfig)
0179     : jet_radius_(iConfig.getParameter<double>("jet_radius")),
0180       min_jet_pt_(iConfig.getParameter<double>("min_jet_pt")),
0181       max_jet_eta_(iConfig.getParameter<double>("max_jet_eta")),
0182       min_pt_for_track_properties_(iConfig.getParameter<double>("min_pt_for_track_properties")),
0183       min_pt_for_pfcandidates_(iConfig.getParameter<double>("min_pt_for_pfcandidates")),
0184       use_puppiP4_(iConfig.getParameter<bool>("use_puppiP4")),
0185       min_puppi_wgt_(iConfig.getParameter<double>("min_puppi_wgt")),
0186       include_neutrals_(iConfig.getParameter<bool>("include_neutrals")),
0187       sort_by_sip2dsig_(iConfig.getParameter<bool>("sort_by_sip2dsig")),
0188       flip_ip_sign_(iConfig.getParameter<bool>("flip_ip_sign")),
0189       max_sip3dsig_(iConfig.getParameter<double>("sip3dSigMax")),
0190       use_hlt_features_(iConfig.getParameter<bool>("use_hlt_features")),
0191       jet_token_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"))),
0192       vtx_token_(consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0193       sv_token_(consumes<SVCollection>(iConfig.getParameter<edm::InputTag>("secondary_vertices"))),
0194       pfcand_token_(consumes<CandidateView>(iConfig.getParameter<edm::InputTag>("pf_candidates"))),
0195       use_puppi_value_map_(false),
0196       use_pvasq_value_map_(false),
0197       track_builder_token_(
0198           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))) {
0199   const auto &puppi_value_map_tag = iConfig.getParameter<edm::InputTag>("puppi_value_map");
0200   if (!puppi_value_map_tag.label().empty()) {
0201     puppi_value_map_token_ = consumes<edm::ValueMap<float>>(puppi_value_map_tag);
0202     use_puppi_value_map_ = true;
0203   }
0204 
0205   const auto &pvas_tag = iConfig.getParameter<edm::InputTag>("vertex_associator");
0206   if (!pvas_tag.label().empty()) {
0207     pvasq_value_map_token_ = consumes<edm::ValueMap<int>>(pvas_tag);
0208     pvas_token_ = consumes<edm::Association<VertexCollection>>(pvas_tag);
0209     use_pvasq_value_map_ = true;
0210   }
0211 
0212   produces<DeepBoostedJetTagInfoCollection>();
0213 }
0214 
0215 DeepBoostedJetTagInfoProducer::~DeepBoostedJetTagInfoProducer() {}
0216 
0217 void DeepBoostedJetTagInfoProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0218   // pfDeepBoostedJetTagInfos
0219   edm::ParameterSetDescription desc;
0220   desc.add<double>("jet_radius", 0.8);
0221   desc.add<double>("min_jet_pt", 150);
0222   desc.add<double>("max_jet_eta", 99);
0223   desc.add<double>("min_pt_for_track_properties", -1);
0224   desc.add<double>("min_pt_for_pfcandidates", -1);
0225   desc.add<bool>("use_puppiP4", true);
0226   desc.add<bool>("include_neutrals", true);
0227   desc.add<bool>("sort_by_sip2dsig", false);
0228   desc.add<double>("min_puppi_wgt", 0.01);
0229   desc.add<bool>("flip_ip_sign", false);
0230   desc.add<double>("sip3dSigMax", -1);
0231   desc.add<bool>("use_hlt_features", false);
0232   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0233   desc.add<edm::InputTag>("secondary_vertices", edm::InputTag("inclusiveCandidateSecondaryVertices"));
0234   desc.add<edm::InputTag>("pf_candidates", edm::InputTag("particleFlow"));
0235   desc.add<edm::InputTag>("jets", edm::InputTag("ak8PFJetsPuppi"));
0236   desc.add<edm::InputTag>("puppi_value_map", edm::InputTag("puppi"));
0237   desc.add<edm::InputTag>("vertex_associator", edm::InputTag("primaryVertexAssociation", "original"));
0238   descriptions.add("pfDeepBoostedJetTagInfos", desc);
0239 }
0240 
0241 void DeepBoostedJetTagInfoProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0242   // output collection
0243   auto output_tag_infos = std::make_unique<DeepBoostedJetTagInfoCollection>();
0244   // Input jets
0245   auto jets = iEvent.getHandle(jet_token_);
0246   // Primary vertexes
0247   iEvent.getByToken(vtx_token_, vtxs_);
0248   if (vtxs_->empty()) {
0249     // produce empty TagInfos in case no primary vertex
0250     iEvent.put(std::move(output_tag_infos));
0251     return;  // exit event
0252   }
0253   // Leading vertex
0254   pv_ = &vtxs_->at(0);
0255   // Secondary vertexs
0256   iEvent.getByToken(sv_token_, svs_);
0257   // PF candidates
0258   iEvent.getByToken(pfcand_token_, pfcands_);
0259   is_packed_pf_candidate_collection_ = false;
0260   if (dynamic_cast<const pat::PackedCandidateCollection *>(pfcands_.product()))
0261     is_packed_pf_candidate_collection_ = true;
0262   // Track builder
0263   track_builder_ = iSetup.getHandle(track_builder_token_);
0264   // Puppi weight value map
0265   if (use_puppi_value_map_)
0266     iEvent.getByToken(puppi_value_map_token_, puppi_value_map_);
0267   // Vertex associator map
0268   if (use_pvasq_value_map_) {
0269     iEvent.getByToken(pvasq_value_map_token_, pvasq_value_map_);
0270     iEvent.getByToken(pvas_token_, pvas_);
0271   }
0272 
0273   // Loop over jet
0274   for (std::size_t jet_n = 0; jet_n < jets->size(); jet_n++) {
0275     const auto &jet = (*jets)[jet_n];
0276     edm::RefToBase<reco::Jet> jet_ref(jets, jet_n);
0277 
0278     // create jet features
0279     DeepBoostedJetFeatures features;
0280     if (not use_hlt_features_) {
0281       for (const auto &name : particle_features_) {
0282         features.add(name);
0283       }
0284       for (const auto &name : sv_features_) {
0285         features.add(name);
0286       }
0287     } else {
0288       // declare all the feature variables (init as empty vector)
0289       for (const auto &name : particle_features_hlt_) {
0290         features.add(name);
0291       }
0292       for (const auto &name : sv_features_hlt_) {
0293         features.add(name);
0294       }
0295     }
0296 
0297     // fill values only if above pt threshold and has daughters, otherwise left
0298     bool fill_vars = true;
0299     if (jet.pt() < min_jet_pt_ or std::abs(jet.eta()) > max_jet_eta_)
0300       fill_vars = false;
0301     if (jet.numberOfDaughters() == 0)
0302       fill_vars = false;
0303 
0304     // fill features
0305     if (fill_vars) {
0306       fillParticleFeatures(features, jet);
0307       fillSVFeatures(features, jet);
0308       if (use_hlt_features_) {
0309         features.check_consistency(particle_features_hlt_);
0310         features.check_consistency(sv_features_hlt_);
0311       } else {
0312         features.check_consistency(particle_features_);
0313         features.check_consistency(sv_features_);
0314       }
0315     }
0316     // this should always be done even if features are not filled
0317     output_tag_infos->emplace_back(features, jet_ref);
0318   }
0319   // move output collection
0320   iEvent.put(std::move(output_tag_infos));
0321 }
0322 
0323 float DeepBoostedJetTagInfoProducer::puppiWgt(const reco::CandidatePtr &cand) {
0324   const auto *pack_cand = dynamic_cast<const pat::PackedCandidate *>(&(*cand));
0325   const auto *reco_cand = dynamic_cast<const reco::PFCandidate *>(&(*cand));
0326   float wgt = 1.;
0327   if (pack_cand)
0328     wgt = pack_cand->puppiWeight();
0329   else if (reco_cand) {
0330     if (use_puppi_value_map_)
0331       wgt = (*puppi_value_map_)[cand];
0332   } else
0333     throw edm::Exception(edm::errors::InvalidReference)
0334         << "Cannot convert to either pat::PackedCandidate or reco::PFCandidate";
0335   puppi_wgt_cache[cand.key()] = wgt;
0336   return wgt;
0337 }
0338 
0339 bool DeepBoostedJetTagInfoProducer::useTrackProperties(const reco::PFCandidate *reco_cand) {
0340   const auto *track = reco_cand->bestTrack();
0341   return track != nullptr and track->pt() > min_pt_for_track_properties_;
0342 };
0343 
0344 void DeepBoostedJetTagInfoProducer::fillParticleFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet) {
0345   // some jet properties
0346   math::XYZVector jet_dir = jet.momentum().Unit();
0347   TVector3 jet_direction(jet.momentum().Unit().x(), jet.momentum().Unit().y(), jet.momentum().Unit().z());
0348   GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz());
0349   const float etasign = jet.eta() > 0 ? 1 : -1;
0350   // vertexes
0351   reco::VertexRefProd PVRefProd(vtxs_);
0352   // track builder
0353   TrackInfoBuilder trackinfo(track_builder_);
0354 
0355   // make list of pf-candidates to be considered
0356   std::vector<reco::CandidatePtr> daughters;
0357   for (const auto &dau : jet.daughterPtrVector()) {
0358     // remove particles w/ extremely low puppi weights
0359     // [Note] use jet daughters here to get the puppiWgt correctly
0360     if ((puppiWgt(dau)) < min_puppi_wgt_)
0361       continue;
0362     // from here: get the original reco/packed candidate not scaled by the puppi weight
0363     auto cand = pfcands_->ptrAt(dau.key());
0364     // base requirements on PF candidates
0365     if (use_hlt_features_ and cand->pt() < min_pt_for_pfcandidates_)
0366       continue;
0367     // charged candidate selection (for Higgs Interaction Net)
0368     if (!include_neutrals_ and (cand->charge() == 0 or cand->pt() < min_pt_for_track_properties_))
0369       continue;
0370     // only when computing the nagative tagger: remove charged candidates with high sip3d
0371     if (flip_ip_sign_ and cand->charge()) {
0372       trackinfo.buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_);
0373       if (trackinfo.getTrackSip3dSig() > max_sip3dsig_)
0374         continue;
0375     }
0376     // filling daughters
0377     daughters.push_back(cand);
0378   }
0379 
0380   // Sorting of PF-candidates
0381   std::vector<btagbtvdeep::SortingClass<reco::CandidatePtr>> c_sorted;
0382   if (sort_by_sip2dsig_) {
0383     // sort charged pf candidates by 2d impact parameter significance
0384     for (const auto &cand : daughters) {
0385       trackinfo.buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_);
0386       c_sorted.emplace_back(cand,
0387                             trackinfo.getTrackSip2dSig(),
0388                             -btagbtvdeep::mindrsvpfcand(*svs_, &(*cand), jet_radius_),
0389                             cand->pt() / jet.pt());
0390     }
0391     std::sort(c_sorted.begin(), c_sorted.end(), btagbtvdeep::SortingClass<reco::CandidatePtr>::compareByABCInv);
0392     for (unsigned int i = 0; i < c_sorted.size(); i++) {
0393       const auto &c = c_sorted.at(i);
0394       const auto &cand = c.get();
0395       daughters.at(i) = cand;
0396     }
0397   } else {
0398     // sort by Puppi-weighted pt
0399     if (use_puppiP4_)
0400       std::sort(daughters.begin(), daughters.end(), [&](const reco::CandidatePtr &a, const reco::CandidatePtr &b) {
0401         return puppi_wgt_cache.at(a.key()) * a->pt() > puppi_wgt_cache.at(b.key()) * b->pt();
0402       });
0403     // sort by original pt (not Puppi-weighted)
0404     else
0405       std::sort(daughters.begin(), daughters.end(), [](const auto &a, const auto &b) { return a->pt() > b->pt(); });
0406   }
0407 
0408   // reserve space
0409   if (use_hlt_features_) {
0410     for (const auto &name : particle_features_hlt_)
0411       fts.reserve(name, daughters.size());
0412   } else {
0413     for (const auto &name : particle_features_)
0414       fts.reserve(name, daughters.size());
0415   }
0416 
0417   // build white list of candidates i.e. particles and tracks belonging to SV. Needed when input particles are not packed PF candidates
0418   std::vector<unsigned int> whiteListSV;
0419   std::vector<reco::TrackRef> whiteListTk;
0420   if (not is_packed_pf_candidate_collection_) {
0421     for (size_t isv = 0; isv < svs_->size(); isv++) {
0422       for (size_t icand = 0; icand < svs_->at(isv).numberOfSourceCandidatePtrs(); icand++) {
0423         const edm::Ptr<reco::Candidate> &cand = svs_->at(isv).sourceCandidatePtr(icand);
0424         if (cand.id() == pfcands_.id())
0425           whiteListSV.push_back(cand.key());
0426       }
0427       for (auto cand = svs_->at(isv).begin(); cand != svs_->at(isv).end(); cand++) {
0428         const reco::RecoChargedCandidate *chCand = dynamic_cast<const reco::RecoChargedCandidate *>(&(*cand));
0429         if (chCand != nullptr) {
0430           whiteListTk.push_back(chCand->track());
0431         }
0432       }
0433     }
0434   }
0435 
0436   // Build observables
0437   size_t icand = 0;
0438   for (const auto &cand : daughters) {
0439     const auto *packed_cand = dynamic_cast<const pat::PackedCandidate *>(&(*cand));
0440     const auto *reco_cand = dynamic_cast<const reco::PFCandidate *>(&(*cand));
0441 
0442     if (not packed_cand and not reco_cand)
0443       throw edm::Exception(edm::errors::InvalidReference)
0444           << "Cannot convert to either reco::PFCandidate or pat::PackedCandidate";
0445 
0446     if (!include_neutrals_ and
0447         ((packed_cand and !packed_cand->hasTrackDetails()) or (reco_cand and !useTrackProperties(reco_cand)))) {
0448       icand++;
0449       continue;
0450     }
0451 
0452     const float ip_sign = flip_ip_sign_ ? -1 : 1;
0453 
0454     // input particle is a packed PF candidate
0455     auto candP4 = use_puppiP4_ ? puppi_wgt_cache.at(cand.key()) * cand->p4() : cand->p4();
0456     auto candP3 = use_puppiP4_ ? puppi_wgt_cache.at(cand.key()) * cand->momentum() : cand->momentum();
0457 
0458     // candidate track
0459     const reco::Track *track = nullptr;
0460     if (packed_cand)
0461       track = packed_cand->bestTrack();
0462     else if (reco_cand and useTrackProperties(reco_cand))
0463       track = reco_cand->bestTrack();
0464 
0465     // reco-vertex association
0466     int pv_ass_quality = 0;
0467     reco::VertexRef pv_ass;
0468     float vtx_ass = 0;
0469     if (reco_cand) {
0470       if (use_pvasq_value_map_) {
0471         pv_ass_quality = (*pvasq_value_map_)[cand];
0472         pv_ass = (*pvas_)[cand];
0473         vtx_ass = vtx_ass_from_pfcand(*reco_cand, pv_ass_quality, pv_ass);
0474       } else
0475         throw edm::Exception(edm::errors::InvalidReference) << "Vertex association missing";
0476     }
0477 
0478     // Building offline features
0479     if (not use_hlt_features_) {
0480       // in case of packed candidate
0481       if (packed_cand) {
0482         float hcal_fraction = 0.;
0483         if (packed_cand->pdgId() == 1 or packed_cand->pdgId() == 130)
0484           hcal_fraction = packed_cand->hcalFraction();
0485         else if (packed_cand->isIsolatedChargedHadron())
0486           hcal_fraction = packed_cand->rawHcalFraction();
0487 
0488         fts.fill("pfcand_hcalFrac", hcal_fraction);
0489         fts.fill("pfcand_VTX_ass", packed_cand->pvAssociationQuality());
0490         fts.fill("pfcand_lostInnerHits", packed_cand->lostInnerHits());
0491         fts.fill("pfcand_quality", track ? track->qualityMask() : 0);
0492         fts.fill("pfcand_charge", packed_cand->charge());
0493         fts.fill("pfcand_isEl", std::abs(packed_cand->pdgId()) == 11);
0494         fts.fill("pfcand_isMu", std::abs(packed_cand->pdgId()) == 13);
0495         fts.fill("pfcand_isChargedHad", std::abs(packed_cand->pdgId()) == 211);
0496         fts.fill("pfcand_isGamma", std::abs(packed_cand->pdgId()) == 22);
0497         fts.fill("pfcand_isNeutralHad", std::abs(packed_cand->pdgId()) == 130);
0498         fts.fill("pfcand_dz", ip_sign * packed_cand->dz());
0499         fts.fill("pfcand_dxy", ip_sign * packed_cand->dxy());
0500         fts.fill("pfcand_dzsig", track ? ip_sign * packed_cand->dz() / packed_cand->dzError() : 0);
0501         fts.fill("pfcand_dxysig", track ? ip_sign * packed_cand->dxy() / packed_cand->dxyError() : 0);
0502 
0503       }
0504       // in the case of reco candidate
0505       else if (reco_cand) {
0506         fts.fill("pfcand_hcalFrac", reco_cand->hcalEnergy() / (reco_cand->ecalEnergy() + reco_cand->hcalEnergy()));
0507         fts.fill("pfcand_VTX_ass", vtx_ass);
0508         fts.fill("pfcand_lostInnerHits", useTrackProperties(reco_cand) ? lost_inner_hits_from_pfcand(*reco_cand) : 0);
0509         fts.fill("pfcand_quality", useTrackProperties(reco_cand) ? quality_from_pfcand(*reco_cand) : 0);
0510         fts.fill("pfcand_charge", reco_cand->charge());
0511         fts.fill("pfcand_isEl", std::abs(reco_cand->pdgId()) == 11);
0512         fts.fill("pfcand_isMu", std::abs(reco_cand->pdgId()) == 13);
0513         fts.fill("pfcand_isChargedHad", std::abs(reco_cand->pdgId()) == 211);
0514         fts.fill("pfcand_isGamma", std::abs(reco_cand->pdgId()) == 22);
0515         fts.fill("pfcand_isNeutralHad", std::abs(reco_cand->pdgId()) == 130);
0516         fts.fill("pfcand_dz", track ? ip_sign * track->dz(pv_->position()) : 0);
0517         fts.fill("pfcand_dzsig", track ? ip_sign * track->dz(pv_->position()) / track->dzError() : 0);
0518         fts.fill("pfcand_dxy", track ? ip_sign * track->dxy(pv_->position()) : 0);
0519         fts.fill("pfcand_dxysig", track ? ip_sign * track->dxy(pv_->position()) / track->dxyError() : 0);
0520       }
0521 
0522       // generic candidate observables
0523       fts.fill("pfcand_puppiw", puppi_wgt_cache.at(cand.key()));
0524       fts.fill("pfcand_phirel", reco::deltaPhi(candP4, jet));
0525       fts.fill("pfcand_etarel", etasign * (candP4.eta() - jet.eta()));
0526       fts.fill("pfcand_deltaR", reco::deltaR(candP4, jet));
0527       fts.fill("pfcand_abseta", std::abs(candP4.eta()));
0528 
0529       fts.fill("pfcand_ptrel_log", std::log(candP4.pt() / jet.pt()));
0530       fts.fill("pfcand_ptrel", candP4.pt() / jet.pt());
0531       fts.fill("pfcand_erel_log", std::log(candP4.energy() / jet.energy()));
0532       fts.fill("pfcand_erel", candP4.energy() / jet.energy());
0533       fts.fill("pfcand_pt_log", std::log(candP4.pt()));
0534 
0535       fts.fill("pfcand_mask", 1);
0536       fts.fill("pfcand_pt_log_nopuppi", std::log(cand->pt()));
0537       fts.fill("pfcand_e_log_nopuppi", std::log(cand->energy()));
0538 
0539       float drminpfcandsv = btagbtvdeep::mindrsvpfcand(*svs_, &(*cand), std::numeric_limits<float>::infinity());
0540       fts.fill("pfcand_drminsv", drminpfcandsv);
0541 
0542       if (track) {
0543         auto cov = [&](unsigned i, unsigned j) { return track->covariance(i, j); };
0544         fts.fill("pfcand_dptdpt", cov(0, 0));
0545         fts.fill("pfcand_detadeta", cov(1, 1));
0546         fts.fill("pfcand_dphidphi", cov(2, 2));
0547         fts.fill("pfcand_dxydxy", cov(3, 3));
0548         fts.fill("pfcand_dzdz", cov(4, 4));
0549         fts.fill("pfcand_dxydz", cov(3, 4));
0550         fts.fill("pfcand_dphidxy", cov(2, 3));
0551         fts.fill("pfcand_dlambdadz", cov(1, 4));
0552 
0553         fts.fill("pfcand_normchi2", std::floor(track->normalizedChi2()));
0554 
0555         trackinfo.buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_);
0556         fts.fill("pfcand_btagEtaRel", trackinfo.getTrackEtaRel());
0557         fts.fill("pfcand_btagPtRatio", trackinfo.getTrackPtRatio());
0558         fts.fill("pfcand_btagPParRatio", trackinfo.getTrackPParRatio());
0559         fts.fill("pfcand_btagSip2dVal", ip_sign * trackinfo.getTrackSip2dVal());
0560         fts.fill("pfcand_btagSip2dSig", ip_sign * trackinfo.getTrackSip2dSig());
0561         fts.fill("pfcand_btagSip3dVal", ip_sign * trackinfo.getTrackSip3dVal());
0562         fts.fill("pfcand_btagSip3dSig", ip_sign * trackinfo.getTrackSip3dSig());
0563         fts.fill("pfcand_btagJetDistVal", trackinfo.getTrackJetDistVal());
0564       } else {
0565         fts.fill("pfcand_normchi2", 999);
0566         fts.fill("pfcand_dptdpt", 0);
0567         fts.fill("pfcand_detadeta", 0);
0568         fts.fill("pfcand_dphidphi", 0);
0569         fts.fill("pfcand_dxydxy", 0);
0570         fts.fill("pfcand_dzdz", 0);
0571         fts.fill("pfcand_dxydz", 0);
0572         fts.fill("pfcand_dphidxy", 0);
0573         fts.fill("pfcand_dlambdadz", 0);
0574         fts.fill("pfcand_btagEtaRel", 0);
0575         fts.fill("pfcand_btagPtRatio", 0);
0576         fts.fill("pfcand_btagPParRatio", 0);
0577         fts.fill("pfcand_btagSip2dVal", 0);
0578         fts.fill("pfcand_btagSip2dSig", 0);
0579         fts.fill("pfcand_btagSip3dVal", 0);
0580         fts.fill("pfcand_btagSip3dSig", 0);
0581         fts.fill("pfcand_btagJetDistVal", 0);
0582       }
0583 
0584       // subjets only if the incomming jets is a PAT one
0585       const auto *patJet = dynamic_cast<const pat::Jet *>(&jet);
0586       if (patJet and patJet->nSubjetCollections() > 0) {
0587         auto subjets = patJet->subjets();
0588         std::sort(subjets.begin(), subjets.end(), [](const edm::Ptr<pat::Jet> &p1, const edm::Ptr<pat::Jet> &p2) {
0589           return p1->pt() > p2->pt();
0590         });
0591         fts.fill("pfcand_drsubjet1", !subjets.empty() ? reco::deltaR(*cand, *subjets.at(0)) : -1);
0592         fts.fill("pfcand_drsubjet2", subjets.size() > 1 ? reco::deltaR(*cand, *subjets.at(1)) : -1);
0593       } else {
0594         fts.fill("pfcand_drsubjet1", -1);
0595         fts.fill("pfcand_drsubjet2", -1);
0596       }
0597     }
0598     // using HLT features
0599     else {
0600       pat::PackedCandidate candidate;
0601       math::XYZPoint pv_ass_pos;
0602       // In case input is a packed candidate (evaluate HLT network on offline)
0603       if (packed_cand) {
0604         candidate = *packed_cand;
0605         pv_ass = reco::VertexRef(vtxs_, 0);
0606         pv_ass_pos = pv_ass->position();
0607       }
0608       // In case input is a reco::PFCandidate
0609       else if (reco_cand) {
0610         // follow what is done in PhysicsTools/PatAlgos/plugins/PATPackedCandidateProducer.cc to minimize differences between HLT and RECO choices of input observables
0611         // no reference to vertex, take closest in dz or position 0
0612         if (not pv_ass.isNonnull()) {
0613           if (track) {
0614             float z_dist = 99999;
0615             int pv_pos = -1;
0616             for (size_t iv = 0; iv < vtxs_->size(); iv++) {
0617               float dz = std::abs(track->dz(((*vtxs_)[iv]).position()));
0618               if (dz < z_dist) {
0619                 z_dist = dz;
0620                 pv_pos = iv;
0621               }
0622             }
0623             pv_ass = reco::VertexRef(vtxs_, pv_pos);
0624           } else
0625             pv_ass = reco::VertexRef(vtxs_, 0);
0626         }
0627         pv_ass_pos = pv_ass->position();
0628 
0629         // create a transient packed candidate
0630         if (track) {
0631           candidate = pat::PackedCandidate(cand->polarP4(),
0632                                            track->referencePoint(),
0633                                            track->pt(),
0634                                            track->eta(),
0635                                            track->phi(),
0636                                            cand->pdgId(),
0637                                            PVRefProd,
0638                                            pv_ass.key());
0639           candidate.setAssociationQuality(pat::PackedCandidate::PVAssociationQuality(
0640               btagbtvdeep::vtx_ass_from_pfcand(*reco_cand, pv_ass_quality, pv_ass)));
0641           candidate.setCovarianceVersion(0);
0642           pat::PackedCandidate::LostInnerHits lostHits = pat::PackedCandidate::noLostInnerHits;
0643           int nlost = track->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0644           if (nlost == 0) {
0645             if (track->hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1))
0646               lostHits = pat::PackedCandidate::validHitInFirstPixelBarrelLayer;
0647           } else
0648             lostHits = (nlost == 1 ? pat::PackedCandidate::oneLostInnerHit : pat::PackedCandidate::moreLostInnerHits);
0649           candidate.setLostInnerHits(lostHits);
0650 
0651           if (useTrackProperties(reco_cand) or
0652               std::find(whiteListSV.begin(), whiteListSV.end(), icand) != whiteListSV.end() or
0653               std::find(whiteListTk.begin(), whiteListTk.end(), reco_cand->trackRef()) != whiteListTk.end()) {
0654             candidate.setFirstHit(track->hitPattern().getHitPattern(reco::HitPattern::TRACK_HITS, 0));
0655             if (abs(cand->pdgId()) == 22)
0656               candidate.setTrackProperties(*track, 0, 0);
0657             else {
0658               if (track->hitPattern().numberOfValidPixelHits() > min_valid_pixel_hits_)
0659                 candidate.setTrackProperties(*track, 8, 0);
0660               else
0661                 candidate.setTrackProperties(*track, 264, 0);
0662             }
0663           } else {
0664             if (candidate.pt() > min_track_pt_property_) {
0665               if (track->hitPattern().numberOfValidPixelHits() > 0)
0666                 candidate.setTrackProperties(*track, 520, 0);
0667               else
0668                 candidate.setTrackProperties(*track, 776, 0);
0669             }
0670           }
0671           candidate.setTrackHighPurity(reco_cand->trackRef().isNonnull() and
0672                                        reco_cand->trackRef()->quality(reco::Track::highPurity));
0673         } else {
0674           candidate = pat::PackedCandidate(
0675               cand->polarP4(), pv_ass_pos, cand->pt(), cand->eta(), cand->phi(), cand->pdgId(), PVRefProd, pv_ass.key());
0676           candidate.setAssociationQuality(
0677               pat::PackedCandidate::PVAssociationQuality(pat::PackedCandidate::UsedInFitTight));
0678         }
0679         /// override track
0680         track = candidate.bestTrack();
0681       }
0682 
0683       TVector3 cand_direction(candP3.x(), candP3.y(), candP3.z());
0684 
0685       fts.fill("jet_pfcand_pt_log", std::log(candP4.pt()));
0686       fts.fill("jet_pfcand_energy_log", std::log(candP4.energy()));
0687       fts.fill("jet_pfcand_eta", candP4.eta());
0688       fts.fill("jet_pfcand_deta", jet_direction.Eta() - cand_direction.Eta());
0689       fts.fill("jet_pfcand_dphi", jet_direction.DeltaPhi(cand_direction));
0690       fts.fill("jet_pfcand_charge", cand->charge());
0691       fts.fill("jet_pfcand_etarel", reco::btau::etaRel(jet_dir, candP3));
0692       fts.fill("jet_pfcand_pperp_ratio", jet_direction.Perp(cand_direction) / cand_direction.Mag());
0693       fts.fill("jet_pfcand_ppara_ratio", jet_direction.Dot(cand_direction) / cand_direction.Mag());
0694       fts.fill("jet_pfcand_frompv", candidate.fromPV());
0695       fts.fill("jet_pfcand_dz", candidate.dz(pv_ass_pos));
0696       fts.fill("jet_pfcand_dxy", candidate.dxy(pv_ass_pos));
0697       fts.fill("jet_pfcand_puppiw", puppi_wgt_cache.at(cand.key()));
0698       fts.fill("jet_pfcand_nlostinnerhits", candidate.lostInnerHits());
0699       fts.fill("jet_pfcand_nhits", candidate.numberOfHits());
0700       fts.fill("jet_pfcand_npixhits", candidate.numberOfPixelHits());
0701       fts.fill("jet_pfcand_nstriphits", candidate.stripLayersWithMeasurement());
0702       fts.fill("pfcand_mask", 1);
0703 
0704       if (track) {
0705         fts.fill("jet_pfcand_dzsig", fabs(candidate.dz(pv_ass_pos)) / candidate.dzError());
0706         fts.fill("jet_pfcand_dxysig", fabs(candidate.dxy(pv_ass_pos)) / candidate.dxyError());
0707         fts.fill("jet_pfcand_track_chi2", track->normalizedChi2());
0708         fts.fill("jet_pfcand_track_qual", track->qualityMask());
0709 
0710         reco::TransientTrack transientTrack = track_builder_->build(*track);
0711         Measurement1D meas_ip2d =
0712             IPTools::signedTransverseImpactParameter(transientTrack, jet_ref_track_dir, *pv_).second;
0713         Measurement1D meas_ip3d = IPTools::signedImpactParameter3D(transientTrack, jet_ref_track_dir, *pv_).second;
0714         Measurement1D meas_jetdist = IPTools::jetTrackDistance(transientTrack, jet_ref_track_dir, *pv_).second;
0715         Measurement1D meas_decayl = IPTools::signedDecayLength3D(transientTrack, jet_ref_track_dir, *pv_).second;
0716 
0717         fts.fill("jet_pfcand_trackjet_d3d", meas_ip3d.value());
0718         fts.fill("jet_pfcand_trackjet_d3dsig", fabs(meas_ip3d.significance()));
0719         fts.fill("jet_pfcand_trackjet_dist", -meas_jetdist.value());
0720         fts.fill("jet_pfcand_trackjet_decayL", meas_decayl.value());
0721       } else {
0722         fts.fill("jet_pfcand_dzsig", 0);
0723         fts.fill("jet_pfcand_dxysig", 0);
0724         fts.fill("jet_pfcand_track_chi2", 0);
0725         fts.fill("jet_pfcand_track_qual", 0);
0726         fts.fill("jet_pfcand_trackjet_d3d", 0);
0727         fts.fill("jet_pfcand_trackjet_d3dsig", 0);
0728         fts.fill("jet_pfcand_trackjet_dist", 0);
0729         fts.fill("jet_pfcand_trackjet_decayL", 0);
0730       }
0731     }
0732     icand++;
0733   }
0734 }
0735 
0736 void DeepBoostedJetTagInfoProducer::fillSVFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet) {
0737   // secondary vertexes matching jet
0738   std::vector<const reco::VertexCompositePtrCandidate *> jetSVs;
0739   for (const auto &sv : *svs_) {
0740     if (reco::deltaR2(sv, jet) < jet_radius_ * jet_radius_) {
0741       jetSVs.push_back(&sv);
0742     }
0743   }
0744   // sort by dxy significance
0745   std::sort(jetSVs.begin(),
0746             jetSVs.end(),
0747             [&](const reco::VertexCompositePtrCandidate *sva, const reco::VertexCompositePtrCandidate *svb) {
0748               return sv_vertex_comparator(*sva, *svb, *pv_);
0749             });
0750 
0751   // reserve space
0752   if (not use_hlt_features_) {
0753     for (const auto &name : sv_features_) {
0754       fts.reserve(name, jetSVs.size());
0755     }
0756   } else {
0757     for (const auto &name : sv_features_hlt_) {
0758       fts.reserve(name, jetSVs.size());
0759     }
0760   }
0761 
0762   const float etasign = jet.eta() > 0 ? 1 : -1;
0763 
0764   GlobalVector jet_global_vec(jet.px(), jet.py(), jet.pz());
0765 
0766   for (const auto *sv : jetSVs) {
0767     // features for reco
0768     if (not use_hlt_features_) {
0769       fts.fill("sv_mask", 1);
0770       fts.fill("sv_phirel", reco::deltaPhi(*sv, jet));
0771       fts.fill("sv_etarel", etasign * (sv->eta() - jet.eta()));
0772       fts.fill("sv_deltaR", reco::deltaR(*sv, jet));
0773       fts.fill("sv_abseta", std::abs(sv->eta()));
0774       fts.fill("sv_mass", sv->mass());
0775 
0776       fts.fill("sv_ptrel_log", std::log(sv->pt() / jet.pt()));
0777       fts.fill("sv_ptrel", sv->pt() / jet.pt());
0778       fts.fill("sv_erel_log", std::log(sv->energy() / jet.energy()));
0779       fts.fill("sv_erel", sv->energy() / jet.energy());
0780       fts.fill("sv_pt_log", std::log(sv->pt()));
0781       fts.fill("sv_pt", sv->pt());
0782 
0783       fts.fill("sv_ntracks", sv->numberOfDaughters());
0784       fts.fill("sv_normchi2", sv->vertexNormalizedChi2());
0785       const auto &dxy = vertexDxy(*sv, *pv_);
0786       fts.fill("sv_dxy", dxy.value());
0787       fts.fill("sv_dxysig", dxy.significance());
0788       const auto &d3d = vertexD3d(*sv, *pv_);
0789       fts.fill("sv_d3d", d3d.value());
0790       fts.fill("sv_d3dsig", d3d.significance());
0791 
0792       fts.fill("sv_costhetasvpv", (flip_ip_sign_ ? -1.f : 1.f) * vertexDdotP(*sv, *pv_));
0793     } else {
0794       fts.fill("sv_mask", 1);
0795       fts.fill("jet_sv_pt_log", log(sv->pt()));
0796       fts.fill("jet_sv_eta", sv->eta());
0797       fts.fill("jet_sv_mass", sv->mass());
0798       fts.fill("jet_sv_deta", sv->eta() - jet.eta());
0799       fts.fill("jet_sv_dphi", sv->phi() - jet.phi());
0800       fts.fill("jet_sv_ntrack", sv->numberOfDaughters());
0801       fts.fill("jet_sv_chi2", sv->vertexNormalizedChi2());
0802 
0803       reco::Vertex::CovarianceMatrix csv;
0804       sv->fillVertexCovariance(csv);
0805       reco::Vertex svtx(sv->vertex(), csv);
0806 
0807       VertexDistanceXY dxy;
0808       auto valxy = dxy.signedDistance(svtx, *pv_, jet_global_vec);
0809       fts.fill("jet_sv_dxy", valxy.value());
0810       fts.fill("jet_sv_dxysig", fabs(valxy.significance()));
0811 
0812       VertexDistance3D d3d;
0813       auto val3d = d3d.signedDistance(svtx, *pv_, jet_global_vec);
0814       fts.fill("jet_sv_d3d", val3d.value());
0815       fts.fill("jet_sv_d3dsig", fabs(val3d.significance()));
0816     }
0817   }
0818 }
0819 
0820 // define this as a plug-in
0821 DEFINE_FWK_MODULE(DeepBoostedJetTagInfoProducer);