File indexing completed on 2021-02-14 13:33:38
0001
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/stream/EDProducer.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/StreamID.h"
0010 #include "FWCore/Utilities/interface/ESGetToken.h"
0011
0012 #include "DataFormats/PatCandidates/interface/Jet.h"
0013 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0014
0015 #include "DataFormats/BTauReco/interface/ShallowTagInfo.h"
0016
0017 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0018 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0019
0020 #include "DataFormats/BTauReco/interface/DeepFlavourTagInfo.h"
0021 #include "DataFormats/BTauReco/interface/DeepFlavourFeatures.h"
0022
0023 #include "RecoBTag/FeatureTools/interface/JetConverter.h"
0024 #include "RecoBTag/FeatureTools/interface/ShallowTagInfoConverter.h"
0025 #include "RecoBTag/FeatureTools/interface/SecondaryVertexConverter.h"
0026 #include "RecoBTag/FeatureTools/interface/NeutralCandidateConverter.h"
0027 #include "RecoBTag/FeatureTools/interface/ChargedCandidateConverter.h"
0028
0029 #include "RecoBTag/FeatureTools/interface/TrackInfoBuilder.h"
0030 #include "RecoBTag/FeatureTools/interface/sorting_modules.h"
0031
0032 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0033 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0034
0035 #include "RecoBTag/FeatureTools/interface/deep_helpers.h"
0036
0037 #include "FWCore/ParameterSet/interface/Registry.h"
0038 #include "FWCore/Common/interface/Provenance.h"
0039 #include "DataFormats/Provenance/interface/ProductProvenance.h"
0040
0041 #include "DataFormats/BTauReco/interface/SeedingTrackFeatures.h"
0042 #include "DataFormats/BTauReco/interface/TrackPairFeatures.h"
0043 #include "RecoBTag/FeatureTools/interface/TrackPairInfoBuilder.h"
0044 #include "RecoBTag/FeatureTools/interface/SeedingTrackInfoBuilder.h"
0045 #include "RecoBTag/FeatureTools/interface/SeedingTracksConverter.h"
0046
0047 #include "FWCore/Framework/interface/ESHandle.h"
0048 #include "FWCore/Framework/interface/EventSetup.h"
0049 #include "RecoBTag/TrackProbability/interface/HistogramProbabilityEstimator.h"
0050 class HistogramProbabilityEstimator;
0051 #include <memory>
0052
0053 #include <typeinfo>
0054 #include "CondFormats/BTauObjects/interface/TrackProbabilityCalibration.h"
0055 #include "CondFormats/DataRecord/interface/BTagTrackProbability2DRcd.h"
0056 #include "CondFormats/DataRecord/interface/BTagTrackProbability3DRcd.h"
0057 #include "FWCore/Framework/interface/EventSetupRecord.h"
0058 #include "FWCore/Framework/interface/EventSetupRecordImplementation.h"
0059 #include "FWCore/Framework/interface/EventSetupRecordKey.h"
0060
0061 class DeepFlavourTagInfoProducer : public edm::stream::EDProducer<> {
0062 public:
0063 explicit DeepFlavourTagInfoProducer(const edm::ParameterSet&);
0064 ~DeepFlavourTagInfoProducer() override;
0065
0066 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0067
0068 private:
0069 typedef std::vector<reco::DeepFlavourTagInfo> DeepFlavourTagInfoCollection;
0070 typedef reco::VertexCompositePtrCandidateCollection SVCollection;
0071 typedef reco::VertexCollection VertexCollection;
0072 typedef edm::View<reco::ShallowTagInfo> ShallowTagInfoCollection;
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
0082 edm::EDGetTokenT<edm::View<reco::Jet>> jet_token_;
0083 edm::EDGetTokenT<VertexCollection> vtx_token_;
0084 edm::EDGetTokenT<SVCollection> sv_token_;
0085 edm::EDGetTokenT<ShallowTagInfoCollection> shallow_tag_info_token_;
0086 edm::EDGetTokenT<edm::ValueMap<float>> puppi_value_map_token_;
0087 edm::EDGetTokenT<edm::ValueMap<int>> pvasq_value_map_token_;
0088 edm::EDGetTokenT<edm::Association<VertexCollection>> pvas_token_;
0089 edm::EDGetTokenT<edm::View<reco::Candidate>> candidateToken_;
0090 edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> track_builder_token_;
0091 edm::ESGetToken<TrackProbabilityCalibration, BTagTrackProbability2DRcd> calib2d_token_;
0092 edm::ESGetToken<TrackProbabilityCalibration, BTagTrackProbability3DRcd> calib3d_token_;
0093
0094 bool use_puppi_value_map_;
0095 bool use_pvasq_value_map_;
0096
0097 bool fallback_puppi_weight_;
0098 bool fallback_vertex_association_;
0099
0100 bool run_deepVertex_;
0101
0102
0103 void checkEventSetup(const edm::EventSetup& iSetup);
0104 std::unique_ptr<HistogramProbabilityEstimator> probabilityEstimator_;
0105 bool compute_probabilities_;
0106 unsigned long long calibrationCacheId2D_;
0107 unsigned long long calibrationCacheId3D_;
0108
0109 const double min_jet_pt_;
0110 const double max_jet_eta_;
0111 };
0112
0113 DeepFlavourTagInfoProducer::DeepFlavourTagInfoProducer(const edm::ParameterSet& iConfig)
0114 : jet_radius_(iConfig.getParameter<double>("jet_radius")),
0115 min_candidate_pt_(iConfig.getParameter<double>("min_candidate_pt")),
0116 flip_(iConfig.getParameter<bool>("flip")),
0117 jet_token_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"))),
0118 vtx_token_(consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0119 sv_token_(consumes<SVCollection>(iConfig.getParameter<edm::InputTag>("secondary_vertices"))),
0120 shallow_tag_info_token_(
0121 consumes<ShallowTagInfoCollection>(iConfig.getParameter<edm::InputTag>("shallow_tag_infos"))),
0122 candidateToken_(consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("candidates"))),
0123 track_builder_token_(
0124 esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
0125 use_puppi_value_map_(false),
0126 use_pvasq_value_map_(false),
0127 fallback_puppi_weight_(iConfig.getParameter<bool>("fallback_puppi_weight")),
0128 fallback_vertex_association_(iConfig.getParameter<bool>("fallback_vertex_association")),
0129 run_deepVertex_(iConfig.getParameter<bool>("run_deepVertex")),
0130 compute_probabilities_(iConfig.getParameter<bool>("compute_probabilities")),
0131 min_jet_pt_(iConfig.getParameter<double>("min_jet_pt")),
0132 max_jet_eta_(iConfig.getParameter<double>("max_jet_eta")) {
0133 produces<DeepFlavourTagInfoCollection>();
0134
0135 const auto& puppi_value_map_tag = iConfig.getParameter<edm::InputTag>("puppi_value_map");
0136 if (!puppi_value_map_tag.label().empty()) {
0137 puppi_value_map_token_ = consumes<edm::ValueMap<float>>(puppi_value_map_tag);
0138 use_puppi_value_map_ = true;
0139 }
0140
0141 const auto& pvas_tag = iConfig.getParameter<edm::InputTag>("vertex_associator");
0142 if (!pvas_tag.label().empty()) {
0143 pvasq_value_map_token_ = consumes<edm::ValueMap<int>>(pvas_tag);
0144 pvas_token_ = consumes<edm::Association<VertexCollection>>(pvas_tag);
0145 use_pvasq_value_map_ = true;
0146 }
0147 if (compute_probabilities_) {
0148 calib2d_token_ = esConsumes<TrackProbabilityCalibration, BTagTrackProbability2DRcd>();
0149 calib3d_token_ = esConsumes<TrackProbabilityCalibration, BTagTrackProbability3DRcd>();
0150 }
0151 }
0152
0153 DeepFlavourTagInfoProducer::~DeepFlavourTagInfoProducer() {}
0154
0155 void DeepFlavourTagInfoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0156
0157 edm::ParameterSetDescription desc;
0158 desc.add<edm::InputTag>("shallow_tag_infos", edm::InputTag("pfDeepCSVTagInfos"));
0159 desc.add<double>("jet_radius", 0.4);
0160 desc.add<double>("min_candidate_pt", 0.95);
0161 desc.add<bool>("flip", false);
0162 desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0163 desc.add<edm::InputTag>("puppi_value_map", edm::InputTag("puppi"));
0164 desc.add<edm::InputTag>("secondary_vertices", edm::InputTag("inclusiveCandidateSecondaryVertices"));
0165 desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
0166 desc.add<edm::InputTag>("candidates", edm::InputTag("packedPFCandidates"));
0167 desc.add<edm::InputTag>("vertex_associator", edm::InputTag("primaryVertexAssociation", "original"));
0168 desc.add<bool>("fallback_puppi_weight", false);
0169 desc.add<bool>("fallback_vertex_association", false);
0170 desc.add<bool>("run_deepVertex", false);
0171 desc.add<bool>("compute_probabilities", false);
0172 desc.add<double>("min_jet_pt", 15.0);
0173 desc.add<double>("max_jet_eta", 2.5);
0174 descriptions.add("pfDeepFlavourTagInfos", desc);
0175 }
0176
0177 void DeepFlavourTagInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0178 auto output_tag_infos = std::make_unique<DeepFlavourTagInfoCollection>();
0179 if (compute_probabilities_)
0180 checkEventSetup(iSetup);
0181
0182 edm::Handle<edm::View<reco::Jet>> jets;
0183 iEvent.getByToken(jet_token_, jets);
0184
0185 edm::Handle<VertexCollection> vtxs;
0186 iEvent.getByToken(vtx_token_, vtxs);
0187 if (vtxs->empty()) {
0188
0189 iEvent.put(std::move(output_tag_infos));
0190 return;
0191 }
0192
0193 const auto& pv = vtxs->at(0);
0194
0195 edm::Handle<edm::View<reco::Candidate>> tracks;
0196 iEvent.getByToken(candidateToken_, tracks);
0197
0198 edm::Handle<SVCollection> svs;
0199 iEvent.getByToken(sv_token_, svs);
0200
0201 edm::Handle<ShallowTagInfoCollection> shallow_tag_infos;
0202 iEvent.getByToken(shallow_tag_info_token_, shallow_tag_infos);
0203 double negative_cut = 0;
0204 if (flip_) {
0205 const edm::Provenance* prov = shallow_tag_infos.provenance();
0206 const edm::ParameterSet& psetFromProvenance = edm::parameterSet(prov->stable(), iEvent.processHistory());
0207 negative_cut = ((psetFromProvenance.getParameter<edm::ParameterSet>("computer"))
0208 .getParameter<edm::ParameterSet>("trackSelection"))
0209 .getParameter<double>("sip3dSigMax");
0210 }
0211
0212 edm::Handle<edm::ValueMap<float>> puppi_value_map;
0213 if (use_puppi_value_map_) {
0214 iEvent.getByToken(puppi_value_map_token_, puppi_value_map);
0215 }
0216
0217 edm::Handle<edm::ValueMap<int>> pvasq_value_map;
0218 edm::Handle<edm::Association<VertexCollection>> pvas;
0219 if (use_pvasq_value_map_) {
0220 iEvent.getByToken(pvasq_value_map_token_, pvasq_value_map);
0221 iEvent.getByToken(pvas_token_, pvas);
0222 }
0223
0224 edm::ESHandle<TransientTrackBuilder> track_builder = iSetup.getHandle(track_builder_token_);
0225
0226 std::vector<reco::TransientTrack> selectedTracks;
0227 std::vector<float> masses;
0228
0229 if (run_deepVertex_)
0230 {
0231 for (typename edm::View<reco::Candidate>::const_iterator track = tracks->begin(); track != tracks->end(); ++track) {
0232 unsigned int k = track - tracks->begin();
0233 if (track->bestTrack() != nullptr && track->pt() > 0.5) {
0234 if (std::fabs(track->vz() - pv.z()) < 0.5) {
0235 selectedTracks.push_back(track_builder->build(tracks->ptrAt(k)));
0236 masses.push_back(track->mass());
0237 }
0238 }
0239 }
0240 }
0241
0242 for (std::size_t jet_n = 0; jet_n < jets->size(); jet_n++) {
0243
0244 btagbtvdeep::DeepFlavourFeatures features;
0245
0246
0247 const auto& jet = jets->at(jet_n);
0248
0249 const auto* pf_jet = dynamic_cast<const reco::PFJet*>(&jet);
0250 const auto* pat_jet = dynamic_cast<const pat::Jet*>(&jet);
0251 edm::RefToBase<reco::Jet> jet_ref(jets, jet_n);
0252
0253 const edm::View<reco::ShallowTagInfo>& taginfos = *shallow_tag_infos;
0254 edm::Ptr<reco::ShallowTagInfo> match;
0255
0256 if ((jet_n < taginfos.size()) && (taginfos[jet_n].jet() == jet_ref)) {
0257 match = taginfos.ptrAt(jet_n);
0258 } else {
0259
0260 for (auto itTI = taginfos.begin(), edTI = taginfos.end(); itTI != edTI; ++itTI) {
0261 if (itTI->jet() == jet_ref) {
0262 match = taginfos.ptrAt(itTI - taginfos.begin());
0263 break;
0264 }
0265 }
0266 }
0267 reco::ShallowTagInfo tag_info;
0268 if (match.isNonnull()) {
0269 tag_info = *match;
0270 }
0271
0272
0273 btagbtvdeep::JetConverter::jetToFeatures(jet, features.jet_features);
0274
0275
0276 features.npv = vtxs->size();
0277 math::XYZVector jet_dir = jet.momentum().Unit();
0278 GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz());
0279
0280
0281 const auto& tag_info_vars = tag_info.taggingVariables();
0282 btagbtvdeep::bTagToFeatures(tag_info_vars, features.tag_info_features);
0283
0284
0285 auto svs_sorted = *svs;
0286
0287 std::sort(svs_sorted.begin(), svs_sorted.end(), [&pv](const auto& sva, const auto& svb) {
0288 return btagbtvdeep::sv_vertex_comparator(sva, svb, pv);
0289 });
0290
0291 for (const auto& sv : svs_sorted) {
0292 if (reco::deltaR2(sv, jet_dir) > (jet_radius_ * jet_radius_))
0293 continue;
0294 else {
0295 features.sv_features.emplace_back();
0296
0297 auto& sv_features = features.sv_features.back();
0298 btagbtvdeep::svToFeatures(sv, pv, jet, sv_features, flip_);
0299 }
0300 }
0301
0302
0303
0304 std::vector<btagbtvdeep::SortingClass<size_t>> c_sorted, n_sorted;
0305
0306
0307 std::map<unsigned int, btagbtvdeep::TrackInfoBuilder> trackinfos;
0308
0309
0310 const auto& svs_unsorted = *svs;
0311
0312 for (unsigned int i = 0; i < jet.numberOfDaughters(); i++) {
0313 auto cand = jet.daughter(i);
0314 if (cand) {
0315
0316
0317 if (cand->pt() < min_candidate_pt_)
0318 continue;
0319 if (cand->charge() != 0) {
0320 auto& trackinfo = trackinfos.emplace(i, track_builder).first->second;
0321 trackinfo.buildTrackInfo(cand, jet_dir, jet_ref_track_dir, pv);
0322 c_sorted.emplace_back(
0323 i, trackinfo.getTrackSip2dSig(), -btagbtvdeep::mindrsvpfcand(svs_unsorted, cand), cand->pt() / jet.pt());
0324 } else {
0325 n_sorted.emplace_back(i, -1, -btagbtvdeep::mindrsvpfcand(svs_unsorted, cand), cand->pt() / jet.pt());
0326 }
0327 }
0328 }
0329
0330
0331 std::sort(c_sorted.begin(), c_sorted.end(), btagbtvdeep::SortingClass<std::size_t>::compareByABCInv);
0332 std::sort(n_sorted.begin(), n_sorted.end(), btagbtvdeep::SortingClass<std::size_t>::compareByABCInv);
0333
0334 std::vector<size_t> c_sortedindices, n_sortedindices;
0335
0336
0337 c_sortedindices = btagbtvdeep::invertSortingVector(c_sorted);
0338 n_sortedindices = btagbtvdeep::invertSortingVector(n_sorted);
0339
0340
0341 features.c_pf_features.clear();
0342 features.c_pf_features.resize(c_sorted.size());
0343 features.n_pf_features.clear();
0344 features.n_pf_features.resize(n_sorted.size());
0345
0346 for (unsigned int i = 0; i < jet.numberOfDaughters(); i++) {
0347
0348 auto cand = dynamic_cast<const reco::Candidate*>(jet.daughter(i));
0349 if (!cand)
0350 continue;
0351
0352
0353 if (cand->pt() < 0.95)
0354 continue;
0355
0356 auto packed_cand = dynamic_cast<const pat::PackedCandidate*>(cand);
0357 auto reco_cand = dynamic_cast<const reco::PFCandidate*>(cand);
0358
0359
0360 reco::PFCandidatePtr reco_ptr;
0361 if (pf_jet) {
0362 reco_ptr = pf_jet->getPFConstituent(i);
0363 } else if (pat_jet && reco_cand) {
0364 reco_ptr = pat_jet->getPFConstituent(i);
0365 }
0366
0367 float puppiw = 1.0;
0368 if (reco_cand && use_puppi_value_map_) {
0369 puppiw = (*puppi_value_map)[reco_ptr];
0370 } else if (reco_cand && !fallback_puppi_weight_) {
0371 throw edm::Exception(edm::errors::InvalidReference, "PUPPI value map missing")
0372 << "use fallback_puppi_weight option to use " << puppiw << "as default";
0373 }
0374
0375 float drminpfcandsv = btagbtvdeep::mindrsvpfcand(svs_unsorted, cand);
0376
0377 if (cand->charge() != 0) {
0378
0379 auto entry = c_sortedindices.at(i);
0380
0381 auto& trackinfo = trackinfos.at(i);
0382 if (flip_ && (trackinfo.getTrackSip3dSig() > negative_cut)) {
0383 continue;
0384 }
0385
0386 auto& c_pf_features = features.c_pf_features.at(entry);
0387
0388 if (packed_cand) {
0389 btagbtvdeep::packedCandidateToFeatures(
0390 packed_cand, jet, trackinfo, drminpfcandsv, static_cast<float>(jet_radius_), c_pf_features, flip_);
0391 } else if (reco_cand) {
0392
0393 int pv_ass_quality = 0;
0394 if (use_pvasq_value_map_) {
0395 pv_ass_quality = (*pvasq_value_map)[reco_ptr];
0396 } else if (!fallback_vertex_association_) {
0397 throw edm::Exception(edm::errors::InvalidReference, "vertex association missing")
0398 << "use fallback_vertex_association option to use" << pv_ass_quality
0399 << "as default quality and closest dz PV as criteria";
0400 }
0401
0402
0403 auto ctrack = reco_cand->bestTrack();
0404 int pvi = -1;
0405 float dist = 1e99;
0406 for (size_t ii = 0; ii < vtxs->size(); ii++) {
0407 float dz = (ctrack) ? std::abs(ctrack->dz(((*vtxs)[ii]).position())) : 0;
0408 if (dz < dist) {
0409 pvi = ii;
0410 dist = dz;
0411 }
0412 }
0413 auto PV = reco::VertexRef(vtxs, pvi);
0414 if (use_pvasq_value_map_) {
0415 const reco::VertexRef& PV_orig = (*pvas)[reco_ptr];
0416 if (PV_orig.isNonnull())
0417 PV = reco::VertexRef(vtxs, PV_orig.key());
0418 }
0419 btagbtvdeep::recoCandidateToFeatures(reco_cand,
0420 jet,
0421 trackinfo,
0422 drminpfcandsv,
0423 static_cast<float>(jet_radius_),
0424 puppiw,
0425 pv_ass_quality,
0426 PV,
0427 c_pf_features,
0428 flip_);
0429 }
0430 } else {
0431
0432 auto entry = n_sortedindices.at(i);
0433
0434 auto& n_pf_features = features.n_pf_features.at(entry);
0435
0436 if (packed_cand) {
0437 btagbtvdeep::packedCandidateToFeatures(
0438 packed_cand, jet, drminpfcandsv, static_cast<float>(jet_radius_), n_pf_features);
0439 } else if (reco_cand) {
0440 btagbtvdeep::recoCandidateToFeatures(
0441 reco_cand, jet, drminpfcandsv, static_cast<float>(jet_radius_), puppiw, n_pf_features);
0442 }
0443 }
0444 }
0445
0446 if (run_deepVertex_) {
0447 if (jet.pt() > min_jet_pt_ && std::fabs(jet.eta()) < max_jet_eta_)
0448 btagbtvdeep::seedingTracksToFeatures(selectedTracks,
0449 masses,
0450 jet,
0451 pv,
0452 probabilityEstimator_.get(),
0453 compute_probabilities_,
0454 features.seed_features);
0455 }
0456
0457 output_tag_infos->emplace_back(features, jet_ref);
0458 }
0459
0460 iEvent.put(std::move(output_tag_infos));
0461 }
0462
0463 void DeepFlavourTagInfoProducer::checkEventSetup(const edm::EventSetup& iSetup) {
0464 using namespace edm;
0465 using namespace edm::eventsetup;
0466
0467 const EventSetupRecord& re2D = iSetup.get<BTagTrackProbability2DRcd>();
0468 const EventSetupRecord& re3D = iSetup.get<BTagTrackProbability3DRcd>();
0469 unsigned long long cacheId2D = re2D.cacheIdentifier();
0470 unsigned long long cacheId3D = re3D.cacheIdentifier();
0471 if (cacheId2D != calibrationCacheId2D_ || cacheId3D != calibrationCacheId3D_)
0472 {
0473 ESHandle<TrackProbabilityCalibration> calib2DHandle = iSetup.getHandle(calib2d_token_);
0474 ESHandle<TrackProbabilityCalibration> calib3DHandle = iSetup.getHandle(calib3d_token_);
0475 probabilityEstimator_ =
0476 std::make_unique<HistogramProbabilityEstimator>(calib3DHandle.product(), calib2DHandle.product());
0477 }
0478
0479 calibrationCacheId3D_ = cacheId3D;
0480 calibrationCacheId2D_ = cacheId2D;
0481 }
0482
0483
0484 DEFINE_FWK_MODULE(DeepFlavourTagInfoProducer);