diff --git a/RecoTauTag/RecoTau/interface/DeepTauBase.h b/RecoTauTag/RecoTau/interface/DeepTauBase.h index b140baad733b4..259bd1baaf123 100644 --- a/RecoTauTag/RecoTau/interface/DeepTauBase.h +++ b/RecoTauTag/RecoTau/interface/DeepTauBase.h @@ -96,6 +96,8 @@ class DeepTauBase : public edm::stream::EDProducer tausToken_; + edm::EDGetTokenT pfcand_token_; + edm::EDGetTokenT vtx_token_; std::map workingPoints_; OutputCollection outputs_; const DeepTauCache* cache_; diff --git a/RecoTauTag/RecoTau/plugins/DPFIsolation.cc b/RecoTauTag/RecoTau/plugins/DPFIsolation.cc index 541e1954f8509..ef8adce50a59d 100644 --- a/RecoTauTag/RecoTau/plugins/DPFIsolation.cc +++ b/RecoTauTag/RecoTau/plugins/DPFIsolation.cc @@ -67,18 +67,15 @@ class DPFIsolation : public deep_tau::DeepTauBase { explicit DPFIsolation(const edm::ParameterSet& cfg,const deep_tau::DeepTauCache* cache) : DeepTauBase(cfg, GetOutputs(), cache), - pfcand_token(consumes(cfg.getParameter("pfcands"))), - vtx_token(consumes(cfg.getParameter("vertices"))), graphVersion(cfg.getParameter("version")) { const auto& shape = cache_->getGraph().node(0).attr().at("shape").shape(); if(!(graphVersion == 1 || graphVersion == 0 )) - throw cms::Exception("DPFIsolation") << "unknown version of the graph_ file."; + throw cms::Exception("DPFIsolation") << "unknown version of the graph file."; if(!(shape.dim(1).size() == getNumberOfParticles(graphVersion) && shape.dim(2).size() == GetNumberOfFeatures(graphVersion))) throw cms::Exception("DPFIsolation") << "number of inputs does not match the expected inputs for the given version"; - } private: @@ -86,10 +83,10 @@ class DPFIsolation : public deep_tau::DeepTauBase { edm::Handle taus) override { edm::Handle pfcands; - event.getByToken(pfcand_token, pfcands); + event.getByToken(pfcand_token_, pfcands); edm::Handle vertices; - event.getByToken(vtx_token, vertices); + event.getByToken(vtx_token_, vertices); tensorflow::Tensor tensor(tensorflow::DT_FLOAT, {1, static_cast(getNumberOfParticles(graphVersion)), static_cast(GetNumberOfFeatures(graphVersion))}); @@ -394,8 +391,6 @@ class DPFIsolation : public deep_tau::DeepTauBase { } private: - edm::EDGetTokenT pfcand_token; - edm::EDGetTokenT vtx_token; unsigned graphVersion; }; diff --git a/RecoTauTag/RecoTau/plugins/DeepTauId.cc b/RecoTauTag/RecoTau/plugins/DeepTauId.cc index e8f4fd8fc5b0d..f1902e1c7563b 100644 --- a/RecoTauTag/RecoTau/plugins/DeepTauId.cc +++ b/RecoTauTag/RecoTau/plugins/DeepTauId.cc @@ -82,7 +82,7 @@ namespace MuonSubdetId { enum { DT = 1, CSC = 2, RPC = 3, GEM = 4, ME0 = 5 }; } -struct MuonHitMatch { +struct MuonHitMatchV1 { static constexpr int n_muon_stations = 4; std::map> n_matches, n_hits; @@ -90,7 +90,7 @@ struct MuonHitMatch { const pat::Muon* best_matched_muon{nullptr}; double deltaR2_best_match{-1}; - MuonHitMatch() + MuonHitMatchV1() { n_matches[MuonSubdetId::DT].assign(n_muon_stations, 0); n_matches[MuonSubdetId::CSC].assign(n_muon_stations, 0); @@ -107,8 +107,8 @@ struct MuonHitMatch { ++n_muons; const double dR2 = reco::deltaR2(tau.p4(), muon.p4()); if(!best_matched_muon || dR2 < deltaR2_best_match) { - best_matched_muon = &muon; - deltaR2_best_match = dR2; + best_matched_muon = &muon; + deltaR2_best_match = dR2; } for(const auto& segment : muon.matches()) { @@ -225,6 +225,99 @@ struct MuonHitMatch { } }; + +enum class CellObjectType { PfCand_electron, PfCand_muon, PfCand_chargedHadron, PfCand_neutralHadron, + PfCand_gamma, Electron, Muon }; + +template +CellObjectType GetCellObjectType(const Object&); +template<> +CellObjectType GetCellObjectType(const pat::Electron&) { return CellObjectType::Electron; } +template<> +CellObjectType GetCellObjectType(const pat::Muon&) { return CellObjectType::Muon; } + +template<> +CellObjectType GetCellObjectType(const pat::PackedCandidate& cand) +{ + static const std::map obj_types = { + { 11, CellObjectType::PfCand_electron }, + { 13, CellObjectType::PfCand_muon }, + { 22, CellObjectType::PfCand_gamma }, + { 130, CellObjectType::PfCand_neutralHadron }, + { 211, CellObjectType::PfCand_chargedHadron } + }; + + auto iter = obj_types.find(std::abs(cand.pdgId())); + if(iter == obj_types.end()) + throw cms::Exception("DeepTauId") << "Unknown object pdg id = " << cand.pdgId(); + return iter->second; +} + + +using Cell = std::map; +struct CellIndex { + int eta, phi; + + bool operator<(const CellIndex& other) const + { + if(eta != other.eta) return eta < other.eta; + return phi < other.phi; + } +}; + +class CellGrid { +public: + CellGrid(unsigned _nCellsEta, unsigned _nCellsPhi, double _cellSizeEta, double _cellSizePhi) : + nCellsEta(_nCellsEta), nCellsPhi(_nCellsPhi), nTotal(nCellsEta * nCellsPhi), + cellSizeEta(_cellSizeEta), cellSizePhi(_cellSizePhi), cells(nTotal) + { + if(nCellsEta % 2 != 1 || nCellsEta < 1) + throw cms::Exception("DeepTauId") << "Invalid number of eta cells."; + if(nCellsPhi % 2 != 1 || nCellsPhi < 1) + throw cms::Exception("DeepTauId") << "Invalid number of phi cells."; + if(cellSizeEta <= 0 || cellSizePhi <= 0) + throw cms::Exception("DeepTauId") << "Invalid cell size."; + } + + int MaxEtaIndex() const { return static_cast((nCellsEta - 1) / 2); } + int MaxPhiIndex() const { return static_cast((nCellsPhi - 1) / 2); } + double MaxDeltaEta() const { return cellSizeEta * (0.5 + MaxEtaIndex()); } + double MaxDeltaPhi() const { return cellSizePhi * (0.5 + MaxPhiIndex()); } + + bool TryGetCellIndex(double deltaEta, double deltaPhi, CellIndex& cellIndex) const + { + static auto getCellIndex = [](double x, double maxX, double size, int& index) { + const double absX = std::abs(x); + if(absX > maxX) return false; + const double absIndex = std::floor(std::abs(absX / size - 0.5)); + index = static_cast(std::copysign(absIndex, x)); + return true; + }; + + return getCellIndex(deltaEta, MaxDeltaEta(), cellSizeEta, cellIndex.eta) + && getCellIndex(deltaPhi, MaxDeltaPhi(), cellSizePhi, cellIndex.phi); + } + + Cell& at(const CellIndex& cellIndex) { return cells.at(GetFlatIndex(cellIndex)); } + const Cell& at(const CellIndex& cellIndex) const { return cells.at(GetFlatIndex(cellIndex)); } + bool IsEmpty(const CellIndex& cellIndex) const { return at(cellIndex).empty(); } + +private: + size_t GetFlatIndex(const CellIndex& cellIndex) const + { + if(std::abs(cellIndex.eta) > MaxEtaIndex() || std::abs(cellIndex.phi) > MaxPhiIndex()) + throw cms::Exception("DeepTauId") << "Cell index is out of range"; + const unsigned shiftedEta = static_cast(cellIndex.eta + MaxEtaIndex()); + const unsigned shiftedPhi = static_cast(cellIndex.phi + MaxPhiIndex()); + return shiftedEta * nCellsPhi + shiftedPhi; + } + +private: + const unsigned nCellsEta, nCellsPhi, nTotal; + const double cellSizeEta, cellSizePhi; + std::vector cells; +}; + } // anonymous namespace class DeepTauId : public deep_tau::DeepTauBase { @@ -249,8 +342,12 @@ class DeepTauId : public deep_tau::DeepTauBase { desc.add("electrons", edm::InputTag("slimmedElectrons")); desc.add("muons", edm::InputTag("slimmedMuons")); desc.add("taus", edm::InputTag("slimmedTaus")); - desc.add("graph_file", "RecoTauTag/TrainingFiles/data/DeepTauId/deepTau_2017v1_20L1024N_quantized.pb"); + desc.add("pfcands", edm::InputTag("packedPFCandidates")); + desc.add("vertices", edm::InputTag("offlineSlimmedPrimaryVertices")); + desc.add("rho", edm::InputTag("fixedGridRhoAll")); + desc.add("graph_file", "RecoTauTag/TrainingFiles/data/DeepTauId/deepTau_2017v2.pb"); desc.add("mem_mapped", false); + desc.add("version", 2); edm::ParameterSetDescription descWP; descWP.add("VVVLoose", "0"); @@ -265,21 +362,28 @@ class DeepTauId : public deep_tau::DeepTauBase { desc.add("VSeWP", descWP); desc.add("VSmuWP", descWP); desc.add("VSjetWP", descWP); - descriptions.add("DeepTau2017v1", desc); + descriptions.add("DeepTau", desc); } public: explicit DeepTauId(const edm::ParameterSet& cfg, const deep_tau::DeepTauCache* cache) : DeepTauBase(cfg, GetOutputs(), cache), - electrons_token(consumes(cfg.getParameter("electrons"))), - muons_token(consumes(cfg.getParameter("muons"))), - input_layer(cache_->getGraph().node(0).name()), - output_layer(cache_->getGraph().node(cache_->getGraph().node_size() - 1).name()) + electrons_token_(consumes(cfg.getParameter("electrons"))), + muons_token_(consumes(cfg.getParameter("muons"))), + rho_token_(consumes(cfg.getParameter("rho"))), + version(cfg.getParameter("version")) { - const auto& shape = cache_->getGraph().node(0).attr().at("shape").shape(); - if(shape.dim(1).size() != dnn_inputs_2017v1::NumberOfInputs) - throw cms::Exception("DeepTauId") << "number of inputs does not match the expected inputs for the given version"; + if(version == 1) { + input_layer_ = cache_->getGraph().node(0).name(); + output_layer_ = cache_->getGraph().node(cache_->getGraph().node_size() - 1).name(); + const auto& shape = cache_->getGraph().node(0).attr().at("shape").shape(); + if(shape.dim(1).size() != dnn_inputs_2017v1::NumberOfInputs) + throw cms::Exception("DeepTauId") << "number of inputs does not match the expected inputs for the given version"; + } else if(version == 2) { + } else { + throw cms::Exception("DeepTauId") << "version " << version << " is not supported."; + } } static std::unique_ptr initializeGlobalCache(const edm::ParameterSet& cfg) @@ -294,29 +398,140 @@ class DeepTauId : public deep_tau::DeepTauBase { private: tensorflow::Tensor getPredictions(edm::Event& event, const edm::EventSetup& es, - edm::Handle taus) override + edm::Handle taus) override { edm::Handle electrons; - event.getByToken(electrons_token, electrons); + event.getByToken(electrons_token_, electrons); edm::Handle muons; - event.getByToken(muons_token, muons); + event.getByToken(muons_token_, muons); + + edm::Handle pfCands; + event.getByToken(pfcand_token_, pfCands); + + edm::Handle vertices; + event.getByToken(vtx_token_, vertices); + + edm::Handle rho; + event.getByToken(rho_token_, rho); tensorflow::Tensor predictions(tensorflow::DT_FLOAT, { static_cast(taus->size()), dnn_inputs_2017v1::NumberOfOutputs}); for(size_t tau_index = 0; tau_index < taus->size(); ++tau_index) { - const tensorflow::Tensor& inputs = createInputs(taus->at(tau_index), *electrons, *muons); std::vector pred_vector; - tensorflow::run(&(cache_->getSession()), { { input_layer, inputs } }, { output_layer }, &pred_vector); + if(version == 1) + getPredictionsV1(taus->at(tau_index), *electrons, *muons, pred_vector); + else if(version == 2) + getPredictionsV2(taus->at(tau_index), *electrons, *muons, *pfCands, vertices->at(0), *rho, pred_vector); + else + throw cms::Exception("DeepTauId") << "version " << version << " is not supported."; for(int k = 0; k < dnn_inputs_2017v1::NumberOfOutputs; ++k) predictions.matrix()(tau_index, k) = pred_vector[0].flat()(k); } return predictions; } + void getPredictionsV1(const TauType& tau, const pat::ElectronCollection& electrons, + const pat::MuonCollection& muons, std::vector& pred_vector) + { + const tensorflow::Tensor& inputs = createInputsV1(tau, electrons, muons); + tensorflow::run(&(cache_->getSession()), { { input_layer_, inputs } }, { output_layer_ }, &pred_vector); + } + + void getPredictionsV2(const TauType& tau, const pat::ElectronCollection& electrons, + const pat::MuonCollection& muons, const pat::PackedCandidateCollection& pfCands, + const reco::Vertex& pv, double rho, std::vector& pred_vector) + { + CellGrid inner_grid(11, 11, 0.02, 0.02); + CellGrid outer_grid(21, 21, 0.05, 0.05); + fillGrids(tau, electrons, inner_grid, outer_grid); + fillGrids(tau, muons, inner_grid, outer_grid); + fillGrids(tau, pfCands, inner_grid, outer_grid); + + const auto input_tau = createTauBlockInputs(tau, pv, rho); + const auto input_inner_egamma = createEgammaBlockInputs(tau, pv, rho, electrons, pfCands, inner_grid, true); + const auto input_inner_muon = createMuonBlockInputs(tau, pv, rho, muons, pfCands, inner_grid, true); + const auto input_inner_hadrons = createHadronsBlockInputs(tau, pv, rho, pfCands, inner_grid, true); + const auto input_outer_egamma = createEgammaBlockInputs(tau, pv, rho, electrons, pfCands, outer_grid, false); + const auto input_outer_muon = createMuonBlockInputs(tau, pv, rho, muons, pfCands, outer_grid, false); + const auto input_outer_hadrons = createHadronsBlockInputs(tau, pv, rho, pfCands, outer_grid, false); + + tensorflow::run(&(cache_->getSession()), + { { "input_tau", input_tau }, + { "input_inner_egamma", input_inner_egamma}, { "input_outer_egamma", input_outer_egamma }, + { "input_inner_muon", input_inner_muon }, { "input_outer_muon", input_outer_muon }, + { "input_inner_hadrons", input_inner_hadrons }, { "input_outer_hadrons", input_outer_hadrons } }, + { "main_output" }, &pred_vector); + } + + template + void fillGrids(const TauType& tau, const Collection& objects, CellGrid& inner_grid, CellGrid& outer_grid) + { + static constexpr double outer_dR2 = std::pow(0.5, 2); + const double inner_radius = getInnerSignalConeRadius(tau.polarP4().pt()); + const double inner_dR2 = std::pow(inner_radius, 2); + + const auto addObject = [&](size_t n, double deta, double dphi, CellGrid& grid) { + const auto& obj = objects.at(n); + const CellObjectType obj_type = GetCellObjectType(obj); + CellIndex cell_index; + if(grid.TryGetCellIndex(deta, dphi, cell_index)) { + Cell& cell = grid.at(cell_index); + auto iter = cell.find(obj_type); + if(iter != cell.end()) { + const auto& prev_obj = objects.at(iter->second); + if(obj.polarP4().pt() > prev_obj.polarP4().pt()) + iter->second = n; + } else { + cell[obj_type] = n; + } + } + }; + + for(size_t n = 0; n < objects.size(); ++n) { + const auto& obj = objects.at(n); + const double deta = obj.polarP4().eta() - tau.polarP4().eta(); + const double dphi = ROOT::Math::VectorUtil::DeltaPhi(obj.polarP4(), tau.polarP4()); + const double dR2 = std::pow(deta, 2) + std::pow(dphi, 2); + if(dR2 < inner_dR2) + addObject(n, deta, dphi, inner_grid); + if(dR2 < outer_dR2) + addObject(n, deta, dphi, outer_grid); + } + } + + tensorflow::Tensor createTauBlockInputs(const TauType& tau, const reco::Vertex& pv, double rho) + { + return tensorflow::Tensor(); + } + + tensorflow::Tensor createEgammaBlockInputs(const TauType& tau, const reco::Vertex& pv, double rho, + const pat::ElectronCollection& electrons, + const pat::PackedCandidateCollection& pfCands, + const CellGrid& grid, bool is_inner) + { + return tensorflow::Tensor(); + } + + tensorflow::Tensor createMuonBlockInputs(const TauType& tau, const reco::Vertex& pv, double rho, + const pat::MuonCollection& electrons, + const pat::PackedCandidateCollection& pfCands, + const CellGrid& grid, bool is_inner) + { + return tensorflow::Tensor(); + } + + tensorflow::Tensor createHadronsBlockInputs(const TauType& tau, const reco::Vertex& pv, double rho, + const pat::PackedCandidateCollection& pfCands, + const CellGrid& grid, bool is_inner) + { + return tensorflow::Tensor(); + } + + template - tensorflow::Tensor createInputs(const TauType& tau, const ElectronCollection& electrons, - const MuonCollection& muons) const + tensorflow::Tensor createInputsV1(const TauType& tau, const ElectronCollection& electrons, + const MuonCollection& muons) const { static constexpr bool check_all_set = false; static constexpr float default_value_for_set_check = -42; @@ -421,7 +636,7 @@ class DeepTauId : public deep_tau::DeepTauBase { / tau.leadChargedHadrCand()->pt(); } - MuonHitMatch muon_hit_match; + MuonHitMatchV1 muon_hit_match; if(tau.leadPFChargedHadrCand().isNonnull() && tau.leadPFChargedHadrCand()->muonRef().isNonnull()) muon_hit_match.addMatchedMuon(*tau.leadPFChargedHadrCand()->muonRef(), tau); @@ -628,9 +843,11 @@ class DeepTauId : public deep_tau::DeepTauBase { } private: - edm::EDGetTokenT electrons_token; - edm::EDGetTokenT muons_token; - std::string input_layer, output_layer; + edm::EDGetTokenT electrons_token_; + edm::EDGetTokenT muons_token_; + edm::EDGetTokenT rho_token_; + std::string input_layer_, output_layer_; + const unsigned version; }; #include "FWCore/Framework/interface/MakerMacros.h" diff --git a/RecoTauTag/RecoTau/python/tools/runTauIdMVA.py b/RecoTauTag/RecoTau/python/tools/runTauIdMVA.py index 3b7f1a62d2df8..def331fad2e10 100644 --- a/RecoTauTag/RecoTau/python/tools/runTauIdMVA.py +++ b/RecoTauTag/RecoTau/python/tools/runTauIdMVA.py @@ -633,8 +633,12 @@ def runTauID(self): electrons = self.cms.InputTag('slimmedElectrons'), muons = self.cms.InputTag('slimmedMuons'), taus = self.cms.InputTag('slimmedTaus'), + pfcands = self.cms.InputTag('packedPFCandidates'), + vertices = self.cms.InputTag('offlineSlimmedPrimaryVertices'), + rho = self.cms.InputTag('fixedGridRhoAll'), graph_file = self.cms.string(file_name), - mem_mapped = self.cms.bool(False) + mem_mapped = self.cms.bool(False), + version = self.cms.uint32(1) ) self.processDeepProducer('deepTau2017v1', tauIDSources, workingPoints_) diff --git a/RecoTauTag/RecoTau/src/DeepTauBase.cc b/RecoTauTag/RecoTau/src/DeepTauBase.cc index 459ce0ebbfc08..0918bfbf206e9 100644 --- a/RecoTauTag/RecoTau/src/DeepTauBase.cc +++ b/RecoTauTag/RecoTau/src/DeepTauBase.cc @@ -76,8 +76,11 @@ DeepTauBase::Output::ResultMap DeepTauBase::Output::get_value(const edm::Handle< return output; } -DeepTauBase::DeepTauBase(const edm::ParameterSet& cfg, const OutputCollection& outputCollection, const DeepTauCache* cache) : +DeepTauBase::DeepTauBase(const edm::ParameterSet& cfg, const OutputCollection& outputCollection, + const DeepTauCache* cache) : tausToken_(consumes(cfg.getParameter("taus"))), + pfcand_token_(consumes(cfg.getParameter("pfcands"))), + vtx_token_(consumes(cfg.getParameter("vertices"))), outputs_(outputCollection), cache_(cache) {