Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-11-05 05:20:43 

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/Candidate/interface/Candidate.h"
0011 
0012 #include "DataFormats/PatCandidates/interface/Jet.h"
0013 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0014 
0015 #include "RecoBTag/FeatureTools/interface/TrackInfoBuilder.h"
0016 #include "RecoBTag/FeatureTools/interface/deep_helpers.h"
0017 #include "RecoBTag/FeatureTools/interface/sorting_modules.h"
0018 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0019 
0020 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0021 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0022 
0023 #include "DataFormats/BTauReco/interface/DeepBoostedJetTagInfo.h"
0024 
0025 #include "Math/PtEtaPhiE4D.h"
0026 #include "Math/LorentzVector.h"
0027 
0028 using namespace btagbtvdeep;
0029 
0030 class GlobalParticleTransformerAK8TagInfoProducer : public edm::stream::EDProducer<> {
0031 public:
0032   explicit GlobalParticleTransformerAK8TagInfoProducer(const edm::ParameterSet &);
0033   ~GlobalParticleTransformerAK8TagInfoProducer() 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 fillJetFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet);
0048   void fillParticleFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet);
0049   void fillSVFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet);
0050 
0051   const double jet_radius_;
0052   const double min_jet_pt_;
0053   const double max_jet_eta_;
0054   const double min_pt_for_track_properties_;
0055   const bool use_puppiP4_;
0056   const bool include_neutrals_;
0057   const bool sort_by_sip2dsig_;
0058   const double min_puppi_wgt_;
0059   const bool flip_ip_sign_;
0060   const double max_sip3dsig_;
0061 
0062   edm::EDGetTokenT<edm::View<reco::Jet>> jet_token_;
0063   edm::EDGetTokenT<VertexCollection> vtx_token_;
0064   edm::EDGetTokenT<SVCollection> sv_token_;
0065   edm::EDGetTokenT<CandidateView> pfcand_token_;
0066   edm::EDGetTokenT<CandidateView> lt_token_;
0067 
0068   bool use_puppi_value_map_;
0069   bool use_pvasq_value_map_;
0070 
0071   edm::EDGetTokenT<edm::ValueMap<float>> puppi_value_map_token_;
0072   edm::EDGetTokenT<edm::ValueMap<int>> pvasq_value_map_token_;
0073   edm::EDGetTokenT<edm::Association<VertexCollection>> pvas_token_;
0074   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> track_builder_token_;
0075 
0076   edm::Handle<VertexCollection> vtxs_;
0077   edm::Handle<SVCollection> svs_;
0078   edm::Handle<CandidateView> pfcands_;
0079   edm::Handle<CandidateView> lts_;
0080   edm::ESHandle<TransientTrackBuilder> track_builder_;
0081   edm::Handle<edm::ValueMap<float>> puppi_value_map_;
0082   edm::Handle<edm::ValueMap<int>> pvasq_value_map_;
0083   edm::Handle<edm::Association<VertexCollection>> pvas_;
0084 
0085   const static std::vector<std::string> jet_features_;
0086   const static std::vector<std::string> charged_particle_features_;
0087   const static std::vector<std::string> neutral_particle_features_;
0088   const static std::vector<std::string> sv_features_;
0089   const reco::Vertex *pv_ = nullptr;
0090 };
0091 
0092 const std::vector<std::string> GlobalParticleTransformerAK8TagInfoProducer::jet_features_{
0093     "jet_pt_log",
0094     "jet_eta",
0095     "jet_mass_log",
0096     "jet_energy_log",
0097 };
0098 const std::vector<std::string> GlobalParticleTransformerAK8TagInfoProducer::charged_particle_features_{
0099     "cpfcandlt_puppiw",
0100     "cpfcandlt_hcalFrac",
0101     "cpfcandlt_VTX_ass",
0102     "cpfcandlt_lostInnerHits",
0103     "cpfcandlt_quality",
0104     "cpfcandlt_charge",
0105     "cpfcandlt_isEl",
0106     "cpfcandlt_isMu",
0107     "cpfcandlt_isChargedHad",
0108     "cpfcandlt_phirel",
0109     "cpfcandlt_etarel",
0110     "cpfcandlt_deltaR",
0111     "cpfcandlt_abseta",
0112     "cpfcandlt_ptrel_log",
0113     "cpfcandlt_erel_log",
0114     "cpfcandlt_pt_log",
0115     "cpfcandlt_drminsv",
0116     "cpfcandlt_drsubjet1",
0117     "cpfcandlt_drsubjet2",
0118     "cpfcandlt_normchi2",
0119     "cpfcandlt_dz",
0120     "cpfcandlt_dzsig",
0121     "cpfcandlt_dxy",
0122     "cpfcandlt_dxysig",
0123     "cpfcandlt_dptdpt",
0124     "cpfcandlt_detadeta",
0125     "cpfcandlt_dphidphi",
0126     "cpfcandlt_dxydxy",
0127     "cpfcandlt_dzdz",
0128     "cpfcandlt_dxydz",
0129     "cpfcandlt_dphidxy",
0130     "cpfcandlt_dlambdadz",
0131     "cpfcandlt_btagEtaRel",
0132     "cpfcandlt_btagPtRatio",
0133     "cpfcandlt_btagPParRatio",
0134     "cpfcandlt_btagSip2dVal",
0135     "cpfcandlt_btagSip2dSig",
0136     "cpfcandlt_btagSip3dVal",
0137     "cpfcandlt_btagSip3dSig",
0138     "cpfcandlt_btagJetDistVal",
0139     "cpfcandlt_mask",
0140     "cpfcandlt_pt_log_nopuppi",
0141     "cpfcandlt_e_log_nopuppi",
0142     "cpfcandlt_ptrel",
0143     "cpfcandlt_erel",
0144     "cpfcandlt_isLostTrack",
0145     "cpfcandlt_pixelBarrelLayersWithMeasurement",
0146     "cpfcandlt_pixelEndcapLayersWithMeasurement",
0147     "cpfcandlt_stripTECLayersWithMeasurement",
0148     "cpfcandlt_stripTIBLayersWithMeasurement",
0149     "cpfcandlt_stripTIDLayersWithMeasurement",
0150     "cpfcandlt_stripTOBLayersWithMeasurement",
0151     "cpfcandlt_px",
0152     "cpfcandlt_py",
0153     "cpfcandlt_pz",
0154     "cpfcandlt_energy"};
0155 
0156 const std::vector<std::string> GlobalParticleTransformerAK8TagInfoProducer::neutral_particle_features_{
0157     "npfcand_puppiw",
0158     "npfcand_hcalFrac",
0159     "npfcand_isGamma",
0160     "npfcand_isNeutralHad",
0161     "npfcand_phirel",
0162     "npfcand_etarel",
0163     "npfcand_deltaR",
0164     "npfcand_abseta",
0165     "npfcand_ptrel_log",
0166     "npfcand_erel_log",
0167     "npfcand_pt_log",
0168     "npfcand_mask",
0169     "npfcand_pt_log_nopuppi",
0170     "npfcand_e_log_nopuppi",
0171     "npfcand_ptrel",
0172     "npfcand_erel",
0173     "npfcand_px",
0174     "npfcand_py",
0175     "npfcand_pz",
0176     "npfcand_energy"};
0177 
0178 const std::vector<std::string> GlobalParticleTransformerAK8TagInfoProducer::sv_features_{
0179     "sv_mask",         "sv_ptrel",    "sv_erel",      "sv_phirel",   "sv_etarel", "sv_deltaR",
0180     "sv_abseta",       "sv_mass",     "sv_ptrel_log", "sv_erel_log", "sv_pt_log", "sv_pt",
0181     "sv_ntracks",      "sv_normchi2", "sv_dxy",       "sv_dxysig",   "sv_d3d",    "sv_d3dsig",
0182     "sv_costhetasvpv", "sv_px",       "sv_py",        "sv_pz",       "sv_energy"};
0183 
0184 GlobalParticleTransformerAK8TagInfoProducer::GlobalParticleTransformerAK8TagInfoProducer(
0185     const edm::ParameterSet &iConfig)
0186     : jet_radius_(iConfig.getParameter<double>("jet_radius")),
0187       min_jet_pt_(iConfig.getParameter<double>("min_jet_pt")),
0188       max_jet_eta_(iConfig.getParameter<double>("max_jet_eta")),
0189       min_pt_for_track_properties_(iConfig.getParameter<double>("min_pt_for_track_properties")),
0190       use_puppiP4_(iConfig.getParameter<bool>("use_puppiP4")),
0191       include_neutrals_(iConfig.getParameter<bool>("include_neutrals")),
0192       sort_by_sip2dsig_(iConfig.getParameter<bool>("sort_by_sip2dsig")),
0193       min_puppi_wgt_(iConfig.getParameter<double>("min_puppi_wgt")),
0194       flip_ip_sign_(iConfig.getParameter<bool>("flip_ip_sign")),
0195       max_sip3dsig_(iConfig.getParameter<double>("sip3dSigMax")),
0196       jet_token_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"))),
0197       vtx_token_(consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0198       sv_token_(consumes<SVCollection>(iConfig.getParameter<edm::InputTag>("secondary_vertices"))),
0199       pfcand_token_(consumes<CandidateView>(iConfig.getParameter<edm::InputTag>("pf_candidates"))),
0200       lt_token_(consumes<CandidateView>(iConfig.getParameter<edm::InputTag>("lost_tracks"))),
0201       use_puppi_value_map_(false),
0202       use_pvasq_value_map_(false),
0203       track_builder_token_(
0204           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))) {
0205   const auto &puppi_value_map_tag = iConfig.getParameter<edm::InputTag>("puppi_value_map");
0206   if (!puppi_value_map_tag.label().empty()) {
0207     puppi_value_map_token_ = consumes<edm::ValueMap<float>>(puppi_value_map_tag);
0208     use_puppi_value_map_ = true;
0209   }
0210 
0211   const auto &pvas_tag = iConfig.getParameter<edm::InputTag>("vertex_associator");
0212   if (!pvas_tag.label().empty()) {
0213     pvasq_value_map_token_ = consumes<edm::ValueMap<int>>(pvas_tag);
0214     pvas_token_ = consumes<edm::Association<VertexCollection>>(pvas_tag);
0215     use_pvasq_value_map_ = true;
0216   }
0217 
0218   produces<DeepBoostedJetTagInfoCollection>();
0219 }
0220 
0221 GlobalParticleTransformerAK8TagInfoProducer::~GlobalParticleTransformerAK8TagInfoProducer() {}
0222 
0223 void GlobalParticleTransformerAK8TagInfoProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0224   // pfGlobalParticleTransformerAK8TagInfos
0225   edm::ParameterSetDescription desc;
0226   desc.add<double>("jet_radius", 0.8);
0227   desc.add<double>("min_jet_pt", 150);
0228   desc.add<double>("max_jet_eta", 99);
0229   desc.add<double>("min_pt_for_track_properties", -1);
0230   desc.add<bool>("use_puppiP4", true);
0231   desc.add<bool>("include_neutrals", true);
0232   desc.add<bool>("sort_by_sip2dsig", false);
0233   desc.add<double>("min_puppi_wgt", 0.01);
0234   desc.add<bool>("flip_ip_sign", false);
0235   desc.add<double>("sip3dSigMax", -1);
0236   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0237   desc.add<edm::InputTag>("secondary_vertices", edm::InputTag("inclusiveCandidateSecondaryVertices"));
0238   desc.add<edm::InputTag>("pf_candidates", edm::InputTag("particleFlow"));
0239   desc.add<edm::InputTag>("lost_tracks", edm::InputTag("lostTracks"));
0240   desc.add<edm::InputTag>("jets", edm::InputTag("ak8PFJetsPuppi"));
0241   desc.add<edm::InputTag>("puppi_value_map", edm::InputTag("puppi"));
0242   desc.add<edm::InputTag>("vertex_associator", edm::InputTag("primaryVertexAssociation", "original"));
0243   descriptions.add("pfGlobalParticleTransformerAK8TagInfos", desc);
0244 }
0245 
0246 void GlobalParticleTransformerAK8TagInfoProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0247   auto output_tag_infos = std::make_unique<DeepBoostedJetTagInfoCollection>();
0248 
0249   auto jets = iEvent.getHandle(jet_token_);
0250 
0251   iEvent.getByToken(vtx_token_, vtxs_);
0252   if (vtxs_->empty()) {
0253     // produce empty TagInfos in case no primary vertex
0254     iEvent.put(std::move(output_tag_infos));
0255     return;  // exit event
0256   }
0257   // primary vertex
0258   pv_ = &vtxs_->at(0);
0259 
0260   iEvent.getByToken(sv_token_, svs_);
0261   iEvent.getByToken(pfcand_token_, pfcands_);
0262   iEvent.getByToken(lt_token_, lts_);
0263 
0264   track_builder_ = iSetup.getHandle(track_builder_token_);
0265 
0266   if (use_puppi_value_map_) {
0267     iEvent.getByToken(puppi_value_map_token_, puppi_value_map_);
0268   }
0269 
0270   if (use_pvasq_value_map_) {
0271     iEvent.getByToken(pvasq_value_map_token_, pvasq_value_map_);
0272     iEvent.getByToken(pvas_token_, pvas_);
0273   }
0274 
0275   for (std::size_t jet_n = 0; jet_n < jets->size(); jet_n++) {
0276     const auto &jet = (*jets)[jet_n];
0277     edm::RefToBase<reco::Jet> jet_ref(jets, jet_n);
0278 
0279     // create jet features
0280     DeepBoostedJetFeatures features;
0281     // declare all the feature variables (init as empty vector)
0282     for (const auto &name : jet_features_) {
0283       features.add(name);
0284     }
0285     for (const auto &name : charged_particle_features_) {
0286       features.add(name);
0287     }
0288     for (const auto &name : neutral_particle_features_) {
0289       features.add(name);
0290     }
0291     for (const auto &name : sv_features_) {
0292       features.add(name);
0293     }
0294 
0295     // fill values only if above pt threshold and has daughters, otherwise left
0296     // empty
0297     bool fill_vars = true;
0298     if (jet.pt() < min_jet_pt_ || std::abs(jet.eta()) > max_jet_eta_)
0299       fill_vars = false;
0300     if (jet.numberOfDaughters() == 0)
0301       fill_vars = false;
0302 
0303     if (fill_vars) {
0304       fillJetFeatures(features, jet);
0305       fillParticleFeatures(features, jet);
0306       fillSVFeatures(features, jet);
0307 
0308       features.check_consistency(charged_particle_features_);
0309       features.check_consistency(neutral_particle_features_);
0310       features.check_consistency(sv_features_);
0311     }
0312 
0313     // this should always be done even if features are not filled
0314     output_tag_infos->emplace_back(features, jet_ref);
0315   }
0316 
0317   iEvent.put(std::move(output_tag_infos));
0318 }
0319 
0320 void GlobalParticleTransformerAK8TagInfoProducer::fillJetFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet) {
0321   // reserve space
0322   for (const auto &name : jet_features_) {
0323     fts.reserve(name, 1);
0324   }
0325   fts.fill("jet_pt_log", std::log(jet.pt()));
0326   fts.fill("jet_eta", jet.eta());
0327   fts.fill("jet_mass_log", std::log(jet.mass()));
0328   fts.fill("jet_energy_log", std::log(jet.energy()));
0329 }
0330 
0331 void GlobalParticleTransformerAK8TagInfoProducer::fillParticleFeatures(DeepBoostedJetFeatures &fts,
0332                                                                        const reco::Jet &jet) {
0333   // require the input to be a pat::Jet
0334   const auto *patJet = dynamic_cast<const pat::Jet *>(&jet);
0335   if (!patJet) {
0336     throw edm::Exception(edm::errors::InvalidReference) << "Input is not a pat::Jet.";
0337   }
0338 
0339   // do nothing if jet does not have constituents
0340   if (jet.numberOfDaughters() == 0)
0341     return;
0342 
0343   // some jet properties
0344   math::XYZVector jet_dir = jet.momentum().Unit();
0345   GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz());
0346   const float etasign = jet.eta() > 0 ? 1 : -1;
0347 
0348   TrackInfoBuilder trkinfo(track_builder_);
0349 
0350   std::map<reco::CandidatePtr, float> puppi_wgt_cache;
0351   auto puppiWgt = [&](const reco::CandidatePtr &cand) {
0352     const auto *pack_cand = dynamic_cast<const pat::PackedCandidate *>(&(*cand));
0353     const auto *reco_cand = dynamic_cast<const reco::PFCandidate *>(&(*cand));
0354     float wgt = 1.;
0355     if (pack_cand) {
0356       wgt = pack_cand->puppiWeight();
0357     } else if (reco_cand) {
0358       if (use_puppi_value_map_) {
0359         wgt = (*puppi_value_map_)[cand];
0360       } else {
0361         throw edm::Exception(edm::errors::InvalidReference) << "Puppi value map is missing";
0362       }
0363     } else {
0364       throw edm::Exception(edm::errors::InvalidReference) << "Cannot convert to either pat::PackedCandidate or "
0365                                                              "reco::PFCandidate";
0366     }
0367     puppi_wgt_cache[cand] = wgt;
0368     return wgt;
0369   };
0370 
0371   std::vector<reco::CandidatePtr> cpfPtrs, npfPtrs;
0372   std::map<reco::CandidatePtr, bool> isLostTrackMap;
0373 
0374   for (const auto &dau : jet.daughterPtrVector()) {
0375     // remove particles w/ extremely low puppi weights
0376     // [Note] use jet daughters here to get the puppiWgt correctly
0377     if ((puppiWgt(dau)) < min_puppi_wgt_)
0378       continue;
0379     // from here: get the original reco/packed candidate not scaled by the puppi weight
0380     auto cand = pfcands_->ptrAt(dau.key());
0381     // charged candidate selection (for Higgs Interaction Net)
0382     if (!include_neutrals_ && (cand->charge() == 0 || cand->pt() < min_pt_for_track_properties_))
0383       continue;
0384     // only when computing the nagative tagger: remove charged candidates with high sip3d
0385     if (flip_ip_sign_ && cand->charge()) {
0386       trkinfo.buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_);
0387       if (trkinfo.getTrackSip3dSig() > max_sip3dsig_)
0388         continue;
0389     }
0390     if (cand->charge() != 0) {
0391       cpfPtrs.push_back(cand);
0392       isLostTrackMap[cand] = false;
0393     } else {
0394       npfPtrs.push_back(cand);
0395     }
0396   }
0397   // lost tracks: fill to cpfcands
0398   for (size_t i = 0; i < lts_->size(); ++i) {
0399     auto cand = lts_->ptrAt(i);
0400     if (reco::deltaR(*cand, jet) < jet_radius_) {
0401       cpfPtrs.push_back(cand);
0402       isLostTrackMap[cand] = true;
0403       puppi_wgt_cache[cand] = 1.;  // set puppi weight to 1 for lost tracks
0404     }
0405   }
0406 
0407   std::vector<btagbtvdeep::SortingClass<reco::CandidatePtr>> c_sorted;
0408   if (sort_by_sip2dsig_) {
0409     // sort charged pf candidates by 2d impact parameter significance
0410     for (const auto &cand : cpfPtrs) {
0411       trkinfo.buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_);
0412       c_sorted.emplace_back(cand,
0413                             trkinfo.getTrackSip2dSig(),
0414                             -btagbtvdeep::mindrsvpfcand(*svs_, &(*cand), jet_radius_),
0415                             cand->pt() / jet.pt());
0416       std::sort(c_sorted.begin(), c_sorted.end(), btagbtvdeep::SortingClass<reco::CandidatePtr>::compareByABCInv);
0417     }
0418     for (unsigned int i = 0; i < c_sorted.size(); i++) {
0419       const auto &c = c_sorted.at(i);
0420       const auto &cand = c.get();
0421       cpfPtrs.at(i) = cand;
0422     }
0423   } else {
0424     if (use_puppiP4_) {
0425       // sort by Puppi-weighted pt
0426       std::sort(
0427           cpfPtrs.begin(), cpfPtrs.end(), [&puppi_wgt_cache](const reco::CandidatePtr &a, const reco::CandidatePtr &b) {
0428             return puppi_wgt_cache.at(a) * a->pt() > puppi_wgt_cache.at(b) * b->pt();
0429           });
0430       std::sort(
0431           npfPtrs.begin(), npfPtrs.end(), [&puppi_wgt_cache](const reco::CandidatePtr &a, const reco::CandidatePtr &b) {
0432             return puppi_wgt_cache.at(a) * a->pt() > puppi_wgt_cache.at(b) * b->pt();
0433           });
0434     } else {
0435       // sort by original pt (not Puppi-weighted)
0436       std::sort(cpfPtrs.begin(), cpfPtrs.end(), [](const auto &a, const auto &b) { return a->pt() > b->pt(); });
0437       std::sort(npfPtrs.begin(), npfPtrs.end(), [](const auto &a, const auto &b) { return a->pt() > b->pt(); });
0438     }
0439   }
0440 
0441   // reserve space
0442   for (const auto &name : charged_particle_features_) {
0443     fts.reserve(name, cpfPtrs.size());
0444   }
0445 
0446   auto useTrackProperties = [&](const reco::PFCandidate *reco_cand) {
0447     const auto *trk = reco_cand->bestTrack();
0448     return trk != nullptr && trk->pt() > min_pt_for_track_properties_;
0449   };
0450 
0451   for (const auto &cand : cpfPtrs) {
0452     const auto *packed_cand = dynamic_cast<const pat::PackedCandidate *>(&(*cand));
0453     const auto *reco_cand = dynamic_cast<const reco::PFCandidate *>(&(*cand));
0454 
0455     if (!include_neutrals_ &&
0456         ((packed_cand && !packed_cand->hasTrackDetails()) || (reco_cand && !useTrackProperties(reco_cand))))
0457       continue;
0458 
0459     const float ip_sign = flip_ip_sign_ ? -1 : 1;
0460 
0461     auto candP4 = use_puppiP4_ ? puppi_wgt_cache.at(cand) * cand->p4() : cand->p4();
0462     if (packed_cand) {
0463       float hcal_fraction = 0.;
0464       if (packed_cand->pdgId() == 1 || packed_cand->pdgId() == 130) {
0465         hcal_fraction = packed_cand->hcalFraction();
0466       } else if (packed_cand->isIsolatedChargedHadron()) {
0467         hcal_fraction = packed_cand->rawHcalFraction();
0468       }
0469 
0470       fts.fill("cpfcandlt_hcalFrac", hcal_fraction);
0471       fts.fill("cpfcandlt_VTX_ass", packed_cand->pvAssociationQuality());
0472       fts.fill("cpfcandlt_lostInnerHits", packed_cand->lostInnerHits());
0473       fts.fill("cpfcandlt_quality", packed_cand->bestTrack() ? packed_cand->bestTrack()->qualityMask() : 0);
0474 
0475       fts.fill("cpfcandlt_charge", packed_cand->charge());
0476       fts.fill("cpfcandlt_isEl", std::abs(packed_cand->pdgId()) == 11);
0477       fts.fill("cpfcandlt_isMu", std::abs(packed_cand->pdgId()) == 13);
0478       fts.fill("cpfcandlt_isChargedHad", std::abs(packed_cand->pdgId()) == 211);
0479       fts.fill("cpfcandlt_isLostTrack", isLostTrackMap[cand]);
0480 
0481       // impact parameters
0482       fts.fill("cpfcandlt_dz", ip_sign * packed_cand->dz());
0483       fts.fill("cpfcandlt_dxy", ip_sign * packed_cand->dxy());
0484       fts.fill("cpfcandlt_dzsig", packed_cand->bestTrack() ? ip_sign * packed_cand->dz() / packed_cand->dzError() : 0);
0485       fts.fill("cpfcandlt_dxysig",
0486                packed_cand->bestTrack() ? ip_sign * packed_cand->dxy() / packed_cand->dxyError() : 0);
0487 
0488     } else if (reco_cand) {
0489       // get vertex association quality
0490       int pv_ass_quality = 0;  // fallback value
0491       float vtx_ass = 0;
0492       if (use_pvasq_value_map_) {
0493         pv_ass_quality = (*pvasq_value_map_)[cand];
0494         const reco::VertexRef &PV_orig = (*pvas_)[cand];
0495         vtx_ass = vtx_ass_from_pfcand(*reco_cand, pv_ass_quality, PV_orig);
0496       } else {
0497         throw edm::Exception(edm::errors::InvalidReference) << "Vertex association missing";
0498       }
0499 
0500       fts.fill("cpfcandlt_hcalFrac", reco_cand->hcalEnergy() / (reco_cand->ecalEnergy() + reco_cand->hcalEnergy()));
0501       fts.fill("cpfcandlt_VTX_ass", vtx_ass);
0502       fts.fill("cpfcandlt_lostInnerHits", useTrackProperties(reco_cand) ? lost_inner_hits_from_pfcand(*reco_cand) : 0);
0503       fts.fill("cpfcandlt_quality", useTrackProperties(reco_cand) ? quality_from_pfcand(*reco_cand) : 0);
0504 
0505       fts.fill("cpfcandlt_charge", reco_cand->charge());
0506       fts.fill("cpfcandlt_isEl", std::abs(reco_cand->pdgId()) == 11);
0507       fts.fill("cpfcandlt_isMu", std::abs(reco_cand->pdgId()) == 13);
0508       fts.fill("cpfcandlt_isChargedHad", std::abs(reco_cand->pdgId()) == 211);
0509       fts.fill("cpfcandlt_isLostTrack", isLostTrackMap[cand]);
0510 
0511       // impact parameters
0512       const auto *trk = reco_cand->bestTrack();
0513       float dz = trk ? ip_sign * trk->dz(pv_->position()) : 0;
0514       float dxy = trk ? ip_sign * trk->dxy(pv_->position()) : 0;
0515       fts.fill("cpfcandlt_dz", dz);
0516       fts.fill("cpfcandlt_dzsig", trk ? dz / trk->dzError() : 0);
0517       fts.fill("cpfcandlt_dxy", dxy);
0518       fts.fill("cpfcandlt_dxysig", trk ? dxy / trk->dxyError() : 0);
0519     }
0520 
0521     // basic kinematics
0522     fts.fill("cpfcandlt_px", candP4.px());
0523     fts.fill("cpfcandlt_py", candP4.py());
0524     fts.fill("cpfcandlt_pz", candP4.pz());
0525     fts.fill("cpfcandlt_energy", candP4.energy());
0526 
0527     fts.fill("cpfcandlt_puppiw", puppi_wgt_cache.at(cand));
0528     fts.fill("cpfcandlt_phirel", reco::deltaPhi(candP4, jet));
0529     fts.fill("cpfcandlt_etarel", etasign * (candP4.eta() - jet.eta()));
0530     fts.fill("cpfcandlt_deltaR", reco::deltaR(candP4, jet));
0531     fts.fill("cpfcandlt_abseta", std::abs(candP4.eta()));
0532 
0533     fts.fill("cpfcandlt_ptrel_log", std::log(candP4.pt() / jet.pt()));
0534     fts.fill("cpfcandlt_ptrel", candP4.pt() / jet.pt());
0535     fts.fill("cpfcandlt_erel_log", std::log(candP4.energy() / jet.energy()));
0536     fts.fill("cpfcandlt_erel", candP4.energy() / jet.energy());
0537     fts.fill("cpfcandlt_pt_log", std::log(candP4.pt()));
0538 
0539     fts.fill("cpfcandlt_mask", 1);
0540     fts.fill("cpfcandlt_pt_log_nopuppi", std::log(cand->pt()));
0541     fts.fill("cpfcandlt_e_log_nopuppi", std::log(cand->energy()));
0542 
0543     float drminpfcandsv = btagbtvdeep::mindrsvpfcand(*svs_, &(*cand), std::numeric_limits<float>::infinity());
0544     fts.fill("cpfcandlt_drminsv", drminpfcandsv);
0545 
0546     // subjets
0547     if (patJet->nSubjetCollections() > 0) {
0548       auto subjets = patJet->subjets();
0549       std::sort(subjets.begin(), subjets.end(), [](const edm::Ptr<pat::Jet> &p1, const edm::Ptr<pat::Jet> &p2) {
0550         return p1->pt() > p2->pt();
0551       });  // sort by pt
0552       fts.fill("cpfcandlt_drsubjet1", !subjets.empty() ? reco::deltaR(*cand, *subjets.at(0)) : -1);
0553       fts.fill("cpfcandlt_drsubjet2", subjets.size() > 1 ? reco::deltaR(*cand, *subjets.at(1)) : -1);
0554     } else {
0555       fts.fill("cpfcandlt_drsubjet1", -1);
0556       fts.fill("cpfcandlt_drsubjet2", -1);
0557     }
0558 
0559     const reco::Track *trk = nullptr;
0560     if (packed_cand) {
0561       trk = packed_cand->bestTrack();
0562     } else if (reco_cand && useTrackProperties(reco_cand)) {
0563       trk = reco_cand->bestTrack();
0564     }
0565     if (trk) {
0566       fts.fill("cpfcandlt_normchi2", std::floor(trk->normalizedChi2()));
0567 
0568       // track covariance
0569       auto cov = [&](unsigned i, unsigned j) { return trk->covariance(i, j); };
0570       fts.fill("cpfcandlt_dptdpt", cov(0, 0));
0571       fts.fill("cpfcandlt_detadeta", cov(1, 1));
0572       fts.fill("cpfcandlt_dphidphi", cov(2, 2));
0573       fts.fill("cpfcandlt_dxydxy", cov(3, 3));
0574       fts.fill("cpfcandlt_dzdz", cov(4, 4));
0575       fts.fill("cpfcandlt_dxydz", cov(3, 4));
0576       fts.fill("cpfcandlt_dphidxy", cov(2, 3));
0577       fts.fill("cpfcandlt_dlambdadz", cov(1, 4));
0578 
0579       trkinfo.buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_);
0580       fts.fill("cpfcandlt_btagEtaRel", trkinfo.getTrackEtaRel());
0581       fts.fill("cpfcandlt_btagPtRatio", trkinfo.getTrackPtRatio());
0582       fts.fill("cpfcandlt_btagPParRatio", trkinfo.getTrackPParRatio());
0583       fts.fill("cpfcandlt_btagSip2dVal", ip_sign * trkinfo.getTrackSip2dVal());
0584       fts.fill("cpfcandlt_btagSip2dSig", ip_sign * trkinfo.getTrackSip2dSig());
0585       fts.fill("cpfcandlt_btagSip3dVal", ip_sign * trkinfo.getTrackSip3dVal());
0586       fts.fill("cpfcandlt_btagSip3dSig", ip_sign * trkinfo.getTrackSip3dSig());
0587       fts.fill("cpfcandlt_btagJetDistVal", trkinfo.getTrackJetDistVal());
0588 
0589       fts.fill("cpfcandlt_pixelBarrelLayersWithMeasurement", trk->hitPattern().pixelBarrelLayersWithMeasurement());
0590       fts.fill("cpfcandlt_pixelEndcapLayersWithMeasurement", trk->hitPattern().pixelEndcapLayersWithMeasurement());
0591       fts.fill("cpfcandlt_stripTIBLayersWithMeasurement", trk->hitPattern().stripTIBLayersWithMeasurement());
0592       fts.fill("cpfcandlt_stripTIDLayersWithMeasurement", trk->hitPattern().stripTIDLayersWithMeasurement());
0593       fts.fill("cpfcandlt_stripTOBLayersWithMeasurement", trk->hitPattern().stripTOBLayersWithMeasurement());
0594       fts.fill("cpfcandlt_stripTECLayersWithMeasurement", trk->hitPattern().stripTECLayersWithMeasurement());
0595 
0596     } else {
0597       fts.fill("cpfcandlt_normchi2", 999);
0598 
0599       fts.fill("cpfcandlt_dptdpt", 0);
0600       fts.fill("cpfcandlt_detadeta", 0);
0601       fts.fill("cpfcandlt_dphidphi", 0);
0602       fts.fill("cpfcandlt_dxydxy", 0);
0603       fts.fill("cpfcandlt_dzdz", 0);
0604       fts.fill("cpfcandlt_dxydz", 0);
0605       fts.fill("cpfcandlt_dphidxy", 0);
0606       fts.fill("cpfcandlt_dlambdadz", 0);
0607 
0608       fts.fill("cpfcandlt_btagEtaRel", 0);
0609       fts.fill("cpfcandlt_btagPtRatio", 0);
0610       fts.fill("cpfcandlt_btagPParRatio", 0);
0611       fts.fill("cpfcandlt_btagSip2dVal", 0);
0612       fts.fill("cpfcandlt_btagSip2dSig", 0);
0613       fts.fill("cpfcandlt_btagSip3dVal", 0);
0614       fts.fill("cpfcandlt_btagSip3dSig", 0);
0615       fts.fill("cpfcandlt_btagJetDistVal", 0);
0616 
0617       fts.fill("cpfcandlt_pixelBarrelLayersWithMeasurement", 0);
0618       fts.fill("cpfcandlt_pixelEndcapLayersWithMeasurement", 0);
0619       fts.fill("cpfcandlt_stripTIBLayersWithMeasurement", 0);
0620       fts.fill("cpfcandlt_stripTIDLayersWithMeasurement", 0);
0621       fts.fill("cpfcandlt_stripTOBLayersWithMeasurement", 0);
0622       fts.fill("cpfcandlt_stripTECLayersWithMeasurement", 0);
0623     }
0624 
0625     // pixel hits pattern variables
0626   }
0627 
0628   // fill neutral candidate features
0629   for (const auto &cand : npfPtrs) {
0630     const auto *packed_cand = dynamic_cast<const pat::PackedCandidate *>(&(*cand));
0631     const auto *reco_cand = dynamic_cast<const reco::PFCandidate *>(&(*cand));
0632 
0633     if (!include_neutrals_ &&
0634         ((packed_cand && !packed_cand->hasTrackDetails()) || (reco_cand && !useTrackProperties(reco_cand))))
0635       continue;
0636 
0637     auto candP4 = use_puppiP4_ ? puppi_wgt_cache.at(cand) * cand->p4() : cand->p4();
0638     if (packed_cand) {
0639       float hcal_fraction = 0.;
0640       if (packed_cand->pdgId() == 1 || packed_cand->pdgId() == 130) {
0641         hcal_fraction = packed_cand->hcalFraction();
0642       } else if (packed_cand->isIsolatedChargedHadron()) {
0643         hcal_fraction = packed_cand->rawHcalFraction();
0644       }
0645 
0646       fts.fill("npfcand_hcalFrac", hcal_fraction);
0647 
0648       fts.fill("npfcand_isGamma", std::abs(packed_cand->pdgId()) == 22);
0649       fts.fill("npfcand_isNeutralHad", std::abs(packed_cand->pdgId()) == 130);
0650 
0651     } else if (reco_cand) {
0652       fts.fill("npfcand_hcalFrac", reco_cand->hcalEnergy() / (reco_cand->ecalEnergy() + reco_cand->hcalEnergy()));
0653 
0654       fts.fill("npfcand_isGamma", std::abs(reco_cand->pdgId()) == 22);
0655       fts.fill("npfcand_isNeutralHad", std::abs(reco_cand->pdgId()) == 130);
0656     }
0657 
0658     // basic kinematics
0659     fts.fill("npfcand_px", candP4.px());
0660     fts.fill("npfcand_py", candP4.py());
0661     fts.fill("npfcand_pz", candP4.pz());
0662     fts.fill("npfcand_energy", candP4.energy());
0663 
0664     fts.fill("npfcand_puppiw", puppi_wgt_cache.at(cand));
0665     fts.fill("npfcand_phirel", reco::deltaPhi(candP4, jet));
0666     fts.fill("npfcand_etarel", etasign * (candP4.eta() - jet.eta()));
0667     fts.fill("npfcand_deltaR", reco::deltaR(candP4, jet));
0668     fts.fill("npfcand_abseta", std::abs(candP4.eta()));
0669 
0670     fts.fill("npfcand_ptrel_log", std::log(candP4.pt() / jet.pt()));
0671     fts.fill("npfcand_ptrel", candP4.pt() / jet.pt());
0672     fts.fill("npfcand_erel_log", std::log(candP4.energy() / jet.energy()));
0673     fts.fill("npfcand_erel", candP4.energy() / jet.energy());
0674     fts.fill("npfcand_pt_log", std::log(candP4.pt()));
0675 
0676     fts.fill("npfcand_mask", 1);
0677     fts.fill("npfcand_pt_log_nopuppi", std::log(cand->pt()));
0678     fts.fill("npfcand_e_log_nopuppi", std::log(cand->energy()));
0679   }
0680 }
0681 
0682 void GlobalParticleTransformerAK8TagInfoProducer::fillSVFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet) {
0683   std::vector<const reco::VertexCompositePtrCandidate *> jetSVs;
0684   for (const auto &sv : *svs_) {
0685     if (reco::deltaR2(sv, jet) < jet_radius_ * jet_radius_) {
0686       jetSVs.push_back(&sv);
0687     }
0688   }
0689   // sort by dxy significance
0690   std::sort(jetSVs.begin(),
0691             jetSVs.end(),
0692             [&](const reco::VertexCompositePtrCandidate *sva, const reco::VertexCompositePtrCandidate *svb) {
0693               return sv_vertex_comparator(*sva, *svb, *pv_);
0694             });
0695 
0696   // reserve space
0697   for (const auto &name : sv_features_) {
0698     fts.reserve(name, jetSVs.size());
0699   }
0700 
0701   const float etasign = jet.eta() > 0 ? 1 : -1;
0702 
0703   for (const auto *sv : jetSVs) {
0704     // basic kinematics
0705     fts.fill("sv_mask", 1);
0706 
0707     fts.fill("sv_px", sv->px());
0708     fts.fill("sv_py", sv->py());
0709     fts.fill("sv_pz", sv->pz());
0710     fts.fill("sv_energy", sv->energy());
0711 
0712     fts.fill("sv_phirel", reco::deltaPhi(*sv, jet));
0713     fts.fill("sv_etarel", etasign * (sv->eta() - jet.eta()));
0714     fts.fill("sv_deltaR", reco::deltaR(*sv, jet));
0715     fts.fill("sv_abseta", std::abs(sv->eta()));
0716     fts.fill("sv_mass", sv->mass());
0717 
0718     fts.fill("sv_ptrel_log", std::log(sv->pt() / jet.pt()));
0719     fts.fill("sv_ptrel", sv->pt() / jet.pt());
0720     fts.fill("sv_erel_log", std::log(sv->energy() / jet.energy()));
0721     fts.fill("sv_erel", sv->energy() / jet.energy());
0722     fts.fill("sv_pt_log", std::log(sv->pt()));
0723     fts.fill("sv_pt", sv->pt());
0724 
0725     // sv properties
0726     fts.fill("sv_ntracks", sv->numberOfDaughters());
0727     fts.fill("sv_normchi2", sv->vertexNormalizedChi2());
0728 
0729     const auto &dxy = vertexDxy(*sv, *pv_);
0730     fts.fill("sv_dxy", dxy.value());
0731     fts.fill("sv_dxysig", dxy.significance());
0732 
0733     const auto &d3d = vertexD3d(*sv, *pv_);
0734     fts.fill("sv_d3d", d3d.value());
0735     fts.fill("sv_d3dsig", d3d.significance());
0736 
0737     fts.fill("sv_costhetasvpv", (flip_ip_sign_ ? -1.f : 1.f) * vertexDdotP(*sv, *pv_));
0738   }
0739 }
0740 
0741 // define this as a plug-in
0742 DEFINE_FWK_MODULE(GlobalParticleTransformerAK8TagInfoProducer);
[ 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!