FeatureTools/plugins/ParticleTransformerAK4TagInfoProducer.cc

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 #include "FWCore/Utilities/interface/ESGetToken.h"
0010
0011 #include "DataFormats/PatCandidates/interface/Jet.h"
0012 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0013
0014 #include "DataFormats/BTauReco/interface/ShallowTagInfo.h"
0015
0016 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0017 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0018
0019 #include "DataFormats/BTauReco/interface/ParticleTransformerAK4TagInfo.h"
0020 #include "DataFormats/BTauReco/interface/ParticleTransformerAK4Features.h"
0021
0022 #include "RecoBTag/FeatureTools/interface/JetConverter.h"
0023 #include "RecoBTag/FeatureTools/interface/ShallowTagInfoConverter.h"
0024 #include "RecoBTag/FeatureTools/interface/SecondaryVertexConverter.h"
0025 #include "RecoBTag/FeatureTools/interface/NeutralCandidateConverter.h"
0026 #include "RecoBTag/FeatureTools/interface/ChargedCandidateConverter.h"
0027
0028 #include "RecoBTag/FeatureTools/interface/TrackInfoBuilder.h"
0029 #include "RecoBTag/FeatureTools/interface/sorting_modules.h"
0030
0031 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0032 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0033
0034 #include "RecoBTag/FeatureTools/interface/deep_helpers.h"
0035
0036 #include "FWCore/ParameterSet/interface/Registry.h"
0037 #include "FWCore/Common/interface/Provenance.h"
0038 #include "DataFormats/Provenance/interface/ProductProvenance.h"
0039
0040 #include "DataFormats/BTauReco/interface/SeedingTrackFeatures.h"
0041 #include "DataFormats/BTauReco/interface/TrackPairFeatures.h"
0042 #include "RecoBTag/FeatureTools/interface/TrackPairInfoBuilder.h"
0043 #include "RecoBTag/FeatureTools/interface/SeedingTrackInfoBuilder.h"
0044 #include "RecoBTag/FeatureTools/interface/SeedingTracksConverter.h"
0045
0046 #include "FWCore/Framework/interface/ESHandle.h"
0047 #include "FWCore/Framework/interface/EventSetup.h"
0048 #include "RecoBTag/TrackProbability/interface/HistogramProbabilityEstimator.h"
0049 class HistogramProbabilityEstimator;
0050 #include <memory>
0051
0052 #include <typeinfo>
0053 #include "CondFormats/BTauObjects/interface/TrackProbabilityCalibration.h"
0054 #include "CondFormats/DataRecord/interface/BTagTrackProbability2DRcd.h"
0055 #include "CondFormats/DataRecord/interface/BTagTrackProbability3DRcd.h"
0056 #include "FWCore/Framework/interface/EventSetupRecord.h"
0057 #include "FWCore/Framework/interface/EventSetupRecordImplementation.h"
0058 #include "FWCore/Framework/interface/EventSetupRecordKey.h"
0059 #include "DataFormats/Common/interface/AssociationMap.h"
0060
0061 class ParticleTransformerAK4TagInfoProducer : public edm::stream::EDProducer<> {
0062 public:
0063   explicit ParticleTransformerAK4TagInfoProducer(const edm::ParameterSet&);
0064   ~ParticleTransformerAK4TagInfoProducer() override = default;
0065
0066   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0067
0068 private:
0069   typedef std::vector<reco::ParticleTransformerAK4TagInfo> ParticleTransformerAK4TagInfoCollection;
0070   typedef reco::VertexCompositePtrCandidateCollection SVCollection;
0071   typedef reco::VertexCollection VertexCollection;
0072   typedef edm::AssociationMap<edm::OneToOne<reco::JetView, reco::JetView>> JetMatchMap;
0073
0074   void beginStream(edm::StreamID) override {}
0075   void produce(edm::Event&, const edm::EventSetup&) override;
0076   void endStream() override {}
0077
0078   const double jet_radius_;
0079   const double min_candidate_pt_;
0080   const bool flip_;
0081   const double max_sip3dsig_for_flip_;
0082
0083   const edm::EDGetTokenT<edm::View<reco::Jet>> jet_token_;
0084   const edm::EDGetTokenT<VertexCollection> vtx_token_;
0085   const edm::EDGetTokenT<SVCollection> sv_token_;
0086   edm::EDGetTokenT<JetMatchMap> unsubjet_map_token_;
0087   edm::EDGetTokenT<edm::ValueMap<float>> puppi_value_map_token_;
0088   edm::EDGetTokenT<edm::ValueMap<int>> pvasq_value_map_token_;
0089   edm::EDGetTokenT<edm::Association<VertexCollection>> pvas_token_;
0090   const edm::EDGetTokenT<edm::View<reco::Candidate>> candidateToken_;
0091   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> track_builder_token_;
0092   bool use_puppi_value_map_;
0093   bool use_pvasq_value_map_;
0094   bool use_unsubjet_map_;
0095
0096   const bool fallback_puppi_weight_;
0097   const bool fallback_vertex_association_;
0098
0099   const bool is_weighted_jet_;
0100
0101   const double min_jet_pt_;
0102   const double max_jet_eta_;
0103 };
0104
0105 ParticleTransformerAK4TagInfoProducer::ParticleTransformerAK4TagInfoProducer(const edm::ParameterSet& iConfig)
0106     : jet_radius_(iConfig.getParameter<double>("jet_radius")),
0107       min_candidate_pt_(iConfig.getParameter<double>("min_candidate_pt")),
0108       flip_(iConfig.getParameter<bool>("flip")),
0109       max_sip3dsig_for_flip_(iConfig.getParameter<double>("max_sip3dsig_for_flip")),
0110       jet_token_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"))),
0111       vtx_token_(consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0112       sv_token_(consumes<SVCollection>(iConfig.getParameter<edm::InputTag>("secondary_vertices"))),
0113       candidateToken_(consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("candidates"))),
0114       track_builder_token_(
0115           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
0116       use_puppi_value_map_(false),
0117       use_pvasq_value_map_(false),
0118       use_unsubjet_map_(false),
0119       fallback_puppi_weight_(iConfig.getParameter<bool>("fallback_puppi_weight")),
0120       fallback_vertex_association_(iConfig.getParameter<bool>("fallback_vertex_association")),
0121       is_weighted_jet_(iConfig.getParameter<bool>("is_weighted_jet")),
0122       min_jet_pt_(iConfig.getParameter<double>("min_jet_pt")),
0123       max_jet_eta_(iConfig.getParameter<double>("max_jet_eta")) {
0124   produces<ParticleTransformerAK4TagInfoCollection>();
0125
0126   const auto& puppi_value_map_tag = iConfig.getParameter<edm::InputTag>("puppi_value_map");
0127   if (!puppi_value_map_tag.label().empty()) {
0128     puppi_value_map_token_ = consumes<edm::ValueMap<float>>(puppi_value_map_tag);
0129     use_puppi_value_map_ = true;
0130   } else if (is_weighted_jet_) {
0131     throw edm::Exception(edm::errors::Configuration,
0132                          "puppi_value_map is not set but jet is weighted. Must set puppi_value_map.");
0133   }
0134
0135   const auto& pvas_tag = iConfig.getParameter<edm::InputTag>("vertex_associator");
0136   if (!pvas_tag.label().empty()) {
0137     pvasq_value_map_token_ = consumes<edm::ValueMap<int>>(pvas_tag);
0138     pvas_token_ = consumes<edm::Association<VertexCollection>>(pvas_tag);
0139     use_pvasq_value_map_ = true;
0140   }
0141
0142   const auto& unsubjet_map_tag = iConfig.getParameter<edm::InputTag>("unsubjet_map");
0143   if (!unsubjet_map_tag.label().empty()) {
0144     unsubjet_map_token_ = consumes<JetMatchMap>(unsubjet_map_tag);
0145     use_unsubjet_map_ = true;
0146   }
0147 }
0148
0149 void ParticleTransformerAK4TagInfoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0150   // pfParticleTransformerAK4TagInfos
0151   edm::ParameterSetDescription desc;
0152   desc.add<double>("jet_radius", 0.4);
0153   desc.add<double>("min_candidate_pt", 0.95);
0154   desc.add<bool>("flip", false);
0155   desc.add<double>("max_sip3dsig_for_flip", 99999);
0156   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0157   desc.add<edm::InputTag>("puppi_value_map", edm::InputTag("puppi"));
0158   desc.add<edm::InputTag>("secondary_vertices", edm::InputTag("inclusiveCandidateSecondaryVertices"));
0159   desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
0160   desc.add<edm::InputTag>("unsubjet_map", {});
0161   desc.add<edm::InputTag>("candidates", edm::InputTag("packedPFCandidates"));
0162   desc.add<edm::InputTag>("vertex_associator", edm::InputTag("primaryVertexAssociation", "original"));
0163   desc.add<bool>("fallback_puppi_weight", false);
0164   desc.add<bool>("fallback_vertex_association", false);
0165   desc.add<bool>("is_weighted_jet", false);
0166   desc.add<double>("min_jet_pt", 15.0);
0167   desc.add<double>("max_jet_eta", 2.5);
0168   descriptions.add("pfParticleTransformerAK4TagInfos", desc);
0169 }
0170
0171 void ParticleTransformerAK4TagInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0172   auto output_tag_infos = std::make_unique<ParticleTransformerAK4TagInfoCollection>();
0173   edm::Handle<edm::View<reco::Jet>> jets;
0174   iEvent.getByToken(jet_token_, jets);
0175
0176   edm::Handle<JetMatchMap> unsubjet_map;
0177   if (use_unsubjet_map_)
0178     iEvent.getByToken(unsubjet_map_token_, unsubjet_map);
0179
0180   edm::Handle<VertexCollection> vtxs;
0181   iEvent.getByToken(vtx_token_, vtxs);
0182   if (vtxs->empty()) {
0183     // produce empty TagInfos in case no primary vertex
0184     iEvent.put(std::move(output_tag_infos));
0185     return;  // exit event
0186   }
0187   // reference to primary vertex
0188   const auto& pv = vtxs->at(0);
0189
0190   edm::Handle<edm::View<reco::Candidate>> tracks;
0191   iEvent.getByToken(candidateToken_, tracks);
0192
0193   edm::Handle<SVCollection> svs;
0194   iEvent.getByToken(sv_token_, svs);
0195
0196   edm::Handle<edm::ValueMap<float>> puppi_value_map;
0197   if (use_puppi_value_map_) {
0198     iEvent.getByToken(puppi_value_map_token_, puppi_value_map);
0199   }
0200
0201   edm::Handle<edm::ValueMap<int>> pvasq_value_map;
0202   edm::Handle<edm::Association<VertexCollection>> pvas;
0203   if (use_pvasq_value_map_) {
0204     iEvent.getByToken(pvasq_value_map_token_, pvasq_value_map);
0205     iEvent.getByToken(pvas_token_, pvas);
0206   }
0207
0208   edm::ESHandle<TransientTrackBuilder> track_builder = iSetup.getHandle(track_builder_token_);
0209
0210   for (std::size_t jet_n = 0; jet_n < jets->size(); jet_n++) {
0211     // create data containing structure
0212     btagbtvdeep::ParticleTransformerAK4Features features;
0213
0214     // reco jet reference (use as much as possible)
0215     const auto& jet = jets->at(jet_n);
0216     if (jet.pt() < 15.0) {
0217       features.is_filled = false;
0218     }
0219     if (std::abs(jet.eta()) > 2.5) {
0220       features.is_filled = false;
0221     }
0222     // dynamical casting to pointers, null if not possible
0223     const auto* pf_jet = dynamic_cast<const reco::PFJet*>(&jet);
0224     const auto* pat_jet = dynamic_cast<const pat::Jet*>(&jet);
0225     edm::RefToBase<reco::Jet> jet_ref(jets, jet_n);
0226     const auto& unsubJet =
0227         (use_unsubjet_map_ && (*unsubjet_map)[jet_ref].isNonnull()) ? *(*unsubjet_map)[jet_ref] : jet;
0228
0229     if (features.is_filled) {
0230       math::XYZVector jet_dir = jet.momentum().Unit();
0231       GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz());
0232
0233       // copy which will be sorted
0234       auto svs_sorted = *svs;
0235       // sort by dxy
0236       std::sort(svs_sorted.begin(), svs_sorted.end(), [&pv](const auto& sva, const auto& svb) {
0237         return btagbtvdeep::sv_vertex_comparator(sva, svb, pv);
0238       });
0239       // fill features from secondary vertices
0240       for (const auto& sv : svs_sorted) {
0241         if (reco::deltaR2(sv, jet_dir) > (jet_radius_ * jet_radius_))
0242           continue;
0243         else {
0244           features.sv_features.emplace_back();
0245           // in C++17 could just get from emplace_back output
0246           auto& sv_features = features.sv_features.back();
0247           btagbtvdeep::svToFeatures(sv, pv, jet, sv_features, flip_);
0248         }
0249       }
0250
0251       // stuff required for dealing with pf candidates
0252
0253       std::vector<btagbtvdeep::SortingClass<size_t>> c_sorted, n_sorted;
0254
0255       // to cache the TrackInfo
0256       std::map<unsigned int, btagbtvdeep::TrackInfoBuilder> trackinfos;
0257
0258       // unsorted reference to sv
0259       const auto& svs_unsorted = *svs;
0260       // fill collection, from DeepTNtuples plus some styling
0261       for (unsigned int i = 0; i < unsubJet.numberOfDaughters(); i++) {
0262         auto cand = unsubJet.daughter(i);
0263         if (cand) {
0264           // candidates under 950MeV (configurable) are not considered
0265           // might change if we use also white-listing
0266           if (cand->pt() < min_candidate_pt_)
0267             continue;
0268           if (cand->charge() != 0) {
0269             auto& trackinfo = trackinfos.emplace(i, track_builder).first->second;
0270             trackinfo.buildTrackInfo(cand, jet_dir, jet_ref_track_dir, pv);
0271
0272             c_sorted.emplace_back(i,
0273                                   trackinfo.getTrackSip2dSig(),
0274                                   -btagbtvdeep::mindrsvpfcand(svs_unsorted, cand),
0275                                   cand->pt() / jet.pt());
0276           } else {
0277             n_sorted.emplace_back(i, -1, -btagbtvdeep::mindrsvpfcand(svs_unsorted, cand), cand->pt() / jet.pt());
0278           }
0279         }
0280       }
0281
0282       // sort collections (open the black-box if you please)
0283       std::sort(c_sorted.begin(), c_sorted.end(), btagbtvdeep::SortingClass<std::size_t>::compareByABCInv);
0284       std::sort(n_sorted.begin(), n_sorted.end(), btagbtvdeep::SortingClass<std::size_t>::compareByABCInv);
0285
0286       std::vector<size_t> c_sortedindices, n_sortedindices;
0287
0288       // this puts 0 everywhere and the right position in ind
0289       c_sortedindices = btagbtvdeep::invertSortingVector(c_sorted);
0290       n_sortedindices = btagbtvdeep::invertSortingVector(n_sorted);
0291
0292       // set right size to vectors
0293       features.c_pf_features.clear();
0294       features.c_pf_features.resize(c_sorted.size());
0295       features.n_pf_features.clear();
0296       features.n_pf_features.resize(n_sorted.size());
0297
0298       for (unsigned int i = 0; i < unsubJet.numberOfDaughters(); i++) {
0299         // get pointer and check that is correct
0300         auto cand = dynamic_cast<const reco::Candidate*>(unsubJet.daughter(i));
0301         if (!cand)
0302           continue;
0303         // candidates under 950MeV are not considered
0304         // might change if we use also white-listing
0305         if (cand->pt() < 0.95)
0306           continue;
0307
0308         auto packed_cand = dynamic_cast<const pat::PackedCandidate*>(cand);
0309         auto reco_cand = dynamic_cast<const reco::PFCandidate*>(cand);
0310
0311         // need some edm::Ptr or edm::Ref if reco candidates
0312         reco::PFCandidatePtr reco_ptr;
0313         if (pf_jet) {
0314           reco_ptr = pf_jet->getPFConstituent(i);
0315         } else if (pat_jet && reco_cand) {
0316           reco_ptr = pat_jet->getPFConstituent(i);
0317         }
0318
0319         reco::CandidatePtr cand_ptr;
0320         if (pat_jet) {
0321           cand_ptr = pat_jet->sourceCandidatePtr(i);
0322         }
0323
0324         //
0325         // Access puppi weight from ValueMap.
0326         //
0327         float puppiw = 1.0;  // Set to fallback value
0328
0329         if (reco_cand) {
0330           if (use_puppi_value_map_)
0331             puppiw = (*puppi_value_map)[reco_ptr];
0332           else if (!fallback_puppi_weight_) {
0333             throw edm::Exception(edm::errors::InvalidReference, "PUPPI value map missing")
0334                 << "use fallback_puppi_weight option to use " << puppiw << " for reco_cand as default";
0335           }
0336         } else if (packed_cand) {
0337           if (use_puppi_value_map_)
0338             puppiw = (*puppi_value_map)[cand_ptr];
0339           else if (!fallback_puppi_weight_) {
0340             throw edm::Exception(edm::errors::InvalidReference, "PUPPI value map missing")
0341                 << "use fallback_puppi_weight option to use " << puppiw << " for packed_cand as default";
0342           }
0343         } else {
0344           throw edm::Exception(edm::errors::InvalidReference)
0345               << "Cannot convert to either reco::PFCandidate or pat::PackedCandidate";
0346         }
0347
0348         float drminpfcandsv = btagbtvdeep::mindrsvpfcand(svs_unsorted, cand);
0349         float distminpfcandsv = 0;
0350         if (cand->charge() != 0) {
0351           // is charged candidate
0352           auto entry = c_sortedindices.at(i);
0353
0354           // get cached track info
0355           auto& trackinfo = trackinfos.at(i);
0356           if (flip_ && (trackinfo.getTrackSip3dSig() > max_sip3dsig_for_flip_)) {
0357             continue;
0358           }
0359           // get_ref to vector element
0360           auto& c_pf_features = features.c_pf_features.at(entry);
0361           // fill feature structure
0362           if (packed_cand) {
0363             if (packed_cand->hasTrackDetails()) {
0364               const reco::Track& PseudoTrack = packed_cand->pseudoTrack();
0365               reco::TransientTrack transientTrack;
0366               transientTrack = track_builder->build(PseudoTrack);
0367               distminpfcandsv = btagbtvdeep::mindistsvpfcand(svs_unsorted, transientTrack);
0368             }
0369
0370             btagbtvdeep::packedCandidateToFeatures(packed_cand,
0371                                                    jet,
0372                                                    trackinfo,
0373                                                    is_weighted_jet_,
0374                                                    drminpfcandsv,
0375                                                    static_cast<float>(jet_radius_),
0376                                                    puppiw,
0377                                                    c_pf_features,
0378                                                    flip_,
0379                                                    distminpfcandsv);
0380           } else if (reco_cand) {
0381             // get vertex association quality
0382             int pv_ass_quality = 0;  // fallback value
0383             if (use_pvasq_value_map_) {
0384               pv_ass_quality = (*pvasq_value_map)[reco_ptr];
0385             } else if (!fallback_vertex_association_) {
0386               throw edm::Exception(edm::errors::InvalidReference, "vertex association missing")
0387                   << "use fallback_vertex_association option to use" << pv_ass_quality
0388                   << "as default quality and closest dz PV as criteria";
0389             }
0390             // getting the PV as PackedCandidatesProducer
0391             // but using not the slimmed but original vertices
0392             auto ctrack = reco_cand->bestTrack();
0393             int pvi = -1;
0394             float dist = 1e99;
0395             for (size_t ii = 0; ii < vtxs->size(); ii++) {
0396               float dz = (ctrack) ? std::abs(ctrack->dz(((*vtxs)[ii]).position())) : 0;
0397               if (dz < dist) {
0398                 pvi = ii;
0399                 dist = dz;
0400               }
0401             }
0402             auto PV = reco::VertexRef(vtxs, pvi);
0403             if (use_pvasq_value_map_) {
0404               const reco::VertexRef& PV_orig = (*pvas)[reco_ptr];
0405               if (PV_orig.isNonnull())
0406                 PV = reco::VertexRef(vtxs, PV_orig.key());
0407             }
0408             btagbtvdeep::recoCandidateToFeatures(reco_cand,
0409                                                  jet,
0410                                                  trackinfo,
0411                                                  is_weighted_jet_,
0412                                                  drminpfcandsv,
0413                                                  static_cast<float>(jet_radius_),
0414                                                  puppiw,
0415                                                  pv_ass_quality,
0416                                                  PV,
0417                                                  c_pf_features,
0418                                                  flip_,
0419                                                  distminpfcandsv);
0420           }
0421         } else {
0422           // is neutral candidate
0423           auto entry = n_sortedindices.at(i);
0424           // get_ref to vector element
0425           auto& n_pf_features = features.n_pf_features.at(entry);
0426           // fill feature structure
0427           if (packed_cand) {
0428             btagbtvdeep::packedCandidateToFeatures(packed_cand,
0429                                                    jet,
0430                                                    is_weighted_jet_,
0431                                                    drminpfcandsv,
0432                                                    static_cast<float>(jet_radius_),
0433                                                    puppiw,
0434                                                    n_pf_features);
0435           } else if (reco_cand) {
0436             btagbtvdeep::recoCandidateToFeatures(
0437                 reco_cand, jet, is_weighted_jet_, drminpfcandsv, static_cast<float>(jet_radius_), puppiw, n_pf_features);
0438           }
0439         }
0440       }
0441     }
0442     output_tag_infos->emplace_back(features, jet_ref);
0443   }
0444   iEvent.put(std::move(output_tag_infos));
0445 }
0446
0447 //define this as a plug-in
0448 DEFINE_FWK_MODULE(ParticleTransformerAK4TagInfoProducer);