Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-09 02:22:32

0001 #include <TVector3.h>
0002 #include <TTree.h>
0003 
0004 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0005 
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/StreamID.h"
0013 #include "FWCore/Utilities/interface/ESGetToken.h"
0014 
0015 #include "DataFormats/Candidate/interface/Candidate.h"
0016 #include "DataFormats/PatCandidates/interface/Jet.h"
0017 #include "DataFormats/PatCandidates/interface/Muon.h"
0018 #include "DataFormats/PatCandidates/interface/Electron.h"
0019 #include "DataFormats/PatCandidates/interface/Photon.h"
0020 #include "DataFormats/PatCandidates/interface/Tau.h"
0021 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0022 
0023 #include "RecoBTag/FeatureTools/interface/TrackInfoBuilder.h"
0024 #include "RecoBTag/FeatureTools/interface/deep_helpers.h"
0025 #include "RecoBTag/FeatureTools/interface/sorting_modules.h"
0026 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0027 
0028 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0029 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0030 #include "DataFormats/BTauReco/interface/DeepBoostedJetTagInfo.h"
0031 
0032 using namespace btagbtvdeep;
0033 
0034 class ParticleNetFeatureEvaluator : public edm::stream::EDProducer<> {
0035 public:
0036   explicit ParticleNetFeatureEvaluator(const edm::ParameterSet &);
0037   ~ParticleNetFeatureEvaluator() override;
0038 
0039   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0040 
0041 private:
0042   void beginStream(edm::StreamID) override {}
0043   void produce(edm::Event &, const edm::EventSetup &) override;
0044   void endStream() override {}
0045   void fillParticleFeatures(DeepBoostedJetFeatures &fts,
0046                             const reco::Jet &jet,
0047                             const std::vector<math::XYZTLorentzVector> &tau_pfcandidates,
0048                             const pat::MuonCollection &muons,
0049                             const pat::ElectronCollection &electrons,
0050                             const pat::PhotonCollection &photons);
0051   void fillSVFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet);
0052   void fillLostTrackFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet);
0053   bool useTrackProperties(const pat::PackedCandidate *cand);
0054 
0055   const double jet_radius_;
0056   const double min_jet_pt_;
0057   const double max_jet_eta_;
0058   const double min_jet_eta_;
0059   const double min_pt_for_track_properties_;
0060   const double min_pt_for_pfcandidates_;
0061   const double min_pt_for_losttrack_;
0062   const double max_dr_for_losttrack_;
0063   const double min_pt_for_taus_;
0064   const double max_eta_for_taus_;
0065   const bool include_neutrals_;
0066   const bool flip_ip_sign_;
0067   const double max_sip3dsig_for_flip_;
0068 
0069   edm::EDGetTokenT<pat::MuonCollection> muon_token_;
0070   edm::EDGetTokenT<pat::ElectronCollection> electron_token_;
0071   edm::EDGetTokenT<pat::PhotonCollection> photon_token_;
0072   edm::EDGetTokenT<pat::TauCollection> tau_token_;
0073   edm::EDGetTokenT<edm::View<reco::Jet>> jet_token_;
0074   edm::EDGetTokenT<pat::PackedCandidateCollection> losttrack_token_;
0075   edm::EDGetTokenT<reco::VertexCollection> vtx_token_;
0076   edm::EDGetTokenT<reco::VertexCompositePtrCandidateCollection> sv_token_;
0077   edm::EDGetTokenT<edm::View<reco::Candidate>> pfcand_token_;
0078   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> track_builder_token_;
0079 
0080   edm::Handle<reco::VertexCollection> vtxs_;
0081   edm::Handle<reco::VertexCompositePtrCandidateCollection> svs_;
0082   edm::Handle<edm::View<reco::Candidate>> pfcands_;
0083   edm::Handle<pat::PackedCandidateCollection> losttracks_;
0084   edm::ESHandle<TransientTrackBuilder> track_builder_;
0085 
0086   const static std::vector<std::string> particle_features_;
0087   const static std::vector<std::string> sv_features_;
0088   const static std::vector<std::string> losttrack_features_;
0089   const reco::Vertex *pv_ = nullptr;
0090 
0091   TTree *tree;
0092   unsigned int event;
0093   float jet_pt;
0094   float jet_pt_raw;
0095   float jet_eta;
0096   float jet_phi;
0097   float jet_mass;
0098   unsigned int ijet;
0099   std::vector<float> jet_pfcand_pt_log;
0100   std::vector<float> jet_pfcand_energy_log;
0101   std::vector<float> jet_pfcand_deta;
0102   std::vector<float> jet_pfcand_dphi;
0103   std::vector<float> jet_pfcand_eta;
0104   std::vector<float> jet_pfcand_charge;
0105   std::vector<float> jet_pfcand_frompv;
0106   std::vector<float> jet_pfcand_nlostinnerhits;
0107   std::vector<float> jet_pfcand_track_chi2;
0108   std::vector<float> jet_pfcand_track_qual;
0109   std::vector<float> jet_pfcand_dz;
0110   std::vector<float> jet_pfcand_dzsig;
0111   std::vector<float> jet_pfcand_dxy;
0112   std::vector<float> jet_pfcand_dxysig;
0113   std::vector<float> jet_pfcand_etarel;
0114   std::vector<float> jet_pfcand_pperp_ratio;
0115   std::vector<float> jet_pfcand_ppara_ratio;
0116   std::vector<float> jet_pfcand_trackjet_d3d;
0117   std::vector<float> jet_pfcand_trackjet_d3dsig;
0118   std::vector<float> jet_pfcand_trackjet_dist;
0119   std::vector<float> jet_pfcand_npixhits;
0120   std::vector<float> jet_pfcand_nstriphits;
0121   std::vector<float> jet_pfcand_trackjet_decayL;
0122   std::vector<float> jet_pfcand_id;
0123   std::vector<float> jet_pfcand_calofraction;
0124   std::vector<float> jet_pfcand_hcalfraction;
0125   std::vector<float> jet_pfcand_puppiw;
0126   std::vector<float> jet_pfcand_muon_id;
0127   std::vector<float> jet_pfcand_muon_isglobal;
0128   std::vector<float> jet_pfcand_muon_chi2;
0129   std::vector<float> jet_pfcand_muon_segcomp;
0130   std::vector<float> jet_pfcand_muon_nvalidhit;
0131   std::vector<float> jet_pfcand_muon_nstation;
0132   std::vector<float> jet_pfcand_electron_detaIn;
0133   std::vector<float> jet_pfcand_electron_dphiIn;
0134   std::vector<float> jet_pfcand_electron_sigIetaIeta;
0135   std::vector<float> jet_pfcand_electron_sigIphiIphi;
0136   std::vector<float> jet_pfcand_electron_r9;
0137   std::vector<float> jet_pfcand_electron_convProb;
0138   std::vector<float> jet_pfcand_photon_sigIetaIeta;
0139   std::vector<float> jet_pfcand_photon_r9;
0140   std::vector<float> jet_pfcand_photon_eVeto;
0141   std::vector<float> jet_pfcand_tau_signal;
0142   std::vector<float> jet_sv_pt_log;
0143   std::vector<float> jet_sv_mass;
0144   std::vector<float> jet_sv_deta;
0145   std::vector<float> jet_sv_dphi;
0146   std::vector<float> jet_sv_eta;
0147   std::vector<float> jet_sv_ntrack;
0148   std::vector<float> jet_sv_chi2;
0149   std::vector<float> jet_sv_dxy;
0150   std::vector<float> jet_sv_dxysig;
0151   std::vector<float> jet_sv_d3d;
0152   std::vector<float> jet_sv_d3dsig;
0153   std::vector<float> jet_losttrack_pt_log;
0154   std::vector<float> jet_losttrack_eta;
0155   std::vector<float> jet_losttrack_deta;
0156   std::vector<float> jet_losttrack_dphi;
0157   std::vector<float> jet_losttrack_charge;
0158   std::vector<float> jet_losttrack_frompv;
0159   std::vector<float> jet_losttrack_track_chi2;
0160   std::vector<float> jet_losttrack_track_qual;
0161   std::vector<float> jet_losttrack_dz;
0162   std::vector<float> jet_losttrack_dxy;
0163   std::vector<float> jet_losttrack_dzsig;
0164   std::vector<float> jet_losttrack_dxysig;
0165   std::vector<float> jet_losttrack_etarel;
0166   std::vector<float> jet_losttrack_trackjet_d3d;
0167   std::vector<float> jet_losttrack_trackjet_d3dsig;
0168   std::vector<float> jet_losttrack_trackjet_dist;
0169   std::vector<float> jet_losttrack_trackjet_decayL;
0170   std::vector<float> jet_losttrack_npixhits;
0171   std::vector<float> jet_losttrack_nstriphits;
0172 };
0173 
0174 const std::vector<std::string> ParticleNetFeatureEvaluator::particle_features_{"jet_pfcand_pt_log",
0175                                                                                "jet_pfcand_energy_log",
0176                                                                                "jet_pfcand_deta",
0177                                                                                "jet_pfcand_dphi",
0178                                                                                "jet_pfcand_eta",
0179                                                                                "jet_pfcand_charge",
0180                                                                                "jet_pfcand_frompv",
0181                                                                                "jet_pfcand_nlostinnerhits",
0182                                                                                "jet_pfcand_track_chi2",
0183                                                                                "jet_pfcand_track_qual",
0184                                                                                "jet_pfcand_dz",
0185                                                                                "jet_pfcand_dzsig",
0186                                                                                "jet_pfcand_dxy",
0187                                                                                "jet_pfcand_dxysig",
0188                                                                                "jet_pfcand_etarel",
0189                                                                                "jet_pfcand_pperp_ratio",
0190                                                                                "jet_pfcand_ppara_ratio",
0191                                                                                "jet_pfcand_trackjet_d3d",
0192                                                                                "jet_pfcand_trackjet_d3dsig",
0193                                                                                "jet_pfcand_trackjet_dist",
0194                                                                                "jet_pfcand_nhits",
0195                                                                                "jet_pfcand_npixhits",
0196                                                                                "jet_pfcand_nstriphits",
0197                                                                                "jet_pfcand_trackjet_decayL",
0198                                                                                "jet_pfcand_id",
0199                                                                                "jet_pfcand_calofraction",
0200                                                                                "jet_pfcand_hcalfraction",
0201                                                                                "jet_pfcand_puppiw",
0202                                                                                "jet_pfcand_muon_id",
0203                                                                                "jet_pfcand_muon_isglobal",
0204                                                                                "jet_pfcand_muon_segcomp",
0205                                                                                "jet_pfcand_muon_chi2",
0206                                                                                "jet_pfcand_muon_nvalidhit",
0207                                                                                "jet_pfcand_muon_nstation",
0208                                                                                "jet_pfcand_electron_detaIn",
0209                                                                                "jet_pfcand_electron_dphiIn",
0210                                                                                "jet_pfcand_electron_sigIetaIeta",
0211                                                                                "jet_pfcand_electron_sigIphiIphi",
0212                                                                                "jet_pfcand_electron_r9",
0213                                                                                "jet_pfcand_electron_convProb",
0214                                                                                "jet_pfcand_photon_sigIetaIeta",
0215                                                                                "jet_pfcand_photon_r9",
0216                                                                                "jet_pfcand_photon_eVeto",
0217                                                                                "jet_pfcand_tau_signal",
0218                                                                                "pfcand_mask"};
0219 
0220 const std::vector<std::string> ParticleNetFeatureEvaluator::sv_features_{"jet_sv_pt_log",
0221                                                                          "jet_sv_mass",
0222                                                                          "jet_sv_deta",
0223                                                                          "jet_sv_dphi",
0224                                                                          "jet_sv_eta",
0225                                                                          "jet_sv_ntrack",
0226                                                                          "jet_sv_chi2",
0227                                                                          "jet_sv_dxy",
0228                                                                          "jet_sv_dxysig",
0229                                                                          "jet_sv_d3d",
0230                                                                          "jet_sv_d3dsig",
0231                                                                          "sv_mask"};
0232 
0233 const std::vector<std::string> ParticleNetFeatureEvaluator::losttrack_features_{"jet_losttrack_pt_log",
0234                                                                                 "jet_losttrack_eta",
0235                                                                                 "jet_losttrack_deta",
0236                                                                                 "jet_losttrack_dphi",
0237                                                                                 "jet_losttrack_charge",
0238                                                                                 "jet_losttrack_frompv",
0239                                                                                 "jet_losttrack_track_chi2",
0240                                                                                 "jet_losttrack_track_qual",
0241                                                                                 "jet_losttrack_dz",
0242                                                                                 "jet_losttrack_dxy",
0243                                                                                 "jet_losttrack_dzsig",
0244                                                                                 "jet_losttrack_dxysig",
0245                                                                                 "jet_losttrack_etarel",
0246                                                                                 "jet_losttrack_trackjet_d3d",
0247                                                                                 "jet_losttrack_trackjet_d3dsig",
0248                                                                                 "jet_losttrack_trackjet_dist",
0249                                                                                 "jet_losttrack_trackjet_decayL",
0250                                                                                 "jet_losttrack_npixhits",
0251                                                                                 "jet_losttrack_nstriphits",
0252                                                                                 "lt_mask"};
0253 
0254 ParticleNetFeatureEvaluator::ParticleNetFeatureEvaluator(const edm::ParameterSet &iConfig)
0255     : jet_radius_(iConfig.getParameter<double>("jet_radius")),
0256       min_jet_pt_(iConfig.getParameter<double>("min_jet_pt")),
0257       max_jet_eta_(iConfig.getParameter<double>("max_jet_eta")),
0258       min_jet_eta_(iConfig.getParameter<double>("min_jet_eta")),
0259       min_pt_for_track_properties_(iConfig.getParameter<double>("min_pt_for_track_properties")),
0260       min_pt_for_pfcandidates_(iConfig.getParameter<double>("min_pt_for_pfcandidates")),
0261       min_pt_for_losttrack_(iConfig.getParameter<double>("min_pt_for_losttrack")),
0262       max_dr_for_losttrack_(iConfig.getParameter<double>("max_dr_for_losttrack")),
0263       min_pt_for_taus_(iConfig.getParameter<double>("min_pt_for_taus")),
0264       max_eta_for_taus_(iConfig.getParameter<double>("max_eta_for_taus")),
0265       include_neutrals_(iConfig.getParameter<bool>("include_neutrals")),
0266       flip_ip_sign_(iConfig.getParameter<bool>("flip_ip_sign")),
0267       max_sip3dsig_for_flip_(iConfig.getParameter<double>("max_sip3dsig_for_flip")),
0268       muon_token_(consumes<pat::MuonCollection>(iConfig.getParameter<edm::InputTag>("muons"))),
0269       electron_token_(consumes<pat::ElectronCollection>(iConfig.getParameter<edm::InputTag>("electrons"))),
0270       photon_token_(consumes<pat::PhotonCollection>(iConfig.getParameter<edm::InputTag>("photons"))),
0271       tau_token_(consumes<pat::TauCollection>(iConfig.getParameter<edm::InputTag>("taus"))),
0272       jet_token_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"))),
0273       losttrack_token_(consumes<pat::PackedCandidateCollection>(iConfig.getParameter<edm::InputTag>("losttracks"))),
0274       vtx_token_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0275       sv_token_(consumes<reco::VertexCompositePtrCandidateCollection>(
0276           iConfig.getParameter<edm::InputTag>("secondary_vertices"))),
0277       pfcand_token_(consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("pf_candidates"))),
0278       track_builder_token_(
0279           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))) {
0280   produces<std::vector<reco::DeepBoostedJetTagInfo>>();
0281 }
0282 
0283 ParticleNetFeatureEvaluator::~ParticleNetFeatureEvaluator() {}
0284 
0285 void ParticleNetFeatureEvaluator::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0286   // pfDeepBoostedJetTagInfos
0287   edm::ParameterSetDescription desc;
0288   desc.add<double>("jet_radius", 0.8);
0289   desc.add<double>("min_jet_pt", 150);
0290   desc.add<double>("max_jet_eta", 99);
0291   desc.add<double>("min_jet_eta", 0.0);
0292   desc.add<double>("min_pt_for_track_properties", -1);
0293   desc.add<double>("min_pt_for_pfcandidates", -1);
0294   desc.add<double>("min_pt_for_losttrack", 1);
0295   desc.add<double>("max_dr_for_losttrack", 0.4);
0296   desc.add<double>("min_pt_for_taus", 20.);
0297   desc.add<double>("max_eta_for_taus", 2.5);
0298   desc.add<bool>("include_neutrals", true);
0299   desc.add<bool>("flip_ip_sign", false);
0300   desc.add<double>("max_sip3dsig_for_flip", 99999);
0301   desc.add<edm::InputTag>("vertices", edm::InputTag("offlineSlimmedPrimaryVertices"));
0302   desc.add<edm::InputTag>("secondary_vertices", edm::InputTag("slimmedSecondaryVertices"));
0303   desc.add<edm::InputTag>("pf_candidates", edm::InputTag("packedPFCandidates"));
0304   desc.add<edm::InputTag>("losttracks", edm::InputTag("lostTracks"));
0305   desc.add<edm::InputTag>("jets", edm::InputTag("slimmedJetsAK8"));
0306   desc.add<edm::InputTag>("muons", edm::InputTag("slimmedMuons"));
0307   desc.add<edm::InputTag>("taus", edm::InputTag("slimmedTaus"));
0308   desc.add<edm::InputTag>("electrons", edm::InputTag("slimmedElectrons"));
0309   desc.add<edm::InputTag>("photons", edm::InputTag("slimmedPhotons"));
0310   descriptions.add("ParticleNetFeatureEvaluator", desc);
0311 }
0312 
0313 void ParticleNetFeatureEvaluator::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0314   // output collection
0315   auto output_tag_infos = std::make_unique<std::vector<reco::DeepBoostedJetTagInfo>>();
0316   // Input jets
0317   auto jets = iEvent.getHandle(jet_token_);
0318   // Input muons
0319   auto muons = iEvent.getHandle(muon_token_);
0320   // Input taus
0321   auto taus = iEvent.getHandle(tau_token_);
0322   // Input electrons
0323   auto electrons = iEvent.getHandle(electron_token_);
0324   // Input photons
0325   auto photons = iEvent.getHandle(photon_token_);
0326   // Input lost tracks
0327   iEvent.getByToken(losttrack_token_, losttracks_);
0328   // Primary vertexes
0329   iEvent.getByToken(vtx_token_, vtxs_);
0330   if (vtxs_->empty()) {
0331     // produce empty TagInfos in case no primary vertex
0332     iEvent.put(std::move(output_tag_infos));
0333     return;  // exit event
0334   }
0335   // Leading vertex
0336   pv_ = &vtxs_->at(0);
0337   // Secondary vertexs
0338   iEvent.getByToken(sv_token_, svs_);
0339   // PF candidates
0340   iEvent.getByToken(pfcand_token_, pfcands_);
0341   // Track builder
0342   track_builder_ = iSetup.getHandle(track_builder_token_);
0343 
0344   // tau signal candidates
0345   std::vector<math::XYZTLorentzVector> tau_pfcandidates;
0346   for (size_t itau = 0; itau < taus->size(); itau++) {
0347     if (taus->at(itau).pt() < min_pt_for_taus_)
0348       continue;
0349     if (fabs(taus->at(itau).eta()) > max_eta_for_taus_)
0350       continue;
0351     for (unsigned ipart = 0; ipart < taus->at(itau).signalCands().size(); ipart++) {
0352       const pat::PackedCandidate *pfcand =
0353           dynamic_cast<const pat::PackedCandidate *>(taus->at(itau).signalCands()[ipart].get());
0354       tau_pfcandidates.push_back(pfcand->p4());
0355     }
0356   }
0357 
0358   // Loop over jet
0359   for (std::size_t jet_n = 0; jet_n < jets->size(); jet_n++) {
0360     const auto &jet = (*jets)[jet_n];
0361     edm::RefToBase<reco::Jet> jet_ref(jets, jet_n);
0362 
0363     // create jet features
0364     DeepBoostedJetFeatures features;
0365     for (const auto &name : particle_features_)
0366       features.add(name);
0367     for (const auto &name : sv_features_)
0368       features.add(name);
0369 
0370     // fill values only if above pt threshold and has daughters, otherwise left
0371     bool fill_vars = true;
0372     if ((jet.pt() < min_jet_pt_ and
0373          dynamic_cast<const pat::Jet *>(&jet)->correctedJet("Uncorrected").pt() < min_jet_pt_) or
0374         std::abs(jet.eta()) >= max_jet_eta_ or std::abs(jet.eta()) < min_jet_eta_)
0375       fill_vars = false;
0376     if (jet.numberOfDaughters() == 0)
0377       fill_vars = false;
0378 
0379     // fill features
0380     if (fill_vars) {
0381       fillParticleFeatures(features, jet, tau_pfcandidates, *muons, *electrons, *photons);
0382       fillSVFeatures(features, jet);
0383       fillLostTrackFeatures(features, jet);
0384       features.check_consistency(particle_features_);
0385       features.check_consistency(sv_features_);
0386       features.check_consistency(losttrack_features_);
0387     }
0388 
0389     // this should always be done even if features are not filled
0390     output_tag_infos->emplace_back(features, jet_ref);
0391   }
0392   // move output collection
0393   iEvent.put(std::move(output_tag_infos));
0394 }
0395 
0396 bool ParticleNetFeatureEvaluator::useTrackProperties(const pat::PackedCandidate *cand) {
0397   const auto *track = cand->bestTrack();
0398   return track != nullptr and track->pt() > min_pt_for_track_properties_;
0399 };
0400 
0401 void ParticleNetFeatureEvaluator::fillParticleFeatures(DeepBoostedJetFeatures &fts,
0402                                                        const reco::Jet &jet,
0403                                                        const std::vector<math::XYZTLorentzVector> &tau_pfcandidates,
0404                                                        const pat::MuonCollection &muons,
0405                                                        const pat::ElectronCollection &electrons,
0406                                                        const pat::PhotonCollection &photons) {
0407   // some jet properties
0408   math::XYZVector jet_dir = jet.momentum().Unit();
0409   TVector3 jet_direction(jet.momentum().Unit().x(), jet.momentum().Unit().y(), jet.momentum().Unit().z());
0410   GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz());
0411   // vertexes
0412   reco::VertexRefProd PVRefProd(vtxs_);
0413   // track builder
0414   TrackInfoBuilder trackinfo(track_builder_);
0415 
0416   // make list of pf-candidates to be considered
0417   std::vector<const pat::PackedCandidate *> daughters;
0418   for (const auto &dau : jet.daughterPtrVector()) {
0419     // remove particles w/ extremely low puppi weights
0420     const pat::PackedCandidate *cand = dynamic_cast<const pat::PackedCandidate *>(&(*dau));
0421     if (not cand)
0422       throw edm::Exception(edm::errors::InvalidReference) << "Cannot convert to either pat::PackedCandidate";
0423     // base requirements on PF candidates
0424     if (cand->pt() < min_pt_for_pfcandidates_)
0425       continue;
0426     // charged candidate selection (for Higgs Interaction Net)
0427     if (!include_neutrals_ and (cand->charge() == 0 or cand->pt() < min_pt_for_track_properties_))
0428       continue;
0429     // filling daughters
0430     daughters.push_back(cand);
0431   }
0432 
0433   // sort by original pt (not Puppi-weighted)
0434   std::sort(daughters.begin(), daughters.end(), [](const auto &a, const auto &b) { return a->pt() > b->pt(); });
0435 
0436   // reserve space
0437   for (const auto &name : particle_features_)
0438     fts.reserve(name, daughters.size());
0439 
0440   // Build observables
0441   for (const auto &cand : daughters) {
0442     if (!include_neutrals_ and !useTrackProperties(cand))
0443       continue;
0444 
0445     // input particle is a packed PF candidate
0446     auto candP4 = cand->p4();
0447     auto candP3 = cand->momentum();
0448 
0449     // candidate track
0450     const reco::Track *track = nullptr;
0451     if (useTrackProperties(cand))
0452       track = cand->bestTrack();
0453 
0454     // reco-vertex association
0455     reco::VertexRef pv_ass = reco::VertexRef(vtxs_, 0);
0456     math::XYZPoint pv_ass_pos = pv_ass->position();
0457 
0458     TVector3 cand_direction(candP3.x(), candP3.y(), candP3.z());
0459 
0460     // only when computing the nagative tagger: remove charged candidates with high sip3d
0461     if (flip_ip_sign_ && track) {
0462       reco::TransientTrack transientTrack = track_builder_->build(*track);
0463       Measurement1D meas_ip3d = IPTools::signedImpactParameter3D(transientTrack, jet_ref_track_dir, *pv_).second;
0464       if (!std::isnan(meas_ip3d.significance()) && meas_ip3d.significance() > max_sip3dsig_for_flip_) {
0465         continue;
0466       }
0467     }
0468     float ip_sign = flip_ip_sign_ ? -1.f : 1.f;
0469 
0470     fts.fill("jet_pfcand_pt_log", std::isnan(std::log(candP4.pt())) ? 0 : std::log(candP4.pt()));
0471     fts.fill("jet_pfcand_energy_log", std::isnan(std::log(candP4.energy())) ? 0 : std::log(candP4.energy()));
0472     fts.fill("jet_pfcand_eta", candP4.eta());
0473     fts.fill("jet_pfcand_deta", jet_direction.Eta() - cand_direction.Eta());
0474     fts.fill("jet_pfcand_dphi", jet_direction.DeltaPhi(cand_direction));
0475     fts.fill("jet_pfcand_charge", cand->charge());
0476     fts.fill("jet_pfcand_etarel",
0477              std::isnan(reco::btau::etaRel(jet_dir, candP3)) ? 0 : reco::btau::etaRel(jet_dir, candP3));
0478     fts.fill("jet_pfcand_pperp_ratio",
0479              std::isnan(jet_direction.Perp(cand_direction) / cand_direction.Mag())
0480                  ? 0
0481                  : jet_direction.Perp(cand_direction) / cand_direction.Mag());
0482     fts.fill("jet_pfcand_ppara_ratio",
0483              std::isnan(jet_direction.Dot(cand_direction) / cand_direction.Mag())
0484                  ? 0
0485                  : jet_direction.Dot(cand_direction) / cand_direction.Mag());
0486     fts.fill("jet_pfcand_frompv", cand->fromPV());
0487     fts.fill("jet_pfcand_dz", ip_sign * (std::isnan(cand->dz(pv_ass_pos)) ? 0 : cand->dz(pv_ass_pos)));
0488     fts.fill("jet_pfcand_dxy", ip_sign * (std::isnan(cand->dxy(pv_ass_pos)) ? 0 : cand->dxy(pv_ass_pos)));
0489     fts.fill("jet_pfcand_puppiw", cand->puppiWeight());
0490     fts.fill("jet_pfcand_nlostinnerhits", cand->lostInnerHits());
0491     fts.fill("jet_pfcand_nhits", cand->numberOfHits());
0492     fts.fill("jet_pfcand_npixhits", cand->numberOfPixelHits());
0493     fts.fill("jet_pfcand_nstriphits", cand->stripLayersWithMeasurement());
0494 
0495     if (abs(cand->pdgId()) == 11 and cand->charge() != 0)
0496       fts.fill("jet_pfcand_id", 0);
0497     else if (abs(cand->pdgId()) == 13 and cand->charge() != 0)
0498       fts.fill("jet_pfcand_id", 1);
0499     else if (abs(cand->pdgId()) == 22 and cand->charge() == 0)
0500       fts.fill("jet_pfcand_id", 2);
0501     else if (abs(cand->pdgId()) != 22 and cand->charge() == 0 and abs(cand->pdgId()) != 1 and abs(cand->pdgId()) != 2)
0502       fts.fill("jet_pfcand_id", 3);
0503     else if (abs(cand->pdgId()) != 11 and abs(cand->pdgId()) != 13 and cand->charge() != 0)
0504       fts.fill("jet_pfcand_id", 4);
0505     else if (cand->charge() == 0 and abs(cand->pdgId()) == 1)
0506       fts.fill("jet_pfcand_id", 5);
0507     else if (cand->charge() == 0 and abs(cand->pdgId()) == 2)
0508       fts.fill("jet_pfcand_id", 6);
0509     else
0510       fts.fill("jet_pfcand_id", -1);
0511 
0512     fts.fill("jet_pfcand_hcalfraction", std::isnan(cand->hcalFraction()) ? 0 : cand->hcalFraction());
0513     fts.fill("jet_pfcand_calofraction", std::isnan(cand->caloFraction()) ? 0 : cand->caloFraction());
0514     fts.fill("pfcand_mask", 1);
0515 
0516     if (track) {
0517       fts.fill(
0518           "jet_pfcand_dzsig",
0519           std::isnan(fabs(cand->dz(pv_ass_pos)) / cand->dzError()) ? 0 : fabs(cand->dz(pv_ass_pos)) / cand->dzError());
0520       fts.fill("jet_pfcand_dxysig",
0521                std::isnan(fabs(cand->dxy(pv_ass_pos)) / cand->dxyError())
0522                    ? 0
0523                    : fabs(cand->dxy(pv_ass_pos)) / cand->dxyError());
0524       fts.fill("jet_pfcand_track_chi2", track->normalizedChi2());
0525       fts.fill("jet_pfcand_track_qual", track->qualityMask());
0526 
0527       reco::TransientTrack transientTrack = track_builder_->build(*track);
0528       Measurement1D meas_ip2d =
0529           IPTools::signedTransverseImpactParameter(transientTrack, jet_ref_track_dir, *pv_).second;
0530       Measurement1D meas_ip3d = IPTools::signedImpactParameter3D(transientTrack, jet_ref_track_dir, *pv_).second;
0531       Measurement1D meas_jetdist = IPTools::jetTrackDistance(transientTrack, jet_ref_track_dir, *pv_).second;
0532       Measurement1D meas_decayl = IPTools::signedDecayLength3D(transientTrack, jet_ref_track_dir, *pv_).second;
0533 
0534       fts.fill("jet_pfcand_trackjet_d3d", ip_sign * (std::isnan(meas_ip3d.value()) ? 0 : meas_ip3d.value()));
0535       fts.fill("jet_pfcand_trackjet_d3dsig",
0536                std::isnan(fabs(meas_ip3d.significance())) ? 0 : fabs(meas_ip3d.significance()));
0537       fts.fill("jet_pfcand_trackjet_dist", std::isnan(-meas_jetdist.value()) ? 0 : -meas_jetdist.value());
0538       fts.fill("jet_pfcand_trackjet_decayL", std::isnan(meas_decayl.value()) ? 0 : meas_decayl.value());
0539     } else {
0540       fts.fill("jet_pfcand_dzsig", 0);
0541       fts.fill("jet_pfcand_dxysig", 0);
0542       fts.fill("jet_pfcand_track_chi2", 0);
0543       fts.fill("jet_pfcand_track_qual", 0);
0544       fts.fill("jet_pfcand_trackjet_d3d", 0);
0545       fts.fill("jet_pfcand_trackjet_d3dsig", 0);
0546       fts.fill("jet_pfcand_trackjet_dist", 0);
0547       fts.fill("jet_pfcand_trackjet_decayL", 0);
0548     }
0549 
0550     // muons specific
0551     if (abs(cand->pdgId()) == 13) {
0552       std::vector<unsigned int> muonsToSkip;
0553       int ipos = -1;
0554       float minDR = 1000;
0555       for (size_t i = 0; i < muons.size(); i++) {
0556         if (not muons[i].isPFMuon())
0557           continue;
0558         if (std::find(muonsToSkip.begin(), muonsToSkip.end(), i) != muonsToSkip.end())
0559           continue;
0560         float dR = reco::deltaR(muons[i].p4(), candP4);
0561         if (dR < jet_radius_ and dR < minDR) {
0562           minDR = dR;
0563           ipos = i;
0564           muonsToSkip.push_back(i);
0565         }
0566       }
0567       if (ipos >= 0) {
0568         int muonId = 0;
0569         if (muons[ipos].passed(reco::Muon::CutBasedIdLoose))
0570           muonId++;
0571         if (muons[ipos].passed(reco::Muon::CutBasedIdMedium))
0572           muonId++;
0573         if (muons[ipos].passed(reco::Muon::CutBasedIdTight))
0574           muonId++;
0575         if (muons[ipos].passed(reco::Muon::CutBasedIdGlobalHighPt))
0576           muonId++;
0577         if (muons[ipos].passed(reco::Muon::CutBasedIdTrkHighPt))
0578           muonId++;
0579         fts.fill("jet_pfcand_muon_id", muonId);
0580         fts.fill("jet_pfcand_muon_isglobal", muons[ipos].isGlobalMuon());
0581         fts.fill("jet_pfcand_muon_chi2",
0582                  (muons[ipos].isGlobalMuon()) ? muons[ipos].globalTrack()->normalizedChi2() : 0);
0583         fts.fill("jet_pfcand_muon_nvalidhit",
0584                  (muons[ipos].isGlobalMuon()) ? muons[ipos].globalTrack()->hitPattern().numberOfValidMuonHits() : 0);
0585         fts.fill("jet_pfcand_muon_nstation", muons[ipos].numberOfMatchedStations());
0586         fts.fill("jet_pfcand_muon_segcomp", muon::segmentCompatibility(muons[ipos]));
0587       } else {
0588         fts.fill("jet_pfcand_muon_id", 0);
0589         fts.fill("jet_pfcand_muon_isglobal", 0);
0590         fts.fill("jet_pfcand_muon_chi2", 0);
0591         fts.fill("jet_pfcand_muon_nvalidhit", 0);
0592         fts.fill("jet_pfcand_muon_nstation", 0);
0593         fts.fill("jet_pfcand_muon_segcomp", 0);
0594       }
0595     } else {
0596       fts.fill("jet_pfcand_muon_id", 0);
0597       fts.fill("jet_pfcand_muon_isglobal", 0);
0598       fts.fill("jet_pfcand_muon_chi2", 0);
0599       fts.fill("jet_pfcand_muon_nvalidhit", 0);
0600       fts.fill("jet_pfcand_muon_nstation", 0);
0601       fts.fill("jet_pfcand_muon_segcomp", 0);
0602     }
0603 
0604     // electrons specific
0605     if (abs(cand->pdgId()) == 11) {
0606       int ipos = -1;
0607       for (size_t i = 0; i < electrons.size(); i++) {
0608         if (electrons[i].isPF()) {
0609           for (const auto &element : electrons[i].associatedPackedPFCandidates()) {
0610             if (abs(element->pdgId()) == 11 and element->p4() == candP4)
0611               ipos = i;
0612           }
0613         }
0614       }
0615       if (ipos >= 0) {
0616         fts.fill("jet_pfcand_electron_detaIn",
0617                  std::isnan(electrons[ipos].deltaEtaSuperClusterTrackAtVtx())
0618                      ? 0
0619                      : electrons[ipos].deltaEtaSuperClusterTrackAtVtx());
0620         fts.fill("jet_pfcand_electron_dphiIn",
0621                  std::isnan(electrons[ipos].deltaPhiSuperClusterTrackAtVtx())
0622                      ? 0
0623                      : electrons[ipos].deltaPhiSuperClusterTrackAtVtx());
0624         fts.fill("jet_pfcand_electron_sigIetaIeta",
0625                  std::isnan(electrons[ipos].full5x5_sigmaIetaIeta()) ? 0 : electrons[ipos].full5x5_sigmaIetaIeta());
0626         fts.fill("jet_pfcand_electron_sigIphiIphi",
0627                  std::isnan(electrons[ipos].full5x5_sigmaIphiIphi()) ? 0 : electrons[ipos].full5x5_sigmaIphiIphi());
0628         fts.fill("jet_pfcand_electron_r9", std::isnan(electrons[ipos].full5x5_r9()) ? 0 : electrons[ipos].full5x5_r9());
0629         fts.fill("jet_pfcand_electron_convProb",
0630                  std::isnan(electrons[ipos].convVtxFitProb()) ? 0 : electrons[ipos].convVtxFitProb());
0631       } else {
0632         fts.fill("jet_pfcand_electron_detaIn", 0);
0633         fts.fill("jet_pfcand_electron_dphiIn", 0);
0634         fts.fill("jet_pfcand_electron_sigIetaIeta", 0);
0635         fts.fill("jet_pfcand_electron_sigIphiIphi", 0);
0636         fts.fill("jet_pfcand_electron_r9", 0);
0637         fts.fill("jet_pfcand_electron_convProb", 0);
0638       }
0639     } else {
0640       fts.fill("jet_pfcand_electron_detaIn", 0);
0641       fts.fill("jet_pfcand_electron_dphiIn", 0);
0642       fts.fill("jet_pfcand_electron_sigIetaIeta", 0);
0643       fts.fill("jet_pfcand_electron_sigIphiIphi", 0);
0644       fts.fill("jet_pfcand_electron_r9", 0);
0645       fts.fill("jet_pfcand_electron_convProb", 0);
0646     }
0647 
0648     // photons specific
0649     if (abs(cand->pdgId()) == 22) {
0650       int ipos = -1;
0651       for (size_t i = 0; i < photons.size(); i++) {
0652         for (const auto &element : photons[i].associatedPackedPFCandidates()) {
0653           if (abs(element->pdgId()) == 22 and element->p4() == candP4)
0654             ipos = i;
0655         }
0656       }
0657       if (ipos >= 0) {
0658         fts.fill("jet_pfcand_photon_sigIetaIeta",
0659                  std::isnan(photons[ipos].full5x5_sigmaIetaIeta()) ? 0 : photons[ipos].full5x5_sigmaIetaIeta());
0660         fts.fill("jet_pfcand_photon_r9", std::isnan(photons[ipos].full5x5_r9()) ? 0 : photons[ipos].full5x5_r9());
0661         fts.fill("jet_pfcand_photon_eVeto", photons[ipos].passElectronVeto());
0662       } else {
0663         fts.fill("jet_pfcand_photon_sigIetaIeta", 0);
0664         fts.fill("jet_pfcand_photon_r9", 0);
0665         fts.fill("jet_pfcand_photon_eVeto", 0);
0666       }
0667     } else {
0668       fts.fill("jet_pfcand_photon_sigIetaIeta", 0);
0669       fts.fill("jet_pfcand_photon_r9", 0);
0670       fts.fill("jet_pfcand_photon_eVeto", 0);
0671     }
0672 
0673     // tau specific prior to any puppi weight application
0674     if (std::find(tau_pfcandidates.begin(), tau_pfcandidates.end(), cand->p4()) != tau_pfcandidates.end())
0675       fts.fill("jet_pfcand_tau_signal", 1);
0676     else
0677       fts.fill("jet_pfcand_tau_signal", 0);
0678   }
0679 }
0680 
0681 void ParticleNetFeatureEvaluator::fillSVFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet) {
0682   // secondary vertexes matching 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 
0690   // sort by dxy significance
0691   std::sort(jetSVs.begin(),
0692             jetSVs.end(),
0693             [&](const reco::VertexCompositePtrCandidate *sva, const reco::VertexCompositePtrCandidate *svb) {
0694               return sv_vertex_comparator(*sva, *svb, *pv_);
0695             });
0696 
0697   // reserve space
0698   for (const auto &name : sv_features_)
0699     fts.reserve(name, jetSVs.size());
0700 
0701   GlobalVector jet_global_vec(jet.px(), jet.py(), jet.pz());
0702 
0703   float ip_sign = flip_ip_sign_ ? -1.f : 1.f;
0704 
0705   for (const auto *sv : jetSVs) {
0706     fts.fill("sv_mask", 1);
0707     fts.fill("jet_sv_pt_log", std::isnan(std::log(sv->pt())) ? 0 : std::log(sv->pt()));
0708     fts.fill("jet_sv_eta", sv->eta());
0709     fts.fill("jet_sv_mass", sv->mass());
0710     fts.fill("jet_sv_deta", sv->eta() - jet.eta());
0711     fts.fill("jet_sv_dphi", sv->phi() - jet.phi());
0712     fts.fill("jet_sv_ntrack", sv->numberOfDaughters());
0713     fts.fill("jet_sv_chi2", sv->vertexNormalizedChi2());
0714 
0715     reco::Vertex::CovarianceMatrix csv;
0716     sv->fillVertexCovariance(csv);
0717     reco::Vertex svtx(sv->vertex(), csv);
0718 
0719     VertexDistanceXY dxy;
0720     auto valxy = dxy.signedDistance(svtx, *pv_, jet_global_vec);
0721     fts.fill("jet_sv_dxy", ip_sign * (std::isnan(valxy.value()) ? 0 : valxy.value()));
0722     fts.fill("jet_sv_dxysig", std::isnan(fabs(valxy.significance())) ? 0 : fabs(valxy.significance()));
0723 
0724     VertexDistance3D d3d;
0725     auto val3d = d3d.signedDistance(svtx, *pv_, jet_global_vec);
0726     fts.fill("jet_sv_d3d", ip_sign * (std::isnan(val3d.value()) ? 0 : val3d.value()));
0727     fts.fill("jet_sv_d3dsig", std::isnan(fabs(val3d.significance())) ? 0 : fabs(val3d.significance()));
0728   }
0729 }
0730 
0731 void ParticleNetFeatureEvaluator::fillLostTrackFeatures(DeepBoostedJetFeatures &fts, const reco::Jet &jet) {
0732   // some jet properties
0733   TVector3 jet_direction(jet.momentum().Unit().x(), jet.momentum().Unit().y(), jet.momentum().Unit().z());
0734   GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz());
0735   math::XYZVector jet_dir = jet.momentum().Unit();
0736 
0737   std::vector<pat::PackedCandidate> jet_lost_tracks;
0738   for (size_t itrk = 0; itrk < losttracks_->size(); itrk++) {
0739     if (reco::deltaR(losttracks_->at(itrk).p4(), jet.p4()) < max_dr_for_losttrack_ and
0740         losttracks_->at(itrk).pt() > min_pt_for_losttrack_) {
0741       jet_lost_tracks.push_back(losttracks_->at(itrk));
0742     }
0743   }
0744   std::sort(
0745       jet_lost_tracks.begin(), jet_lost_tracks.end(), [](const auto &a, const auto &b) { return a.pt() > b.pt(); });
0746 
0747   // reserve space
0748   for (const auto &name : losttrack_features_)
0749     fts.reserve(name, jet_lost_tracks.size());
0750 
0751   reco::VertexRef pv_ass = reco::VertexRef(vtxs_, 0);
0752   math::XYZPoint pv_ass_pos = pv_ass->position();
0753 
0754   for (auto const &ltrack : jet_lost_tracks) {
0755     const reco::Track *track = ltrack.bestTrack();
0756 
0757     // only when computing the nagative tagger: remove charged candidates with high sip3d
0758     if (flip_ip_sign_) {
0759       reco::TransientTrack transientTrack = track_builder_->build(*track);
0760       Measurement1D meas_ip3d = IPTools::signedImpactParameter3D(transientTrack, jet_ref_track_dir, *pv_).second;
0761       if (!std::isnan(meas_ip3d.significance()) && meas_ip3d.significance() > max_sip3dsig_for_flip_) {
0762         continue;
0763       }
0764     }
0765     float ip_sign = flip_ip_sign_ ? -1.f : 1.f;
0766 
0767     fts.fill("jet_losttrack_pt_log", std::isnan(std::log(ltrack.pt())) ? 0 : std::log(ltrack.pt()));
0768     fts.fill("jet_losttrack_eta", ltrack.eta());
0769     fts.fill("jet_losttrack_charge", ltrack.charge());
0770     fts.fill("jet_losttrack_frompv", ltrack.fromPV());
0771     fts.fill("jet_losttrack_dz", ip_sign * (std::isnan(ltrack.dz(pv_ass_pos)) ? 0 : ltrack.dz(pv_ass_pos)));
0772     fts.fill("jet_losttrack_dxy", ip_sign * (std::isnan(ltrack.dxy(pv_ass_pos)) ? 0 : ltrack.dxy(pv_ass_pos)));
0773     fts.fill("jet_losttrack_npixhits", ltrack.numberOfPixelHits());
0774     fts.fill("jet_losttrack_nstriphits", ltrack.stripLayersWithMeasurement());
0775 
0776     TVector3 ltrack_momentum(ltrack.momentum().x(), ltrack.momentum().y(), ltrack.momentum().z());
0777     fts.fill("jet_losttrack_deta", jet_direction.Eta() - ltrack_momentum.Eta());
0778     fts.fill("jet_losttrack_dphi", jet_direction.DeltaPhi(ltrack_momentum));
0779     fts.fill("jet_losttrack_etarel",
0780              std::isnan(reco::btau::etaRel(jet_dir, ltrack.momentum()))
0781                  ? 0
0782                  : reco::btau::etaRel(jet_dir, ltrack.momentum()));
0783 
0784     if (track) {
0785       fts.fill("jet_losttrack_track_chi2", track->normalizedChi2());
0786       fts.fill("jet_losttrack_track_qual", track->qualityMask());
0787       fts.fill("jet_losttrack_dxysig",
0788                std::isnan(fabs(ltrack.dxy(pv_ass_pos)) / ltrack.dxyError())
0789                    ? 0
0790                    : fabs(ltrack.dxy(pv_ass_pos)) / ltrack.dxyError());
0791       fts.fill("jet_losttrack_dzsig",
0792                std::isnan(fabs(ltrack.dz(pv_ass_pos)) / ltrack.dzError())
0793                    ? 0
0794                    : fabs(ltrack.dz(pv_ass_pos)) / ltrack.dzError());
0795 
0796       reco::TransientTrack transientTrack = track_builder_->build(*track);
0797       Measurement1D meas_ip3d = IPTools::signedImpactParameter3D(transientTrack, jet_ref_track_dir, *pv_).second;
0798       Measurement1D meas_jetdist = IPTools::jetTrackDistance(transientTrack, jet_ref_track_dir, *pv_).second;
0799       Measurement1D meas_decayl = IPTools::signedDecayLength3D(transientTrack, jet_ref_track_dir, *pv_).second;
0800 
0801       fts.fill("jet_losttrack_trackjet_d3d", ip_sign * (std::isnan(meas_ip3d.value()) ? 0 : meas_ip3d.value()));
0802       fts.fill("jet_losttrack_trackjet_d3dsig",
0803                std::isnan(fabs(meas_ip3d.significance())) ? 0 : fabs(meas_ip3d.significance()));
0804       fts.fill("jet_losttrack_trackjet_dist", std::isnan(-meas_jetdist.value()) ? 0 : -meas_jetdist.value());
0805       fts.fill("jet_losttrack_trackjet_decayL", std::isnan(meas_decayl.value()) ? 0 : meas_decayl.value());
0806     } else {
0807       fts.fill("jet_losttrack_track_chi2", 0);
0808       fts.fill("jet_losttrack_track_qual", 0);
0809       fts.fill("jet_losttrack_dxysig", 0);
0810       fts.fill("jet_losttrack_dzsig", 0);
0811       fts.fill("jet_losttrack_trackjet_d3d", 0);
0812       fts.fill("jet_losttrack_trackjet_d3dsig", 0);
0813       fts.fill("jet_losttrack_trackjet_dist", 0);
0814       fts.fill("jet_losttrack_trackjet_decayL", 0);
0815     }
0816 
0817     fts.fill("lt_mask", 1);
0818   }
0819 }
0820 
0821 // define this as a plug-in
0822 DEFINE_FWK_MODULE(ParticleNetFeatureEvaluator);