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
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
0243 auto output_tag_infos = std::make_unique<DeepBoostedJetTagInfoCollection>();
0244
0245 auto jets = iEvent.getHandle(jet_token_);
0246
0247 iEvent.getByToken(vtx_token_, vtxs_);
0248 if (vtxs_->empty()) {
0249
0250 iEvent.put(std::move(output_tag_infos));
0251 return;
0252 }
0253
0254 pv_ = &vtxs_->at(0);
0255
0256 iEvent.getByToken(sv_token_, svs_);
0257
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
0263 track_builder_ = iSetup.getHandle(track_builder_token_);
0264
0265 if (use_puppi_value_map_)
0266 iEvent.getByToken(puppi_value_map_token_, puppi_value_map_);
0267
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
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
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
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
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
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
0317 output_tag_infos->emplace_back(features, jet_ref);
0318 }
0319
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
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
0351 reco::VertexRefProd PVRefProd(vtxs_);
0352
0353 TrackInfoBuilder trackinfo(track_builder_);
0354
0355
0356 std::vector<reco::CandidatePtr> daughters;
0357 for (const auto &dau : jet.daughterPtrVector()) {
0358
0359
0360 if ((puppiWgt(dau)) < min_puppi_wgt_)
0361 continue;
0362
0363 auto cand = pfcands_->ptrAt(dau.key());
0364
0365 if (use_hlt_features_ and cand->pt() < min_pt_for_pfcandidates_)
0366 continue;
0367
0368 if (!include_neutrals_ and (cand->charge() == 0 or cand->pt() < min_pt_for_track_properties_))
0369 continue;
0370
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
0377 daughters.push_back(cand);
0378 }
0379
0380
0381 std::vector<btagbtvdeep::SortingClass<reco::CandidatePtr>> c_sorted;
0382 if (sort_by_sip2dsig_) {
0383
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
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
0404 else
0405 std::sort(daughters.begin(), daughters.end(), [](const auto &a, const auto &b) { return a->pt() > b->pt(); });
0406 }
0407
0408
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
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
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
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
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
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
0479 if (not use_hlt_features_) {
0480
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
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
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
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
0599 else {
0600 pat::PackedCandidate candidate;
0601 math::XYZPoint pv_ass_pos;
0602
0603 if (packed_cand) {
0604 candidate = *packed_cand;
0605 pv_ass = reco::VertexRef(vtxs_, 0);
0606 pv_ass_pos = pv_ass->position();
0607 }
0608
0609 else if (reco_cand) {
0610
0611
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
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
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
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
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
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
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
0821 DEFINE_FWK_MODULE(DeepBoostedJetTagInfoProducer);
