From 016218c356fe44095e4ef7332089034d941d0d87 Mon Sep 17 00:00:00 2001 From: alintulu Date: Tue, 8 Nov 2022 16:39:24 +0100 Subject: [PATCH] Add custom Run 3 scouting NanoAOD --- .../plugins/JetDeltaRValueMapProducer.cc | 37 +- .../NanoAOD/plugins/NanoAODOutputModule.cc | 7 +- ...outingParticleToRecoPFCandidateProducer.cc | 271 ++++++ .../plugins/SimpleFlatTableProducerPlugins.cc | 21 + .../NanoAOD/plugins/TriggerOutputBranches.cc | 6 +- .../NanoAOD/plugins/TriggerOutputBranches.h | 5 +- .../NanoAOD/python/custom_run3scouting_cff.py | 582 +++++++++++++ .../plugins/DeepBoostedJetTagInfoProducer.cc | 800 +++++++++++------- .../plugins/BoostedJetONNXValueMapProducer.cc | 237 ++++++ 9 files changed, 1632 insertions(+), 334 deletions(-) create mode 100644 PhysicsTools/NanoAOD/plugins/Run3ScoutingParticleToRecoPFCandidateProducer.cc create mode 100644 PhysicsTools/NanoAOD/python/custom_run3scouting_cff.py create mode 100644 RecoBTag/ONNXRuntime/plugins/BoostedJetONNXValueMapProducer.cc diff --git a/CommonTools/RecoAlgos/plugins/JetDeltaRValueMapProducer.cc b/CommonTools/RecoAlgos/plugins/JetDeltaRValueMapProducer.cc index f31f88b9a97ce..84ad01353bd7a 100644 --- a/CommonTools/RecoAlgos/plugins/JetDeltaRValueMapProducer.cc +++ b/CommonTools/RecoAlgos/plugins/JetDeltaRValueMapProducer.cc @@ -27,10 +27,12 @@ #include "DataFormats/Common/interface/ValueMap.h" #include "DataFormats/Math/interface/deltaR.h" -template +template class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> { public: - typedef edm::ValueMap JetValueMap; + typedef TN value_type; + typedef edm::ValueMap JetValueMap; + typedef std::vector ValueVector; JetDeltaRValueMapProducer(edm::ParameterSet const& params) : srcToken_(consumes >(params.getParameter("src"))), @@ -46,8 +48,10 @@ class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> { lazyParser_(params.existsAs("lazyParser") ? params.getParameter("lazyParser") : false), multiValue_(false) { if (!value_.empty()) { - evaluationMap_.insert(std::make_pair( - value_, std::unique_ptr >(new StringObjectFunction(value_, lazyParser_)))); + if (value_ != "index") { + evaluationMap_.insert(std::make_pair( + value_, std::unique_ptr >(new StringObjectFunction(value_, lazyParser_)))); + } produces(); } @@ -75,7 +79,8 @@ class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> { edm::Handle > h_jets2; iEvent.getByToken(matchedToken_, h_jets2); - std::vector values(h_jets1->size(), -99999); + ValueVector values(h_jets1->size(), -99999); + std::map > valuesMap; if (multiValue_) { for (size_t i = 0; i < valueLabels_.size(); ++i) @@ -88,19 +93,20 @@ class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> { ++ijet) { float matched_dR2 = 1e9; int matched_index = -1; + int iindex = ijet - ibegin; for (typename edm::View::const_iterator jbegin = h_jets1->begin(), jend = h_jets1->end(), jjet = jbegin; jjet != jend; ++jjet) { - int index = jjet - jbegin; + int jindex = jjet - jbegin; - if (jets1_locks.at(index)) + if (jets1_locks.at(jindex)) continue; // skip jets that have already been matched float temp_dR2 = reco::deltaR2(ijet->eta(), ijet->phi(), jjet->eta(), jjet->phi()); if (temp_dR2 < matched_dR2) { matched_dR2 = temp_dR2; - matched_index = index; + matched_index = jindex; } } // end loop over src jets @@ -109,8 +115,13 @@ class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> { edm::LogInfo("MatchedJetsFarApart") << "Matched jets separated by dR greater than distMax=" << distMax_; else { jets1_locks.at(matched_index) = true; - if (!value_.empty()) - values.at(matched_index) = (*(evaluationMap_.at(value_)))(*ijet); + if (!value_.empty()) { + if (value_ == "index") { + values.at(matched_index) = iindex; + } else { + values.at(matched_index) = (*(evaluationMap_.at(value_)))(*ijet); + } + } if (multiValue_) { for (size_t i = 0; i < valueLabels_.size(); ++i) valuesMap.at(valueLabels_[i]).at(matched_index) = (*(evaluationMap_.at(valueLabels_[i])))(*ijet); @@ -122,7 +133,7 @@ class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> { if (!value_.empty()) { std::unique_ptr jetValueMap(new JetValueMap()); - JetValueMap::Filler filler(*jetValueMap); + typename JetValueMap::Filler filler(*jetValueMap); filler.insert(h_jets1, values.begin(), values.end()); filler.fill(); @@ -133,7 +144,7 @@ class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> { for (size_t i = 0; i < valueLabels_.size(); ++i) { std::unique_ptr jetValueMap(new JetValueMap()); - JetValueMap::Filler filler(*jetValueMap); + typename JetValueMap::Filler filler(*jetValueMap); filler.insert(h_jets1, valuesMap.at(valueLabels_[i]).begin(), valuesMap.at(valueLabels_[i]).end()); filler.fill(); @@ -157,7 +168,9 @@ class JetDeltaRValueMapProducer : public edm::stream::EDProducer<> { typedef JetDeltaRValueMapProducer RecoJetDeltaRValueMapProducer; typedef JetDeltaRValueMapProducer RecoJetToPatJetDeltaRValueMapProducer; typedef JetDeltaRValueMapProducer PatJetDeltaRValueMapProducer; +typedef JetDeltaRValueMapProducer RecoJetToGenJetDeltaRValueMapProducer; DEFINE_FWK_MODULE(RecoJetDeltaRValueMapProducer); DEFINE_FWK_MODULE(RecoJetToPatJetDeltaRValueMapProducer); DEFINE_FWK_MODULE(PatJetDeltaRValueMapProducer); +DEFINE_FWK_MODULE(RecoJetToGenJetDeltaRValueMapProducer); diff --git a/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc b/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc index 7c5ac3cb70cf2..3cb8a0c33d796 100644 --- a/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc +++ b/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc @@ -77,6 +77,7 @@ class NanoAODOutputModule : public edm::one::OutputModule<> { std::unique_ptr m_tree, m_lumiTree, m_runTree, m_metaDataTree, m_parameterSetsTree; static constexpr int m_firstFlush{1000}; + bool m_isPFScouting; class CommonEventBranches { public: @@ -159,7 +160,8 @@ NanoAODOutputModule::NanoAODOutputModule(edm::ParameterSet const& pset) m_writeProvenance(pset.getUntrackedParameter("saveProvenance", true)), m_fakeName(pset.getUntrackedParameter("fakeNameForCrab", false)), m_autoFlush(pset.getUntrackedParameter("autoFlush", -10000000)), - m_processHistoryRegistry() {} + m_processHistoryRegistry(), + m_isPFScouting(pset.getUntrackedParameter("isPFScouting", false)) {} NanoAODOutputModule::~NanoAODOutputModule() {} @@ -322,7 +324,7 @@ void NanoAODOutputModule::openFile(edm::FileBlock const&) { if (keep.first->className() == "nanoaod::FlatTable") m_tables.emplace_back(keep.first, keep.second); else if (keep.first->className() == "edm::TriggerResults") { - m_triggers.emplace_back(keep.first, keep.second); + m_triggers.emplace_back(keep.first, keep.second, m_isPFScouting); } else if (keep.first->className() == "std::basic_string >" && keep.first->productInstanceName() == "genModel") { // friendlyClassName == "String" m_evstrings.emplace_back(keep.first, keep.second, true); // update only at lumiBlock transitions @@ -416,6 +418,7 @@ void NanoAODOutputModule::fillDescriptions(edm::ConfigurationDescriptions& descr "Change the OutputModule name in the fwk job report to fake PoolOutputModule. This is needed to run on cran " "(and publish) till crab is fixed"); desc.addUntracked("autoFlush", -10000000)->setComment("Autoflush parameter for ROOT file"); + desc.addUntracked("isPFScouting", false); //replace with whatever you want to get from the EDM by default const std::vector keep = {"drop *", diff --git a/PhysicsTools/NanoAOD/plugins/Run3ScoutingParticleToRecoPFCandidateProducer.cc b/PhysicsTools/NanoAOD/plugins/Run3ScoutingParticleToRecoPFCandidateProducer.cc new file mode 100644 index 0000000000000..73cf2464b2d66 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/Run3ScoutingParticleToRecoPFCandidateProducer.cc @@ -0,0 +1,271 @@ +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h" +#include "DataFormats/Scouting/interface/Run3ScoutingParticle.h" +#include "DataFormats/Common/interface/OrphanHandle.h" + +#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h" +#include "fastjet/contrib/SoftKiller.hh" + +class Run3ScoutingParticleToRecoPFCandidateProducer : public edm::stream::EDProducer<> { +public: + explicit Run3ScoutingParticleToRecoPFCandidateProducer(const edm::ParameterSet &); + ~Run3ScoutingParticleToRecoPFCandidateProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + void beginStream(edm::StreamID) override {} + void produce(edm::Event & iEvent, edm::EventSetup const & setup) override; + void endStream() override {} + + void createPFCandidates(edm::Handle> scoutingparticleHandle, std::unique_ptr & pfcands); + void createPFCandidatesSK(edm::Handle> scoutingparticleHandle, std::unique_ptr & pfcands); + reco::PFCandidate createPFCand(Run3ScoutingParticle scoutingparticle); + void clearVars(); + +private: + const edm::EDGetTokenT> input_scoutingparticle_token_; + const edm::ESGetToken particletable_token_; + bool use_softKiller_; + bool use_CHS_; + const HepPDT::ParticleDataTable* pdTable_; + + std::vector vertexIndex_; + std::vector normchi2_; + std::vector dz_; + std::vector dxy_; + std::vector dzsig_; + std::vector dxysig_; + std::vector lostInnerHits_; + std::vector quality_; + std::vector trkPt_; + std::vector trkEta_; + std::vector trkPhi_; +}; + +// +// constructors and destructor +// +Run3ScoutingParticleToRecoPFCandidateProducer::Run3ScoutingParticleToRecoPFCandidateProducer(const edm::ParameterSet &iConfig) + : input_scoutingparticle_token_(consumes(iConfig.getParameter("scoutingparticle"))), + particletable_token_(esConsumes()), + use_softKiller_(iConfig.getParameter("softKiller")), + use_CHS_(iConfig.getParameter("CHS")) { + + //register products + produces(); + produces>("vertexIndex"); + produces>("normchi2"); + produces>("dz"); + produces>("dxy"); + produces>("dzsig"); + produces>("dxysig"); + produces>("lostInnerHits"); + produces>("quality"); + produces>("trkPt"); + produces>("trkEta"); + produces>("trkPhi"); +} + +Run3ScoutingParticleToRecoPFCandidateProducer::~Run3ScoutingParticleToRecoPFCandidateProducer() = default; + +reco::PFCandidate Run3ScoutingParticleToRecoPFCandidateProducer::createPFCand(Run3ScoutingParticle scoutingparticle) { + auto m = pdTable_->particle(HepPDT::ParticleID(scoutingparticle.pdgId())) != nullptr + ? pdTable_->particle(HepPDT::ParticleID(scoutingparticle.pdgId()))->mass() + : -99.f; + auto q = pdTable_->particle(HepPDT::ParticleID(scoutingparticle.pdgId())) != nullptr + ? pdTable_->particle(HepPDT::ParticleID(scoutingparticle.pdgId()))->charge() + : -99.f; + if (m < -90 or q < -90) { + LogDebug("createPFCand") << ":" << std::endl + << "Unrecognisable pdgId - skipping particle" << std::endl; + return reco::PFCandidate(); + } + + float px = scoutingparticle.pt() * cos(scoutingparticle.phi()); + float py = scoutingparticle.pt() * sin(scoutingparticle.phi()); + float pz = scoutingparticle.pt() * sinh(scoutingparticle.eta()); + float p = scoutingparticle.pt() * cosh(scoutingparticle.eta()); + float energy = std::sqrt(p*p + m*m); + reco::Particle::LorentzVector p4(px, py, pz, energy); + + static const reco::PFCandidate dummy; + auto pfcand = reco::PFCandidate(q, p4, dummy.translatePdgIdToType(scoutingparticle.pdgId())); + + bool relativeTrackVars = scoutingparticle.relative_trk_vars(); + vertexIndex_.push_back(scoutingparticle.vertex()); + normchi2_.push_back(scoutingparticle.normchi2()); + dz_.push_back(scoutingparticle.dz()); + dxy_.push_back(scoutingparticle.dxy()); + dzsig_.push_back(scoutingparticle.dzsig()); + dxysig_.push_back(scoutingparticle.dxysig()); + lostInnerHits_.push_back(scoutingparticle.lostInnerHits()); + quality_.push_back(scoutingparticle.quality()); + trkPt_.push_back(relativeTrackVars ? scoutingparticle.trk_pt() + scoutingparticle.pt() : scoutingparticle.trk_pt()); + trkEta_.push_back(relativeTrackVars ? scoutingparticle.trk_eta() + scoutingparticle.eta() : scoutingparticle.trk_eta()); + trkPhi_.push_back(relativeTrackVars ? scoutingparticle.trk_phi() + scoutingparticle.phi() : scoutingparticle.trk_phi()); + + return pfcand; +} + +void Run3ScoutingParticleToRecoPFCandidateProducer::createPFCandidates(edm::Handle> scoutingparticleHandle, std::unique_ptr & pfcands) { + for (unsigned int icand = 0; icand < scoutingparticleHandle->size(); ++icand) { + + auto& scoutingparticle = (*scoutingparticleHandle)[icand]; + + if (use_CHS_ and scoutingparticle.vertex() != 0) continue; + + auto pfcand = createPFCand(scoutingparticle); + if (pfcand.energy() != 0) pfcands->push_back(pfcand); + } +} + +void Run3ScoutingParticleToRecoPFCandidateProducer::createPFCandidatesSK(edm::Handle> scoutingparticleHandle, std::unique_ptr & pfcands) { + std::vector fj; + + for (auto iter = scoutingparticleHandle->begin(), ibegin = scoutingparticleHandle->begin(), iend = scoutingparticleHandle->end(); iter != iend; ++iter) { + auto m = pdTable_->particle(HepPDT::ParticleID(iter->pdgId())) != nullptr + ? pdTable_->particle(HepPDT::ParticleID(iter->pdgId()))->mass() + : -99.f; + if (m < -90) { + LogDebug("createPFCandidatesSK") << ":" << std::endl + << "Unrecognisable pdgId - skipping particle" << std::endl; + continue; + } + math::PtEtaPhiMLorentzVector p4(iter->pt(), iter->eta(), iter->phi(), m); + fj.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.energy())); + fj.back().set_user_index(iter - ibegin); + } + + fastjet::contrib::SoftKiller soft_killer(5, 0.4); + std::vector soft_killed_particles = soft_killer(fj); + + for (auto &particle : soft_killed_particles) { + const Run3ScoutingParticle scoutingparticle = scoutingparticleHandle->at(particle.user_index()); + auto pfcand = createPFCand(scoutingparticle); + if (pfcand.energy() != 0) pfcands->push_back(pfcand); + } +} + +// ------------ method called to produce the data ------------ +void Run3ScoutingParticleToRecoPFCandidateProducer::produce(edm::Event & iEvent, edm::EventSetup const & setup) { + using namespace edm; + + auto pdt = setup.getHandle(particletable_token_); + pdTable_ = pdt.product(); + + Handle> scoutingparticleHandle; + iEvent.getByToken(input_scoutingparticle_token_, scoutingparticleHandle); + + auto pfcands = std::make_unique(); + + if (use_softKiller_) { + createPFCandidatesSK(scoutingparticleHandle, pfcands); + } else { + createPFCandidates(scoutingparticleHandle, pfcands); + } + + edm::OrphanHandle oh = iEvent.put(std::move(pfcands)); + + std::unique_ptr> vertexIndex_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_vertexIndex(*vertexIndex_VM); + filler_vertexIndex.insert(oh, vertexIndex_.begin(), vertexIndex_.end()); + filler_vertexIndex.fill(); + iEvent.put(std::move(vertexIndex_VM), "vertexIndex"); + + std::unique_ptr> normchi2_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_normchi2(*normchi2_VM); + filler_normchi2.insert(oh, normchi2_.begin(), normchi2_.end()); + filler_normchi2.fill(); + iEvent.put(std::move(normchi2_VM), "normchi2"); + + std::unique_ptr> dz_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_dz(*dz_VM); + filler_dz.insert(oh, dz_.begin(), dz_.end()); + filler_dz.fill(); + iEvent.put(std::move(dz_VM), "dz"); + + std::unique_ptr> dxy_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_dxy(*dxy_VM); + filler_dxy.insert(oh, dxy_.begin(), dxy_.end()); + filler_dxy.fill(); + iEvent.put(std::move(dxy_VM), "dxy"); + + std::unique_ptr> dzsig_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_dzsig(*dzsig_VM); + filler_dzsig.insert(oh, dzsig_.begin(), dzsig_.end()); + filler_dzsig.fill(); + iEvent.put(std::move(dzsig_VM), "dzsig"); + + std::unique_ptr> dxysig_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_dxysig(*dxysig_VM); + filler_dxysig.insert(oh, dxysig_.begin(), dxysig_.end()); + filler_dxysig.fill(); + iEvent.put(std::move(dxysig_VM), "dxysig"); + + std::unique_ptr> lostInnerHits_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_lostInnerHits(*lostInnerHits_VM); + filler_lostInnerHits.insert(oh, lostInnerHits_.begin(), lostInnerHits_.end()); + filler_lostInnerHits.fill(); + iEvent.put(std::move(lostInnerHits_VM), "lostInnerHits"); + + std::unique_ptr> quality_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_quality(*quality_VM); + filler_quality.insert(oh, quality_.begin(), quality_.end()); + filler_quality.fill(); + iEvent.put(std::move(quality_VM), "quality"); + + std::unique_ptr> trkPt_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_trkPt(*trkPt_VM); + filler_trkPt.insert(oh, trkPt_.begin(), trkPt_.end()); + filler_trkPt.fill(); + iEvent.put(std::move(trkPt_VM), "trkPt"); + + std::unique_ptr> trkEta_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_trkEta(*trkEta_VM); + filler_trkEta.insert(oh, trkEta_.begin(), trkEta_.end()); + filler_trkEta.fill(); + iEvent.put(std::move(trkEta_VM), "trkEta"); + + std::unique_ptr> trkPhi_VM(new edm::ValueMap()); + edm::ValueMap::Filler filler_trkPhi(*trkPhi_VM); + filler_trkPhi.insert(oh, trkPhi_.begin(), trkPhi_.end()); + filler_trkPhi.fill(); + iEvent.put(std::move(trkPhi_VM), "trkPhi"); + + clearVars(); +} + +void Run3ScoutingParticleToRecoPFCandidateProducer::clearVars() { + vertexIndex_.clear(); + normchi2_.clear(); + dz_.clear(); + dxy_.clear(); + dzsig_.clear(); + dxysig_.clear(); + lostInnerHits_.clear(); + quality_.clear(); + trkPt_.clear(); + trkEta_.clear(); + trkPhi_.clear(); +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void Run3ScoutingParticleToRecoPFCandidateProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.add("scoutingparticle", edm::InputTag("hltScoutingPFPacker")); + desc.add("softKiller", false); + desc.add("CHS", false); + descriptions.addWithDefaultLabel(desc); +} + +// declare this class as a framework plugin +DEFINE_FWK_MODULE(Run3ScoutingParticleToRecoPFCandidateProducer); diff --git a/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc b/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc index cf30e060c442b..cc2a102e2763f 100644 --- a/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc +++ b/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc @@ -24,6 +24,21 @@ typedef EventSingletonSimpleFlatTableProducer SimpleXYZPointFla #include "DataFormats/BeamSpot/interface/BeamSpot.h" typedef EventSingletonSimpleFlatTableProducer SimpleBeamspotFlatTableProducer; +#include "DataFormats/Scouting/interface/Run3ScoutingVertex.h" +typedef SimpleFlatTableProducer SimpleRun3ScoutingVertexFlatTableProducer; + +#include "DataFormats/Scouting/interface/Run3ScoutingPhoton.h" +typedef SimpleFlatTableProducer SimpleRun3ScoutingPhotonFlatTableProducer; + +#include "DataFormats/Scouting/interface/Run3ScoutingMuon.h" +typedef SimpleFlatTableProducer SimpleRun3ScoutingMuonFlatTableProducer; + +#include "DataFormats/Scouting/interface/Run3ScoutingElectron.h" +typedef SimpleFlatTableProducer SimpleRun3ScoutingElectronFlatTableProducer; + +#include "DataFormats/Scouting/interface/Run3ScoutingTrack.h" +typedef SimpleFlatTableProducer SimpleRun3ScoutingTrackFlatTableProducer; + #include "FWCore/Framework/interface/MakerMacros.h" DEFINE_FWK_MODULE(SimpleCandidateFlatTableProducer); DEFINE_FWK_MODULE(SimpleGenEventFlatTableProducer); @@ -33,3 +48,9 @@ DEFINE_FWK_MODULE(SimpleProtonTrackFlatTableProducer); DEFINE_FWK_MODULE(SimpleLocalTrackFlatTableProducer); DEFINE_FWK_MODULE(SimpleXYZPointFlatTableProducer); DEFINE_FWK_MODULE(SimpleBeamspotFlatTableProducer); +DEFINE_FWK_MODULE(SimpleRun3ScoutingVertexFlatTableProducer); +DEFINE_FWK_MODULE(SimpleRun3ScoutingPhotonFlatTableProducer); +DEFINE_FWK_MODULE(SimpleRun3ScoutingMuonFlatTableProducer); +DEFINE_FWK_MODULE(SimpleRun3ScoutingElectronFlatTableProducer); +DEFINE_FWK_MODULE(SimpleRun3ScoutingTrackFlatTableProducer); + diff --git a/PhysicsTools/NanoAOD/plugins/TriggerOutputBranches.cc b/PhysicsTools/NanoAOD/plugins/TriggerOutputBranches.cc index 047eec11e560d..91cee050f748b 100644 --- a/PhysicsTools/NanoAOD/plugins/TriggerOutputBranches.cc +++ b/PhysicsTools/NanoAOD/plugins/TriggerOutputBranches.cc @@ -21,7 +21,7 @@ void TriggerOutputBranches::updateTriggerNames(TTree& tree, for (unsigned int j = 0; j < newNames.size(); j++) { std::string name = newNames[j]; // no const & as it will be modified below! std::size_t vfound = name.rfind("_v"); - if (vfound != std::string::npos && (name.compare(0, 3, "HLT") == 0 || name.compare(0, 2, "L1") == 0)) { + if (vfound != std::string::npos && ((name.compare(0, 3, "HLT") == 0 && !(m_isPFScouting)) || name.compare(0, 2, "L1") == 0 || (m_isPFScouting && name.find("PFScouting") != std::string::npos))) { name.replace(vfound, name.size() - vfound, ""); } if (name == existing.name) @@ -32,11 +32,11 @@ void TriggerOutputBranches::updateTriggerNames(TTree& tree, for (unsigned int j = 0; j < newNames.size(); j++) { std::string name = newNames[j]; // no const & as it will be modified below! std::size_t vfound = name.rfind("_v"); - if (vfound != std::string::npos && (name.compare(0, 3, "HLT") == 0 || name.compare(0, 2, "L1") == 0)) { + if (vfound != std::string::npos && ((name.compare(0, 3, "HLT") == 0 && !(m_isPFScouting)) || name.compare(0, 2, "L1") == 0 || (m_isPFScouting && name.find("PFScouting") != std::string::npos))) { name.replace(vfound, name.size() - vfound, ""); } bool found = false; - if (name.compare(0, 3, "HLT") == 0 || name.compare(0, 4, "Flag") == 0 || name.compare(0, 2, "L1") == 0) { + if ((name.compare(0, 3, "HLT") == 0 && !(m_isPFScouting)) || name.compare(0, 4, "Flag") == 0 || name.compare(0, 2, "L1") == 0 || (m_isPFScouting && name.find("PFScouting") != std::string::npos)) { for (auto& existing : m_triggerBranches) { if (name == existing.name) found = true; diff --git a/PhysicsTools/NanoAOD/plugins/TriggerOutputBranches.h b/PhysicsTools/NanoAOD/plugins/TriggerOutputBranches.h index 6b2f3fe031285..48f3c9430c961 100644 --- a/PhysicsTools/NanoAOD/plugins/TriggerOutputBranches.h +++ b/PhysicsTools/NanoAOD/plugins/TriggerOutputBranches.h @@ -12,8 +12,8 @@ class TriggerOutputBranches { public: - TriggerOutputBranches(const edm::BranchDescription *desc, const edm::EDGetToken &token) - : m_token(token), m_lastRun(-1), m_fills(0), m_processName(desc->processName()) { + TriggerOutputBranches(const edm::BranchDescription *desc, const edm::EDGetToken &token, bool isPFScouting) + : m_token(token), m_lastRun(-1), m_fills(0), m_processName(desc->processName()), m_isPFScouting(isPFScouting) { if (desc->className() != "edm::TriggerResults") throw cms::Exception("Configuration", "NanoAODOutputModule/TriggerOutputBranches can only write out edm::TriggerResults objects"); @@ -44,6 +44,7 @@ class TriggerOutputBranches { unsigned long m_fills; std::string m_processName; bool verifyBranchUniqueName(TTree &, std::string) const; + bool m_isPFScouting; template void fillColumn(NamedBranchPtr &nb, const edm::TriggerResults &triggers) { diff --git a/PhysicsTools/NanoAOD/python/custom_run3scouting_cff.py b/PhysicsTools/NanoAOD/python/custom_run3scouting_cff.py new file mode 100644 index 0000000000000..f2dad6c84dcbb --- /dev/null +++ b/PhysicsTools/NanoAOD/python/custom_run3scouting_cff.py @@ -0,0 +1,582 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.common_cff import * +from PhysicsTools.NanoAOD.simpleCandidateFlatTableProducer_cfi import simpleCandidateFlatTableProducer + +def AddRun3Scouting(process): + """ + Add all Run 3 scouting objects except jets and particles + """ + process.photonScoutingTable = cms.EDProducer("SimpleRun3ScoutingPhotonFlatTableProducer", + src = cms.InputTag("hltScoutingEgammaPacker"), + cut = cms.string(""), + name = cms.string("ScoutingPhoton"), + doc = cms.string("Photon Scouting Informations"), + singleton = cms.bool(False), # the number of entries is variable + extension = cms.bool(False), # this is the main table for the muons + variables = cms.PSet( + pt = Var('pt', 'float', precision=14, doc='pt'), + eta = Var('eta', 'float', precision=14, doc='eta'), + phi = Var('phi', 'float', precision=14, doc='phi'), + m = Var('m', 'float', precision=14, doc='m'), + sigmaIetaIeta = Var('sigmaIetaIeta', 'float', precision=14, doc='sigmaIetaIeta'), + hOverE = Var('hOverE', 'float', precision=14, doc='hOverE'), + ecalIso = Var('ecalIso', 'float', precision=14, doc='ecalIso'), + hcalIso = Var('hcalIso', 'float', precision=14, doc='hcalIso'), + trackIso = Var('trkIso', 'float', precision=14, doc='trackIso'), + r9 = Var('r9', 'float', precision=14, doc='r9'), + sMin = Var('sMin', 'float', precision=14, doc='sMin'), + sMaj = Var('sMaj', 'float', precision=14, doc='sMaj'), + seedId = Var('seedId', 'int', doc='seedId'), + ) + ) + + process.electronScoutingTable = cms.EDProducer("SimpleRun3ScoutingElectronFlatTableProducer", + src = cms.InputTag("hltScoutingEgammaPacker"), + cut = cms.string(""), + name = cms.string("ScoutingElectron"), + doc = cms.string("Electron Scouting Informations"), + singleton = cms.bool(False), # the number of entries is variable + extension = cms.bool(False), # this is the main table for the muons + variables = cms.PSet( + pt = Var('pt', 'float', precision=14, doc='pt'), + eta = Var('eta', 'float', precision=14, doc='eta'), + phi = Var('phi', 'float', precision=14, doc='phi'), + m = Var('m', 'float', precision=14, doc='m'), + d0 = Var('d0', 'float', precision=14, doc='d0'), + dz = Var('dz', 'float', precision=14, doc='dz'), + dEtaIn = Var('dEtaIn', 'float', precision=14, doc='dEtaIn'), + dPhiIn = Var('dPhiIn', 'float', precision=14, doc='dPhiIn'), + sigmaIetaIeta = Var('sigmaIetaIeta', 'float', precision=14, doc='sigmaIetaIeta'), + hOverE = Var('hOverE', 'float', precision=14, doc='hOverE'), + ooEMOop = Var('ooEMOop', 'float', precision=14, doc='ooEMOop'), + missingHits = Var('missingHits', 'int', doc='missingHits'), + charge = Var('charge', 'int', doc='charge'), + ecalIso = Var('ecalIso', 'float', precision=14, doc='ecalIso'), + hcalIso = Var('hcalIso', 'float', precision=14, doc='hcalIso'), + trackIso = Var('trackIso', 'float', precision=14, doc='trackIso'), + r9 = Var('r9', 'float', precision=14, doc='r9'), + sMin = Var('sMin', 'float', precision=14, doc='sMin'), + sMaj = Var('sMaj', 'float', precision=14, doc='sMaj'), + seedId = Var('seedId', 'int', doc='seedId'), + ) + ) + + process.muonScoutingTable = cms.EDProducer("SimpleRun3ScoutingMuonFlatTableProducer", + src = cms.InputTag("hltScoutingMuonPacker"), + cut = cms.string(""), + name = cms.string("ScoutingMuon"), + doc = cms.string("Muon Scouting Informations"), + singleton = cms.bool(False), # the number of entries is variable + extension = cms.bool(False), # this is the main table for the muons + variables = cms.PSet( + pt = Var('pt', 'float', precision=14, doc='pt'), + eta = Var('eta', 'float', precision=14, doc='eta'), + phi = Var('phi', 'float', precision=14, doc='phi'), + m = Var('m', 'float', precision=14, doc='m'), + type = Var('type', 'int', doc='type'), + charge = Var('charge', 'int', doc='charge'), + normchi2 = Var('normalizedChi2', 'float', precision=14, doc='normalizedChi2'), + ecalIso = Var('ecalIso', 'float', precision=14, doc='ecalIso'), + hcalIso = Var('hcalIso', 'float', precision=14, doc='hcalIso'), + nValidStandAloneMuonHits = Var('nValidStandAloneMuonHits', 'int', doc='nValidStandAloneMuonHits'), + nStandAloneMuonMatchedStations = Var('nStandAloneMuonMatchedStations', 'int', doc='nStandAloneMuonMatchedStations'), + nValidRecoMuonHits = Var('nValidRecoMuonHits', 'int', doc='nValidRecoMuonHits'), + nRecoMuonChambers = Var('nRecoMuonChambers', 'int', doc='nRecoMuonChambers'), + nRecoMuonChambersCSCorDT = Var('nRecoMuonChambersCSCorDT', 'int', doc='nRecoMuonChambersCSCorDT'), + nRecoMuonMatches = Var('nRecoMuonMatches', 'int', doc='nRecoMuonMatches'), + nRecoMuonMatchedStations = Var('nRecoMuonMatchedStations', 'int', doc='nRecoMuonMatchedStations'), + nRecoMuonExpectedMatchedStations = Var('nRecoMuonExpectedMatchedStations', 'int', doc='nRecoMuonExpectedMatchedStations'), + recoMuonStationMask = Var('recoMuonStationMask', 'int', doc='recoMuonStationMask'), + nRecoMuonMatchedRPCLayers = Var('nRecoMuonMatchedRPCLayers', 'int', doc='nRecoMuonMatchedRPCLayers'), + recoMuonRPClayerMask = Var('recoMuonRPClayerMask', 'int', doc='recoMuonRPClayerMask'), + nValidPixelHits = Var('nValidPixelHits', 'int', doc='nValidPixelHits'), + nValidStripHits = Var('nValidStripHits', 'int', doc='nValidStripHits'), + nPixelLayersWithMeasurement = Var('nPixelLayersWithMeasurement', 'int', doc='nPixelLayersWithMeasurement'), + nTrackerLayersWithMeasurement = Var('nTrackerLayersWithMeasurement', 'int', doc='nTrackerLayersWithMeasurement'), + trk_chi2 = Var('trk_chi2', 'float', precision=14, doc='trk_chi2'), + trk_ndof = Var('trk_ndof', 'float', precision=14, doc='trk_ndof'), + trk_dxy = Var('trk_dxy', 'float', precision=14, doc='trk_dxy'), + trk_dz = Var('trk_dz', 'float', precision=14, doc='trk_dz'), + trk_qoverp = Var('trk_qoverp', 'float', precision=14, doc='trk_qoverp'), + trk_lambda = Var('trk_lambda', 'float', precision=14, doc='trk_lambda'), + trk_pt = Var('trk_pt', 'float', precision=14, doc='trk_pt'), + trk_phi = Var('trk_phi', 'float', precision=14, doc='trk_phi'), + trk_eta = Var('trk_eta', 'float', precision=14, doc='trk_eta'), + trk_dxyError = Var('trk_dxyError', 'float', precision=14, doc='trk_dxyError'), + trk_dzError = Var('trk_dzError', 'float', precision=14, doc='trk_dzError'), + trk_qoverpError = Var('trk_qoverpError', 'float', precision=14, doc='trk_qoverpError'), + trk_lambdaError = Var('trk_lambdaError', 'float', precision=14, doc='trk_lambdaError'), + trk_phiError = Var('trk_phiError', 'float', precision=14, doc='trk_phiError'), + trk_dsz = Var('trk_dsz', 'float', precision=14, doc='trk_dsz'), + trk_dszError = Var('trk_dszError', 'float', precision=14, doc='trk_dszError'), + trk_qoverp_lambda_cov = Var('trk_qoverp_lambda_cov', 'float', precision=14, doc='trk_qoverp_lambda_cov'), + trk_qoverp_phi_cov = Var('trk_qoverp_phi_cov', 'float', precision=14, doc='trk_qoverp_phi_cov'), + trk_qoverp_dxy_cov = Var('trk_qoverp_dxy_cov', 'float', precision=14, doc='trk_qoverp_dxy_cov'), + trk_qoverp_dsz_cov = Var('trk_qoverp_dsz_cov', 'float', precision=14, doc='trk_qoverp_dsz_cov'), + trk_lambda_phi_cov = Var('trk_lambda_phi_cov', 'float', precision=14, doc='trk_lambda_phi_cov'), + trk_lambda_dxy_cov = Var('trk_lambda_dxy_cov', 'float', precision=14, doc='trk_lambda_dxy_cov'), + trk_lambda_dsz_cov = Var('trk_lambda_dsz_cov', 'float', precision=14, doc='trk_lambda_dsz_cov'), + trk_phi_dxy_cov = Var('trk_phi_dxy_cov', 'float', precision=14, doc='trk_phi_dxy_cov'), + trk_phi_dsz_cov = Var('trk_phi_dsz_cov', 'float', precision=14, doc='trk_phi_dsz_cov'), + trk_dxy_dsz_cov = Var('trk_dxy_dsz_cov', 'float', precision=14, doc='trk_dxy_dsz_cov'), + trk_vx = Var('trk_vx', 'float', precision=14, doc='trk_vx'), + trk_vy = Var('trk_vy', 'float', precision=14, doc='trk_vy'), + trk_vz = Var('trk_vz', 'float', precision=14, doc='trk_vz'), + ) + ) + + process.trackScoutingTable = cms.EDProducer("SimpleRun3ScoutingTrackFlatTableProducer", + src = cms.InputTag("hltScoutingTrackPacker"), + cut = cms.string(""), + name = cms.string("ScoutingTrack"), + doc = cms.string("Track Scouting Informations"), + singleton = cms.bool(False), # the number of entries is variable + extension = cms.bool(False), # this is the main table for the muons + variables = cms.PSet( + pt = Var('tk_pt', 'float', precision=14, doc='pt'), + eta = Var('tk_eta', 'float', precision=14, doc='eta'), + phi = Var('tk_phi', 'float', precision=14, doc='phi'), + chi2 = Var('tk_chi2', 'float', precision=14, doc='chi2'), + ndof = Var('tk_ndof', 'float', precision=14, doc='ndof'), + charge = Var('tk_charge', 'int', doc='charge'), + dxy = Var('tk_dxy', 'float', precision=14, doc='dxy'), + dz = Var('tk_dz', 'float', precision=14, doc='dz'), + nValidPixelHits = Var('tk_nValidPixelHits', 'int', doc='nValidPixelHits'), + nValidStripHits = Var('tk_nValidStripHits', 'int', doc='nValidStripHits'), + nTrackerLayersWithMeasurement = Var('tk_nTrackerLayersWithMeasurement', 'int', doc='nTrackerLayersWithMeasurement'), + qoverp = Var('tk_qoverp', 'float', precision=14, doc='qoverp'), + _lambda = Var('tk_lambda', 'float', precision=14, doc='lambda'), + dxyError = Var('tk_dxy_Error', 'float', precision=14, doc='dxyError'), + dzError = Var('tk_dz_Error', 'float', precision=14, doc='dzError'), + qoverpError = Var('tk_qoverp_Error', 'float', precision=14, doc='qoverpError'), + lambdaError = Var('tk_lambda_Error', 'float', precision=14, doc='lambdaError'), + phiError = Var('tk_phi_Error', 'float', precision=14, doc='phiError'), + dsz = Var('tk_dsz', 'float', precision=14, doc='dsz'), + dszError = Var('tk_dsz_Error', 'float', precision=14, doc='dszError'), + qoverp_lambda_cov = Var('tk_qoverp_lambda_cov', 'float', precision=14, doc='qoverp_lambda_cov'), + qoverp_phi_cov = Var('tk_qoverp_phi_cov', 'float', precision=14, doc='qoverp_phi_cov'), + qoverp_dxy_cov = Var('tk_qoverp_dxy_cov', 'float', precision=14, doc='qoverp_dxy_cov'), + qoverp_dsz_cov = Var('tk_qoverp_dsz_cov', 'float', precision=14, doc='qoverp_dsz_cov'), + lambda_phi_cov = Var('tk_lambda_phi_cov', 'float', precision=14, doc='lambda_phi_cov'), + lambda_dxy_cov = Var('tk_lambda_dxy_cov', 'float', precision=14, doc='lambda_dxy_cov'), + lambda_dsz_cov = Var('tk_lambda_dsz_cov', 'float', precision=14, doc='lambda_dsz_cov'), + phi_dxy_cov = Var('tk_phi_dxy_cov', 'float', precision=14, doc='phi_dxy_cov'), + phi_dsz_cov = Var('tk_phi_dsz_cov', 'float', precision=14, doc='phi_dsz_cov'), + dxy_dsz_cov = Var('tk_dxy_dsz_cov', 'float', precision=14, doc='dxy_dsz_cov'), + vtxInd = Var('tk_vtxInd', 'int', doc='vtxInd'), + vx = Var('tk_vx', 'float', precision=14, doc='vx'), + vy = Var('tk_vy', 'float', precision=14, doc='vy'), + vz = Var('tk_vz', 'float', precision=14, doc='vz'), + ) + ) + + process.primaryvertexScoutingTable = cms.EDProducer("SimpleRun3ScoutingVertexFlatTableProducer", + src = cms.InputTag("hltScoutingPrimaryVertexPacker", "primaryVtx"), + cut = cms.string(""), + name = cms.string("ScoutingPrimaryVertex"), + doc = cms.string("PrimaryVertex Scouting Informations"), + singleton = cms.bool(False), # the number of entries is variable + extension = cms.bool(False), # this is the main table for the muons + variables = cms.PSet( + x = Var('x', 'float', precision=14, doc='x'), + y = Var('y', 'float', precision=14, doc='y'), + z = Var('z', 'float', precision=14, doc='z'), + xError = Var('xError', 'float', precision=14, doc='xError'), + yError = Var('yError', 'float', precision=14, doc='yError'), + zError = Var('zError', 'float', precision=14, doc='zError'), + tracksSize = Var('tracksSize', 'int', doc='tracksSize'), + chi2 = Var('chi2', 'float', precision=14, doc='chi2'), + ndof = Var('ndof', 'int', doc='ndof'), + isValidVtx = Var('isValidVtx', 'bool', doc='isValidVtx'), + ) + ) + + process.displacedvertexScoutingTable = cms.EDProducer("SimpleRun3ScoutingVertexFlatTableProducer", + src = cms.InputTag("hltScoutingMuonPacker","displacedVtx"), + cut = cms.string(""), + name = cms.string("ScoutingDisplacedVertex"), + doc = cms.string("DisplacedVertex Scouting Informations"), + singleton = cms.bool(False), # the number of entries is variable + extension = cms.bool(False), # this is the main table for the muons + variables = cms.PSet( + x = Var('x', 'float', precision=14, doc='x'), + y = Var('y', 'float', precision=14, doc='y'), + z = Var('z', 'float', precision=14, doc='z'), + xError = Var('xError', 'float', precision=14, doc='xError'), + yError = Var('yError', 'float', precision=14, doc='yError'), + zError = Var('zError', 'float', precision=14, doc='zError'), + tracksSize = Var('tracksSize', 'int', doc='tracksSize'), + chi2 = Var('chi2', 'float', precision=14, doc='chi2'), + ndof = Var('ndof', 'int', doc='ndof'), + isValidVtx = Var('isValidVtx', 'bool', doc='isValidVtx'), + ) + ) + + process.rhoScoutingTable = cms.EDProducer("GlobalVariablesTableProducer", + name = cms.string(""), + variables = cms.PSet( + ScoutingRho = ExtVar( cms.InputTag("hltScoutingPFPacker", "rho"), "double", doc = "rho from all PF Candidates, no foreground removal (for isolation of prompt photons)" ), + ) + ) + + process.metScoutingTable = cms.EDProducer("GlobalVariablesTableProducer", + name = cms.string("ScoutingMET"), + variables = cms.PSet( + pt = ExtVar( cms.InputTag("hltScoutingPFPacker", "pfMetPt"), "double", doc = "Scouting MET pt"), + phi = ExtVar( cms.InputTag("hltScoutingPFPacker", "pfMetPhi"), "double", doc = "Scouting MET phi"), + ) + ) + + run3ScoutingTaskName = "run3ScoutingTask" + setattr(process, run3ScoutingTaskName, cms.Task( + getattr(process,"photonScoutingTable"), + getattr(process,"muonScoutingTable"), + getattr(process,"electronScoutingTable"), + getattr(process,"trackScoutingTable"), + getattr(process,"primaryvertexScoutingTable"), + getattr(process,"displacedvertexScoutingTable"), + getattr(process,"rhoScoutingTable"), + getattr(process,"metScoutingTable"), + ) + ) + process.customRun3ScoutingNanoAODTask.add(getattr(process,run3ScoutingTaskName)) + return process + +def ConvertScoutingToReco(process): + """ + Convert Run 3 scouting particles to recoPFCandidates + """ + process.pfcands = cms.EDProducer( + "Run3ScoutingParticleToRecoPFCandidateProducer", + scoutingparticle=cms.InputTag("hltScoutingPFPacker"), + ) + + process.customRun3ScoutingNanoAODTask.add(cms.Task(getattr(process,"pfcands"))) + return process + +def AddParticles(process): + """ + Add Run 3 scouting particles + """ + process.particleTable = cms.EDProducer("SimpleCandidateFlatTableProducer", + src = cms.InputTag("pfcands"), + name = cms.string("ScoutingParticle"), + cut = cms.string(""), + doc = cms.string("ScoutingParticle"), + singleton = cms.bool(False), + extension = cms.bool(False), # this is the main table + externalVariables = cms.PSet( + vertexIndex = ExtVar(cms.InputTag("pfcands", "vertexIndex"), int, doc="vertex index"), + trkNormchi2 = ExtVar(cms.InputTag("pfcands", "normchi2"), float, doc="normchi2 of best track", precision=6), + trkDz = ExtVar(cms.InputTag("pfcands", "dz"), float, doc="dz of best track", precision=6), + trkDxy = ExtVar(cms.InputTag("pfcands", "dxy"), float, doc="dxy of best track", precision=6), + trkDzsig = ExtVar(cms.InputTag("pfcands", "dzsig"), float, doc="dzsig of best track", precision=6), + trkDxysig = ExtVar(cms.InputTag("pfcands", "dxysig"), float, doc="dxysig of best track", precision=6), + trkLostInnerHits = ExtVar(cms.InputTag("pfcands", "lostInnerHits"), int, doc="lostInnerHits of best track"), + trkQuality = ExtVar(cms.InputTag("pfcands", "quality"), int, doc="quality of best track"), + trkPt = ExtVar(cms.InputTag("pfcands", "trkPt"), float, doc="pt of best track", precision=6), + trkEta = ExtVar(cms.InputTag("pfcands", "trkEta"), float, doc="eta of best track", precision=6), + trkPhi = ExtVar(cms.InputTag("pfcands", "trkPhi"), float, doc="phi of best track", precision=6), + ), + variables = cms.PSet( + CandVars, + ), + ) + process.customRun3ScoutingNanoAODTask.add(cms.Task(getattr(process,"particleTable"))) + return process + +def AddAK4PFJets(process): + """ + Add Run 3 scouting AK4 PF jets + """ + from RecoJets.JetProducers.ak4PFJets_cfi import ak4PFJets + process.ak4Jets = ak4PFJets.clone( + src = ("pfcands"), + ) + + process.ak4ParticleNetJetTagInfos = cms.EDProducer("DeepBoostedJetTagInfoProducer", + jet_radius = cms.double( 0.4 ), + min_jet_pt = cms.double( 5.0 ), + max_jet_eta = cms.double( 2.5 ), + min_pt_for_track_properties = cms.double( 0.95 ), + min_pt_for_pfcandidates = cms.double( 0.1 ), + use_puppiP4 = cms.bool( False ), + include_neutrals = cms.bool( True ), + sort_by_sip2dsig = cms.bool( False ), + min_puppi_wgt = cms.double( -1.0 ), + flip_ip_sign = cms.bool( False ), + sip3dSigMax = cms.double( -1.0 ), + use_hlt_features = cms.bool( False ), + pf_candidates = cms.InputTag( "pfcands" ), + jets = cms.InputTag( "ak4Jets" ), + puppi_value_map = cms.InputTag( "" ), + use_scouting_features = cms.bool( True ), + normchi2_value_map = cms.InputTag("pfcands", "normchi2"), + dz_value_map = cms.InputTag("pfcands", "dz"), + dxy_value_map = cms.InputTag("pfcands", "dxy"), + dzsig_value_map = cms.InputTag("pfcands", "dzsig"), + dxysig_value_map = cms.InputTag("pfcands", "dxysig"), + lostInnerHits_value_map = cms.InputTag("pfcands", "lostInnerHits"), + quality_value_map = cms.InputTag("pfcands", "quality"), + trkPt_value_map = cms.InputTag("pfcands", "trkPt"), + trkEta_value_map = cms.InputTag("pfcands", "trkEta"), + trkPhi_value_map = cms.InputTag("pfcands", "trkPhi"), + ) + + from RecoBTag.ONNXRuntime.boostedJetONNXJetTagsProducer_cfi import boostedJetONNXJetTagsProducer + process.ak4ParticleNetJetTags = cms.EDProducer("BoostedJetONNXValueMapProducer", + jets = cms.InputTag("ak4Jets"), + src = cms.InputTag("ak4ParticleNetJetTagInfos"), + preprocess_json = cms.string("RecoBTag/Combined/data/Run3Scouting/ParticleNetAK4/V00/preprocess.json"), + model_path = cms.FileInPath("RecoBTag/Combined/data/Run3Scouting/ParticleNetAK4/V00/particle-net.onnx"), + flav_names = cms.vstring(["probb", "probbb","probc", "probcc", "probuds", "probg", "probundef"]), + debugMode = cms.untracked.bool(False), + ) + + process.ak4JetTable = cms.EDProducer("SimpleCandidateFlatTableProducer", + src = cms.InputTag("ak4Jets"), + name = cms.string("ScoutingJet"), + cut = cms.string(""), + doc = cms.string("ScoutingJet"), + singleton = cms.bool(False), + extension = cms.bool(False), # this is the main table + externalVariables = cms.PSet( + particleNet_prob_b = ExtVar(cms.InputTag('ak4ParticleNetJetTags:probb'), float, doc="ParticleNet probability of b", precision=10), + particleNet_prob_bb = ExtVar(cms.InputTag('ak4ParticleNetJetTags:probbb'), float, doc="ParticleNet probability of bb", precision=10), + particleNet_prob_c = ExtVar(cms.InputTag('ak4ParticleNetJetTags:probc'), float, doc="ParticleNet probability of c", precision=10), + particleNet_prob_cc = ExtVar(cms.InputTag('ak4ParticleNetJetTags:probcc'), float, doc="ParticleNet probability of cc", precision=10), + particlenet_prob_uds = ExtVar(cms.InputTag('ak4ParticleNetJetTags:probuds'), float, doc="particlenet probability of uds", precision=10), + particleNet_prob_g = ExtVar(cms.InputTag('ak4ParticleNetJetTags:probg'), float, doc="ParticleNet probability of g", precision=10), + particleNet_prob_undef = ExtVar(cms.InputTag('ak4ParticleNetJetTags:probundef'), float, doc="ParticleNet probability of undef", precision=10), + ), + variables = cms.PSet( + P4Vars, + area = Var("jetArea()", float, doc="jet catchment area, for JECs",precision=10), + chHEF = Var("chargedHadronEnergy()/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="charged Hadron Energy Fraction", precision= 6), + neHEF = Var("neutralHadronEnergy()/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="neutral Hadron Energy Fraction", precision= 6), + chEmEF = Var("(electronEnergy()+muonEnergy())/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="charged Electromagnetic Energy Fraction", precision= 6), + neEmEF = Var("(photonEnergy())/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="neutral Electromagnetic Energy Fraction", precision= 6), + muEmEF = Var("(muonEnergy())/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="muon Energy Fraction", precision= 6), + nCh = Var("chargedHadronMultiplicity()", int, doc="number of charged hadrons in the jet"), + nNh = Var("neutralHadronMultiplicity()", int, doc="number of neutral hadrons in the jet"), + nMuons = Var("muonMultiplicity()", int, doc="number of muons in the jet"), + nElectrons = Var("electronMultiplicity()", int, doc="number of electrons in the jet"), + nPhotons = Var("photonMultiplicity()", int, doc="number of photons in the jet"), + nConstituents = Var("numberOfDaughters()", "uint8", doc="number of particles in the jet") + ), + ) + + ak4PFJetTaskName = "ak4PFJetTask" + setattr(process, ak4PFJetTaskName, cms.Task( + getattr(process,"ak4Jets"), + getattr(process,"ak4ParticleNetJetTagInfos"), + getattr(process,"ak4ParticleNetJetTags"), + getattr(process,"ak4JetTable"), + ) + ) + process.customRun3ScoutingNanoAODTask.add(getattr(process,ak4PFJetTaskName)) + + return process + +def AddAK8PFJets(process): + """ + Add Run 3 scouting AK8 PF jets + """ + from RecoJets.JetProducers.ak4PFJets_cfi import ak4PFJets + process.ak8Jets = ak4PFJets.clone( + src = ("pfcands"), + rParam = 0.8, + jetPtMin = 50.0, + ) + + process.ak8JetsSoftDrop = ak4PFJets.clone( + src = ("pfcands"), + rParam = 0.8, + jetPtMin = 50.0, + useSoftDrop = cms.bool(True), + zcut = cms.double(0.1), + beta = cms.double(0.0), + R0 = cms.double(0.8), + useExplicitGhosts = cms.bool(True), + writeCompound = cms.bool(True), + jetCollInstanceName=cms.string("SubJets"), + ) + + process.ak8JetsSoftDropMass = cms.EDProducer("RecoJetDeltaRValueMapProducer", + src = cms.InputTag("ak8Jets"), + matched = cms.InputTag("ak8JetsSoftDrop"), + distMax = cms.double(0.8), + value = cms.string('mass') + ) + + from RecoJets.JetProducers.ECF_cff import ecfNbeta1 + process.ecfNbeta1 = ecfNbeta1.clone(src = cms.InputTag("ak8Jets"), srcWeights="") + + from RecoJets.JetProducers.nJettinessAdder_cfi import Njettiness + process.Njettiness = Njettiness.clone(src = cms.InputTag("ak8Jets"), srcWeights="") + + process.ak8ParticleNetJetTagInfos = cms.EDProducer("DeepBoostedJetTagInfoProducer", + jet_radius = cms.double( 0.8 ), + min_jet_pt = cms.double( 5.0 ), + max_jet_eta = cms.double( 2.5 ), + min_pt_for_track_properties = cms.double( 0.95 ), + min_pt_for_pfcandidates = cms.double( 0.1 ), + use_puppiP4 = cms.bool( False ), + include_neutrals = cms.bool( True ), + sort_by_sip2dsig = cms.bool( False ), + min_puppi_wgt = cms.double( -1.0 ), + flip_ip_sign = cms.bool( False ), + sip3dSigMax = cms.double( -1.0 ), + use_hlt_features = cms.bool( False ), + pf_candidates = cms.InputTag( "pfcands" ), + jets = cms.InputTag( "ak8Jets" ), + puppi_value_map = cms.InputTag( "" ), + use_scouting_features = cms.bool( True ), + normchi2_value_map = cms.InputTag("pfcands", "normchi2"), + dz_value_map = cms.InputTag("pfcands", "dz"), + dxy_value_map = cms.InputTag("pfcands", "dxy"), + dzsig_value_map = cms.InputTag("pfcands", "dzsig"), + dxysig_value_map = cms.InputTag("pfcands", "dxysig"), + lostInnerHits_value_map = cms.InputTag("pfcands", "lostInnerHits"), + quality_value_map = cms.InputTag("pfcands", "quality"), + trkPt_value_map = cms.InputTag("pfcands", "trkPt"), + trkEta_value_map = cms.InputTag("pfcands", "trkEta"), + trkPhi_value_map = cms.InputTag("pfcands", "trkPhi"), + ) + + from RecoBTag.ONNXRuntime.boostedJetONNXJetTagsProducer_cfi import boostedJetONNXJetTagsProducer + process.ak8ParticleNetJetTags = cms.EDProducer("BoostedJetONNXValueMapProducer", + jets = cms.InputTag("ak8Jets"), + src = cms.InputTag("ak8ParticleNetJetTagInfos"), + preprocess_json = cms.string("RecoBTag/Combined/data/Run3Scouting/ParticleNetAK8/General/V00/preprocess.json"), + model_path = cms.FileInPath("RecoBTag/Combined/data/Run3Scouting/ParticleNetAK8/General/V00/particle-net.onnx"), + flav_names = cms.vstring(["probHbb", "probHcc","probHqq", "probQCDall"]), + debugMode = cms.untracked.bool(False), + ) + + process.ak8ParticleNetMassRegressionJetTags = cms.EDProducer("BoostedJetONNXValueMapProducer", + jets = cms.InputTag("ak8Jets"), + src = cms.InputTag("ak8ParticleNetJetTagInfos"), + preprocess_json = cms.string("RecoBTag/Combined/data/Run3Scouting/ParticleNetAK8/MassRegression/V00/preprocess.json"), + model_path = cms.FileInPath("RecoBTag/Combined/data/Run3Scouting/ParticleNetAK8/MassRegression/V00/particle-net.onnx"), + flav_names = cms.vstring(["mass"]), + debugMode = cms.untracked.bool(False), + ) + + process.ak8JetTable = cms.EDProducer("SimpleCandidateFlatTableProducer", + src = cms.InputTag("ak8Jets"), + name = cms.string("ScoutingFatJet"), + cut = cms.string(""), + doc = cms.string("ScoutingFatJet"), + singleton = cms.bool(False), + extension = cms.bool(False), # this is the main table + externalVariables = cms.PSet( + msoftdrop = ExtVar(cms.InputTag('ak8JetsSoftDropMass'), float, doc="Softdrop mass", precision=10), + n2b1 = ExtVar(cms.InputTag('ecfNbeta1:ecfN2'), float, doc="N2 with beta=1", precision=10), + n3b1 = ExtVar(cms.InputTag('ecfNbeta1:ecfN3'), float, doc="N3 with beta=1", precision=10), + tau1 = ExtVar(cms.InputTag('Njettiness:tau1'), float, doc="Nsubjettiness (1 axis)", precision=10), + tau2 = ExtVar(cms.InputTag('Njettiness:tau2'), float, doc="Nsubjettiness (2 axis)", precision=10), + tau3 = ExtVar(cms.InputTag('Njettiness:tau3'), float, doc="Nsubjettiness (3 axis)", precision=10), + tau4 = ExtVar(cms.InputTag('Njettiness:tau4'), float, doc="Nsubjettiness (4 axis)", precision=10), + particleNet_mass = ExtVar(cms.InputTag('ak8ParticleNetMassRegressionJetTags:mass'), float, doc="ParticleNet regressed mass", precision=10), + particleNet_prob_Hbb = ExtVar(cms.InputTag('ak8ParticleNetJetTags:probHbb'), float, doc="ParticleNet probability of Hbb", precision=10), + particleNet_prob_Hcc = ExtVar(cms.InputTag('ak8ParticleNetJetTags:probHcc'), float, doc="ParticleNet probability of Hcc", precision=10), + particleNet_prob_Hqq = ExtVar(cms.InputTag('ak8ParticleNetJetTags:probHqq'), float, doc="ParticleNet probability of Hqq", precision=10), + particleNet_prob_QCD = ExtVar(cms.InputTag('ak8ParticleNetJetTags:probQCDall'), float, doc="ParticleNet probability of QCD", precision=10), + ), + variables = cms.PSet( + P4Vars, + area = Var("jetArea()", float, doc="jet catchment area, for JECs",precision=10), + chHEF = Var("chargedHadronEnergy()/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="charged Hadron Energy Fraction", precision= 6), + neHEF = Var("neutralHadronEnergy()/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="neutral Hadron Energy Fraction", precision= 6), + chEmEF = Var("(electronEnergy()+muonEnergy())/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="charged Electromagnetic Energy Fraction", precision= 6), + neEmEF = Var("(photonEnergy())/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="neutral Electromagnetic Energy Fraction", precision= 6), + muEmEF = Var("(muonEnergy())/(chargedHadronEnergy()+neutralHadronEnergy()+photonEnergy()+electronEnergy()+muonEnergy())", float, doc="muon Energy Fraction", precision= 6), + nCh = Var("chargedHadronMultiplicity()", int, doc="number of charged hadrons in the jet"), + nNh = Var("neutralHadronMultiplicity()", int, doc="number of neutral hadrons in the jet"), + nMuons = Var("muonMultiplicity()", int, doc="number of muons in the jet"), + nElectrons = Var("electronMultiplicity()", int, doc="number of electrons in the jet"), + nPhotons = Var("photonMultiplicity()", int, doc="number of photons in the jet"), + nConstituents = Var("numberOfDaughters()", "uint8", doc="number of particles in the jet") + ), + ) + + ak8PFJetTaskName = "ak8PFJetTask" + setattr(process, ak8PFJetTaskName, cms.Task( + getattr(process,"ak8Jets"), + getattr(process,"ak8JetsSoftDrop"), + getattr(process,"ak8JetsSoftDropMass"), + getattr(process,"ecfNbeta1"), + getattr(process,"Njettiness"), + getattr(process,"ak8ParticleNetJetTagInfos"), + getattr(process,"ak8ParticleNetJetTags"), + getattr(process,"ak8ParticleNetMassRegressionJetTags"), + getattr(process,"ak8JetTable"), + ) + ) + process.customRun3ScoutingNanoAODTask.add(getattr(process,ak8PFJetTaskName)) + + return process + +def AddAK4MatchToGen(process): + process.ak4MatchGen = cms.EDProducer("RecoJetToGenJetDeltaRValueMapProducer", + src = cms.InputTag("ak4Jets"), + matched = cms.InputTag("slimmedGenJets"), + distMax = cms.double(0.4), + value = cms.string("index"), + ) + ak4MatchGenTaskName = "ak4MatchGenTask" + setattr(process, ak4MatchGenTaskName, cms.Task(process.ak4MatchGen)) + externalVariables = getattr(process.ak4JetTable, 'externalVariables', cms.PSet()) + externalVariables.genJetIdx = ExtVar(cms.InputTag("ak4MatchGen"), int, doc="gen jet idx") + process.ak4JetTable.externalVariables = externalVariables + process.customRun3ScoutingNanoAODTask.add(getattr(process,ak4MatchGenTaskName)) + + return process + +def AddAK8MatchToGen(process): + process.ak8MatchGen = cms.EDProducer("RecoJetToGenJetDeltaRValueMapProducer", + src = cms.InputTag("ak8Jets"), + matched = cms.InputTag("slimmedGenJetsAK8"), + distMax = cms.double(0.8), + value = cms.string("index"), + ) + ak8MatchGenTaskName = "ak8MatchGenTask" + setattr(process, ak8MatchGenTaskName, cms.Task(process.ak8MatchGen)) + externalVariables = getattr(process.ak8JetTable, 'externalVariables', cms.PSet()) + externalVariables.genJetAK8Idx = ExtVar(cms.InputTag("ak8MatchGen"), int, doc="gen jet idx") + process.ak8JetTable.externalVariables = externalVariables + process.customRun3ScoutingNanoAODTask.add(getattr(process,ak8MatchGenTaskName)) + + return process + +#=========================================================================== +# +# CUSTOMISATION function +# +#=========================================================================== +def PrepRun3ScoutingCustomNanoAOD(process,runOnMC): + process.customRun3ScoutingNanoAODTask = cms.Task() + process = ConvertScoutingToReco(process) + process = AddRun3Scouting(process) + process = AddParticles(process) + process = AddAK4PFJets(process) + process = AddAK8PFJets(process) + + if (runOnMC): + process = AddAK4MatchToGen(process) + process = AddAK8MatchToGen(process) + + process.schedule.associate(process.customRun3ScoutingNanoAODTask) + + return process + +def PrepRun3ScoutingCustomNanoAOD_MC(process): + process = PrepRun3ScoutingCustomNanoAOD(process,runOnMC=True) + + return process + +def PrepRun3ScoutingCustomNanoAOD_Data(process): + process = PrepRun3ScoutingCustomNanoAOD(process,runOnMC=False) + + return process diff --git a/RecoBTag/FeatureTools/plugins/DeepBoostedJetTagInfoProducer.cc b/RecoBTag/FeatureTools/plugins/DeepBoostedJetTagInfoProducer.cc index f04ce9e473caa..f83218b91750b 100644 --- a/RecoBTag/FeatureTools/plugins/DeepBoostedJetTagInfoProducer.cc +++ b/RecoBTag/FeatureTools/plugins/DeepBoostedJetTagInfoProducer.cc @@ -96,6 +96,30 @@ class DeepBoostedJetTagInfoProducer : public edm::stream::EDProducer<> { const static int min_valid_pixel_hits_; std::map puppi_wgt_cache; + + const static std::vector particle_features_scouting_; + const bool use_scouting_features_; + edm::EDGetTokenT> normchi2_value_map_token_; + edm::EDGetTokenT> dz_value_map_token_; + edm::EDGetTokenT> dxy_value_map_token_; + edm::EDGetTokenT> dzsig_value_map_token_; + edm::EDGetTokenT> dxysig_value_map_token_; + edm::EDGetTokenT> lostInnerHits_value_map_token_; + edm::EDGetTokenT> quality_value_map_token_; + edm::EDGetTokenT> trkPt_value_map_token_; + edm::EDGetTokenT> trkEta_value_map_token_; + edm::EDGetTokenT> trkPhi_value_map_token_; + + edm::Handle> normchi2_value_map_; + edm::Handle> dz_value_map_; + edm::Handle> dxy_value_map_; + edm::Handle> dzsig_value_map_; + edm::Handle> dxysig_value_map_; + edm::Handle> lostInnerHits_value_map_; + edm::Handle> quality_value_map_; + edm::Handle> trkPt_value_map_; + edm::Handle> trkEta_value_map_; + edm::Handle> trkPhi_value_map_; }; const std::vector DeepBoostedJetTagInfoProducer::particle_features_{ @@ -139,6 +163,15 @@ const std::vector DeepBoostedJetTagInfoProducer::particle_features_ "jet_pfcand_puppiw", "pfcand_mask"}; +const std::vector DeepBoostedJetTagInfoProducer::particle_features_scouting_{ + "pfcand_quality", "pfcand_charge", "pfcand_isEl", "pfcand_isMu", + "pfcand_isChargedHad", "pfcand_isGamma", "pfcand_isNeutralHad", "pfcand_phirel", + "pfcand_etarel", "pfcand_deltaR", "pfcand_abseta", "pfcand_ptrel_log", + "pfcand_erel_log", "pfcand_pt_log", "pfcand_normchi2", "pfcand_dz", + "pfcand_dxy", "pfcand_dxysig", "pfcand_btagEtaRel", "pfcand_btagPtRatio", + "pfcand_btagPParRatio", "pfcand_mask", "pfcand_pt_log_nopuppi", "pfcand_dzsig", + "pfcand_e_log_nopuppi", "pfcand_ptrel", "pfcand_erel", "pfcand_lostInnerHits"}; + const std::vector DeepBoostedJetTagInfoProducer::sv_features_{"sv_mask", "sv_ptrel", "sv_erel", @@ -195,7 +228,8 @@ DeepBoostedJetTagInfoProducer::DeepBoostedJetTagInfoProducer(const edm::Paramete use_puppi_value_map_(false), use_pvasq_value_map_(false), track_builder_token_( - esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) { + esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))), + use_scouting_features_(iConfig.getParameter("use_scouting_features")) { const auto &puppi_value_map_tag = iConfig.getParameter("puppi_value_map"); if (!puppi_value_map_tag.label().empty()) { puppi_value_map_token_ = consumes>(puppi_value_map_tag); @@ -209,6 +243,56 @@ DeepBoostedJetTagInfoProducer::DeepBoostedJetTagInfoProducer(const edm::Paramete use_pvasq_value_map_ = true; } + const auto &normchi2_value_map_tag = iConfig.getParameter("normchi2_value_map"); + if (!normchi2_value_map_tag.label().empty()) { + normchi2_value_map_token_ = consumes>(normchi2_value_map_tag); + } + + const auto &dz_value_map_tag = iConfig.getParameter("dz_value_map"); + if (!dz_value_map_tag.label().empty()) { + dz_value_map_token_ = consumes>(dz_value_map_tag); + } + + const auto &dxy_value_map_tag = iConfig.getParameter("dxy_value_map"); + if (!dxy_value_map_tag.label().empty()) { + dxy_value_map_token_ = consumes>(dxy_value_map_tag); + } + + const auto &dzsig_value_map_tag = iConfig.getParameter("dzsig_value_map"); + if (!dzsig_value_map_tag.label().empty()) { + dzsig_value_map_token_ = consumes>(dzsig_value_map_tag); + } + + const auto &dxysig_value_map_tag = iConfig.getParameter("dxysig_value_map"); + if (!dxysig_value_map_tag.label().empty()) { + dxysig_value_map_token_ = consumes>(dxysig_value_map_tag); + } + + const auto &lostInnerHits_value_map_tag = iConfig.getParameter("lostInnerHits_value_map"); + if (!lostInnerHits_value_map_tag.label().empty()) { + lostInnerHits_value_map_token_ = consumes>(lostInnerHits_value_map_tag); + } + + const auto &quality_value_map_tag = iConfig.getParameter("quality_value_map"); + if (!quality_value_map_tag.label().empty()) { + quality_value_map_token_ = consumes>(quality_value_map_tag); + } + + const auto &trkPt_value_map_tag = iConfig.getParameter("trkPt_value_map"); + if (!trkPt_value_map_tag.label().empty()) { + trkPt_value_map_token_ = consumes>(trkPt_value_map_tag); + } + + const auto &trkEta_value_map_tag = iConfig.getParameter("trkEta_value_map"); + if (!trkEta_value_map_tag.label().empty()) { + trkEta_value_map_token_ = consumes>(trkEta_value_map_tag); + } + + const auto &trkPhi_value_map_tag = iConfig.getParameter("trkPhi_value_map"); + if (!trkPhi_value_map_tag.label().empty()) { + trkPhi_value_map_token_ = consumes>(trkPhi_value_map_tag); + } + produces(); } @@ -235,6 +319,17 @@ void DeepBoostedJetTagInfoProducer::fillDescriptions(edm::ConfigurationDescripti desc.add("jets", edm::InputTag("ak8PFJetsPuppi")); desc.add("puppi_value_map", edm::InputTag("puppi")); desc.add("vertex_associator", edm::InputTag("primaryVertexAssociation", "original")); + desc.add("use_scouting_features", false); + desc.add("normchi2_value_map", edm::InputTag("")); + desc.add("dz_value_map", edm::InputTag("")); + desc.add("dxy_value_map", edm::InputTag("")); + desc.add("dzsig_value_map", edm::InputTag("")); + desc.add("dxysig_value_map", edm::InputTag("")); + desc.add("lostInnerHits_value_map", edm::InputTag("")); + desc.add("quality_value_map", edm::InputTag("")); + desc.add("trkPt_value_map", edm::InputTag("")); + desc.add("trkEta_value_map", edm::InputTag("")); + desc.add("trkPhi_value_map", edm::InputTag("")); descriptions.add("pfDeepBoostedJetTagInfos", desc); } @@ -244,23 +339,25 @@ void DeepBoostedJetTagInfoProducer::produce(edm::Event &iEvent, const edm::Event // Input jets auto jets = iEvent.getHandle(jet_token_); // Primary vertexes - iEvent.getByToken(vtx_token_, vtxs_); - if (vtxs_->empty()) { - // produce empty TagInfos in case no primary vertex - iEvent.put(std::move(output_tag_infos)); - return; // exit event + if (!use_scouting_features_) { + iEvent.getByToken(vtx_token_, vtxs_); + if (vtxs_->empty()) { + // produce empty TagInfos in case no primary vertex + iEvent.put(std::move(output_tag_infos)); + return; // exit event + } + // Leading vertex + pv_ = &vtxs_->at(0); + // Secondary vertexs + iEvent.getByToken(sv_token_, svs_); + // Track builder + track_builder_ = iSetup.getHandle(track_builder_token_); } - // Leading vertex - pv_ = &vtxs_->at(0); - // Secondary vertexs - iEvent.getByToken(sv_token_, svs_); // PF candidates iEvent.getByToken(pfcand_token_, pfcands_); is_packed_pf_candidate_collection_ = false; if (dynamic_cast(pfcands_.product())) is_packed_pf_candidate_collection_ = true; - // Track builder - track_builder_ = iSetup.getHandle(track_builder_token_); // Puppi weight value map if (use_puppi_value_map_) iEvent.getByToken(puppi_value_map_token_, puppi_value_map_); @@ -269,6 +366,18 @@ void DeepBoostedJetTagInfoProducer::produce(edm::Event &iEvent, const edm::Event iEvent.getByToken(pvasq_value_map_token_, pvasq_value_map_); iEvent.getByToken(pvas_token_, pvas_); } + if (use_scouting_features_) { + iEvent.getByToken(normchi2_value_map_token_, normchi2_value_map_); + iEvent.getByToken(dz_value_map_token_, dz_value_map_); + iEvent.getByToken(dxy_value_map_token_, dxy_value_map_); + iEvent.getByToken(dzsig_value_map_token_, dzsig_value_map_); + iEvent.getByToken(dxysig_value_map_token_, dxysig_value_map_); + iEvent.getByToken(lostInnerHits_value_map_token_, lostInnerHits_value_map_); + iEvent.getByToken(quality_value_map_token_, quality_value_map_); + iEvent.getByToken(trkPt_value_map_token_, trkPt_value_map_); + iEvent.getByToken(trkEta_value_map_token_, trkEta_value_map_); + iEvent.getByToken(trkPhi_value_map_token_, trkPhi_value_map_); + } // Loop over jet for (std::size_t jet_n = 0; jet_n < jets->size(); jet_n++) { @@ -277,37 +386,47 @@ void DeepBoostedJetTagInfoProducer::produce(edm::Event &iEvent, const edm::Event // create jet features DeepBoostedJetFeatures features; - if (not use_hlt_features_) { - for (const auto &name : particle_features_) { + if (use_hlt_features_) { + // declare all the feature variables (init as empty vector) + for (const auto &name : particle_features_hlt_) { features.add(name); } - for (const auto &name : sv_features_) { + for (const auto &name : sv_features_hlt_) { + features.add(name); + } + } else if (use_scouting_features_) { + for (const auto &name : particle_features_scouting_) { features.add(name); } } else { - // declare all the feature variables (init as empty vector) - for (const auto &name : particle_features_hlt_) { + for (const auto &name : particle_features_) { features.add(name); } - for (const auto &name : sv_features_hlt_) { + for (const auto &name : sv_features_) { features.add(name); } } // fill values only if above pt threshold and has daughters, otherwise left bool fill_vars = true; - if (jet.pt() < min_jet_pt_ or std::abs(jet.eta()) > max_jet_eta_) + if (jet.pt() < min_jet_pt_ or std::abs(jet.eta()) > max_jet_eta_) { fill_vars = false; - if (jet.numberOfDaughters() == 0) + } + if (jet.numberOfDaughters() == 0 and !use_scouting_features_) { fill_vars = false; + } // fill features if (fill_vars) { fillParticleFeatures(features, jet); - fillSVFeatures(features, jet); + if (!use_scouting_features_) { + fillSVFeatures(features, jet); + } if (use_hlt_features_) { features.check_consistency(particle_features_hlt_); features.check_consistency(sv_features_hlt_); + } else if (use_scouting_features_) { + features.check_consistency(particle_features_scouting_); } else { features.check_consistency(particle_features_); features.check_consistency(sv_features_); @@ -347,14 +466,20 @@ void DeepBoostedJetTagInfoProducer::fillParticleFeatures(DeepBoostedJetFeatures TVector3 jet_direction(jet.momentum().Unit().x(), jet.momentum().Unit().y(), jet.momentum().Unit().z()); GlobalVector jet_ref_track_dir(jet.px(), jet.py(), jet.pz()); const float etasign = jet.eta() > 0 ? 1 : -1; - // vertexes - reco::VertexRefProd PVRefProd(vtxs_); - // track builder - TrackInfoBuilder trackinfo(track_builder_); + std::unique_ptr PVRefProd; + std::unique_ptr trackinfo; + if (not use_scouting_features_) { + // vertexes + PVRefProd = std::make_unique(vtxs_); + // track builder + trackinfo = std::make_unique(track_builder_); + } // make list of pf-candidates to be considered std::vector daughters; + int daui = 0; for (const auto &dau : jet.daughterPtrVector()) { + daui++; // remove particles w/ extremely low puppi weights // [Note] use jet daughters here to get the puppiWgt correctly if ((puppiWgt(dau)) < min_puppi_wgt_) @@ -369,8 +494,8 @@ void DeepBoostedJetTagInfoProducer::fillParticleFeatures(DeepBoostedJetFeatures continue; // only when computing the nagative tagger: remove charged candidates with high sip3d if (flip_ip_sign_ and cand->charge()) { - trackinfo.buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_); - if (trackinfo.getTrackSip3dSig() > max_sip3dsig_) + (*trackinfo).buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_); + if ((*trackinfo).getTrackSip3dSig() > max_sip3dsig_) continue; } // filling daughters @@ -382,9 +507,9 @@ void DeepBoostedJetTagInfoProducer::fillParticleFeatures(DeepBoostedJetFeatures if (sort_by_sip2dsig_) { // sort charged pf candidates by 2d impact parameter significance for (const auto &cand : daughters) { - trackinfo.buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_); + (*trackinfo).buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_); c_sorted.emplace_back(cand, - trackinfo.getTrackSip2dSig(), + (*trackinfo).getTrackSip2dSig(), -btagbtvdeep::mindrsvpfcand(*svs_, &(*cand), jet_radius_), cand->pt() / jet.pt()); } @@ -409,6 +534,9 @@ void DeepBoostedJetTagInfoProducer::fillParticleFeatures(DeepBoostedJetFeatures if (use_hlt_features_) { for (const auto &name : particle_features_hlt_) fts.reserve(name, daughters.size()); + } else if (use_scouting_features_) { + for (const auto &name : particle_features_scouting_) + fts.reserve(name, daughters.size()); } else { for (const auto &name : particle_features_) fts.reserve(name, daughters.size()); @@ -417,20 +545,22 @@ void DeepBoostedJetTagInfoProducer::fillParticleFeatures(DeepBoostedJetFeatures // build white list of candidates i.e. particles and tracks belonging to SV. Needed when input particles are not packed PF candidates std::vector whiteListSV; std::vector whiteListTk; - if (not is_packed_pf_candidate_collection_) { - for (size_t isv = 0; isv < svs_->size(); isv++) { - for (size_t icand = 0; icand < svs_->at(isv).numberOfSourceCandidatePtrs(); icand++) { - const edm::Ptr &cand = svs_->at(isv).sourceCandidatePtr(icand); - if (cand.id() == pfcands_.id()) - whiteListSV.push_back(cand.key()); - } - for (auto cand = svs_->at(isv).begin(); cand != svs_->at(isv).end(); cand++) { - const reco::RecoChargedCandidate *chCand = dynamic_cast(&(*cand)); - if (chCand != nullptr) { - whiteListTk.push_back(chCand->track()); - } - } - } + if (!use_scouting_features_) { + if (not is_packed_pf_candidate_collection_) { + for (size_t isv = 0; isv < svs_->size(); isv++) { + for (size_t icand = 0; icand < svs_->at(isv).numberOfSourceCandidatePtrs(); icand++) { + const edm::Ptr &cand = svs_->at(isv).sourceCandidatePtr(icand); + if (cand.id() == pfcands_.id()) + whiteListSV.push_back(cand.key()); + } + for (auto cand = svs_->at(isv).begin(); cand != svs_->at(isv).end(); cand++) { + const reco::RecoChargedCandidate *chCand = dynamic_cast(&(*cand)); + if (chCand != nullptr) { + whiteListTk.push_back(chCand->track()); + } + } + } + } } // Build observables @@ -455,279 +585,319 @@ void DeepBoostedJetTagInfoProducer::fillParticleFeatures(DeepBoostedJetFeatures auto candP4 = use_puppiP4_ ? puppi_wgt_cache.at(cand.key()) * cand->p4() : cand->p4(); auto candP3 = use_puppiP4_ ? puppi_wgt_cache.at(cand.key()) * cand->momentum() : cand->momentum(); - // candidate track - const reco::Track *track = nullptr; - if (packed_cand) - track = packed_cand->bestTrack(); - else if (reco_cand and useTrackProperties(reco_cand)) - track = reco_cand->bestTrack(); - - // reco-vertex association - int pv_ass_quality = 0; - reco::VertexRef pv_ass; - float vtx_ass = 0; - if (reco_cand) { - if (use_pvasq_value_map_) { - pv_ass_quality = (*pvasq_value_map_)[cand]; - pv_ass = (*pvas_)[cand]; - vtx_ass = vtx_ass_from_pfcand(*reco_cand, pv_ass_quality, pv_ass); - } else - throw edm::Exception(edm::errors::InvalidReference) << "Vertex association missing"; - } - - // Building offline features - if (not use_hlt_features_) { - // in case of packed candidate - if (packed_cand) { - float hcal_fraction = 0.; - if (packed_cand->pdgId() == 1 or packed_cand->pdgId() == 130) - hcal_fraction = packed_cand->hcalFraction(); - else if (packed_cand->isIsolatedChargedHadron()) - hcal_fraction = packed_cand->rawHcalFraction(); - - fts.fill("pfcand_hcalFrac", hcal_fraction); - fts.fill("pfcand_VTX_ass", packed_cand->pvAssociationQuality()); - fts.fill("pfcand_lostInnerHits", packed_cand->lostInnerHits()); - fts.fill("pfcand_quality", track ? track->qualityMask() : 0); - fts.fill("pfcand_charge", packed_cand->charge()); - fts.fill("pfcand_isEl", std::abs(packed_cand->pdgId()) == 11); - fts.fill("pfcand_isMu", std::abs(packed_cand->pdgId()) == 13); - fts.fill("pfcand_isChargedHad", std::abs(packed_cand->pdgId()) == 211); - fts.fill("pfcand_isGamma", std::abs(packed_cand->pdgId()) == 22); - fts.fill("pfcand_isNeutralHad", std::abs(packed_cand->pdgId()) == 130); - fts.fill("pfcand_dz", ip_sign * packed_cand->dz()); - fts.fill("pfcand_dxy", ip_sign * packed_cand->dxy()); - fts.fill("pfcand_dzsig", track ? ip_sign * packed_cand->dz() / packed_cand->dzError() : 0); - fts.fill("pfcand_dxysig", track ? ip_sign * packed_cand->dxy() / packed_cand->dxyError() : 0); + if (use_scouting_features_) { + + fts.fill("pfcand_lostInnerHits", (*lostInnerHits_value_map_)[cand]); + fts.fill("pfcand_quality", (*quality_value_map_)[cand]); + fts.fill("pfcand_charge", reco_cand->charge()); + fts.fill("pfcand_isEl", std::abs(reco_cand->pdgId()) == 11); + fts.fill("pfcand_isMu", std::abs(reco_cand->pdgId()) == 13); + fts.fill("pfcand_isChargedHad", std::abs(reco_cand->pdgId()) == 211); + fts.fill("pfcand_isGamma", std::abs(reco_cand->pdgId()) == 22); + fts.fill("pfcand_isNeutralHad", std::abs(reco_cand->pdgId()) == 130); + fts.fill("pfcand_dz", (*dz_value_map_)[cand]); + fts.fill("pfcand_dzsig", (*dzsig_value_map_)[cand]); + fts.fill("pfcand_dxy", (*dxy_value_map_)[cand]); + fts.fill("pfcand_dxysig", (*dxysig_value_map_)[cand]); + fts.fill("pfcand_phirel", reco::deltaPhi(candP4, jet)); + fts.fill("pfcand_etarel", etasign * (candP4.eta() - jet.eta())); + fts.fill("pfcand_deltaR", reco::deltaR(candP4, jet)); + fts.fill("pfcand_abseta", std::abs(candP4.eta())); + fts.fill("pfcand_ptrel_log", std::log(candP4.pt() / jet.pt())); + fts.fill("pfcand_ptrel", candP4.pt() / jet.pt()); + fts.fill("pfcand_erel_log", std::log(candP4.energy() / jet.energy())); + fts.fill("pfcand_erel", candP4.energy() / jet.energy()); + fts.fill("pfcand_pt_log", std::log(candP4.pt())); + fts.fill("pfcand_mask", 1); + fts.fill("pfcand_pt_log_nopuppi", std::log(cand->pt())); + fts.fill("pfcand_e_log_nopuppi", std::log(cand->energy())); + fts.fill("pfcand_normchi2", (*normchi2_value_map_)[cand]); + + float trk_px = (*trkPt_value_map_)[cand] * std::cos((*trkPhi_value_map_)[cand]); + float trk_py = (*trkPt_value_map_)[cand] * std::sin((*trkPhi_value_map_)[cand]); + float trk_pz = (*trkPt_value_map_)[cand] * std::sinh((*trkEta_value_map_)[cand]); + math::XYZVector track_mom(trk_px, trk_py, trk_pz); + TVector3 track_direction(trk_px, trk_py, trk_pz); + double track_mag = sqrt(trk_px * trk_px + trk_py * trk_py + trk_pz * trk_pz); + fts.fill("pfcand_btagEtaRel", reco::btau::etaRel(jet_dir, track_mom)); + fts.fill("pfcand_btagPtRatio", track_direction.Perp(jet_direction) / track_mag); + fts.fill("pfcand_btagPParRatio", jet_dir.Dot(track_mom) / track_mag); - } - // in the case of reco candidate - else if (reco_cand) { - fts.fill("pfcand_hcalFrac", reco_cand->hcalEnergy() / (reco_cand->ecalEnergy() + reco_cand->hcalEnergy())); - fts.fill("pfcand_VTX_ass", vtx_ass); - fts.fill("pfcand_lostInnerHits", useTrackProperties(reco_cand) ? lost_inner_hits_from_pfcand(*reco_cand) : 0); - fts.fill("pfcand_quality", useTrackProperties(reco_cand) ? quality_from_pfcand(*reco_cand) : 0); - fts.fill("pfcand_charge", reco_cand->charge()); - fts.fill("pfcand_isEl", std::abs(reco_cand->pdgId()) == 11); - fts.fill("pfcand_isMu", std::abs(reco_cand->pdgId()) == 13); - fts.fill("pfcand_isChargedHad", std::abs(reco_cand->pdgId()) == 211); - fts.fill("pfcand_isGamma", std::abs(reco_cand->pdgId()) == 22); - fts.fill("pfcand_isNeutralHad", std::abs(reco_cand->pdgId()) == 130); - fts.fill("pfcand_dz", track ? ip_sign * track->dz(pv_->position()) : 0); - fts.fill("pfcand_dzsig", track ? ip_sign * track->dz(pv_->position()) / track->dzError() : 0); - fts.fill("pfcand_dxy", track ? ip_sign * track->dxy(pv_->position()) : 0); - fts.fill("pfcand_dxysig", track ? ip_sign * track->dxy(pv_->position()) / track->dxyError() : 0); - } - - // generic candidate observables - fts.fill("pfcand_puppiw", puppi_wgt_cache.at(cand.key())); - fts.fill("pfcand_phirel", reco::deltaPhi(candP4, jet)); - fts.fill("pfcand_etarel", etasign * (candP4.eta() - jet.eta())); - fts.fill("pfcand_deltaR", reco::deltaR(candP4, jet)); - fts.fill("pfcand_abseta", std::abs(candP4.eta())); - - fts.fill("pfcand_ptrel_log", std::log(candP4.pt() / jet.pt())); - fts.fill("pfcand_ptrel", candP4.pt() / jet.pt()); - fts.fill("pfcand_erel_log", std::log(candP4.energy() / jet.energy())); - fts.fill("pfcand_erel", candP4.energy() / jet.energy()); - fts.fill("pfcand_pt_log", std::log(candP4.pt())); - - fts.fill("pfcand_mask", 1); - fts.fill("pfcand_pt_log_nopuppi", std::log(cand->pt())); - fts.fill("pfcand_e_log_nopuppi", std::log(cand->energy())); - - float drminpfcandsv = btagbtvdeep::mindrsvpfcand(*svs_, &(*cand), std::numeric_limits::infinity()); - fts.fill("pfcand_drminsv", drminpfcandsv); - - if (track) { - auto cov = [&](unsigned i, unsigned j) { return track->covariance(i, j); }; - fts.fill("pfcand_dptdpt", cov(0, 0)); - fts.fill("pfcand_detadeta", cov(1, 1)); - fts.fill("pfcand_dphidphi", cov(2, 2)); - fts.fill("pfcand_dxydxy", cov(3, 3)); - fts.fill("pfcand_dzdz", cov(4, 4)); - fts.fill("pfcand_dxydz", cov(3, 4)); - fts.fill("pfcand_dphidxy", cov(2, 3)); - fts.fill("pfcand_dlambdadz", cov(1, 4)); - - fts.fill("pfcand_normchi2", std::floor(track->normalizedChi2())); - - trackinfo.buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_); - fts.fill("pfcand_btagEtaRel", trackinfo.getTrackEtaRel()); - fts.fill("pfcand_btagPtRatio", trackinfo.getTrackPtRatio()); - fts.fill("pfcand_btagPParRatio", trackinfo.getTrackPParRatio()); - fts.fill("pfcand_btagSip2dVal", ip_sign * trackinfo.getTrackSip2dVal()); - fts.fill("pfcand_btagSip2dSig", ip_sign * trackinfo.getTrackSip2dSig()); - fts.fill("pfcand_btagSip3dVal", ip_sign * trackinfo.getTrackSip3dVal()); - fts.fill("pfcand_btagSip3dSig", ip_sign * trackinfo.getTrackSip3dSig()); - fts.fill("pfcand_btagJetDistVal", trackinfo.getTrackJetDistVal()); - } else { - fts.fill("pfcand_normchi2", 999); - fts.fill("pfcand_dptdpt", 0); - fts.fill("pfcand_detadeta", 0); - fts.fill("pfcand_dphidphi", 0); - fts.fill("pfcand_dxydxy", 0); - fts.fill("pfcand_dzdz", 0); - fts.fill("pfcand_dxydz", 0); - fts.fill("pfcand_dphidxy", 0); - fts.fill("pfcand_dlambdadz", 0); - fts.fill("pfcand_btagEtaRel", 0); - fts.fill("pfcand_btagPtRatio", 0); - fts.fill("pfcand_btagPParRatio", 0); - fts.fill("pfcand_btagSip2dVal", 0); - fts.fill("pfcand_btagSip2dSig", 0); - fts.fill("pfcand_btagSip3dVal", 0); - fts.fill("pfcand_btagSip3dSig", 0); - fts.fill("pfcand_btagJetDistVal", 0); - } - - // subjets only if the incomming jets is a PAT one - const auto *patJet = dynamic_cast(&jet); - if (patJet and patJet->nSubjetCollections() > 0) { - auto subjets = patJet->subjets(); - std::sort(subjets.begin(), subjets.end(), [](const edm::Ptr &p1, const edm::Ptr &p2) { - return p1->pt() > p2->pt(); - }); - fts.fill("pfcand_drsubjet1", !subjets.empty() ? reco::deltaR(*cand, *subjets.at(0)) : -1); - fts.fill("pfcand_drsubjet2", subjets.size() > 1 ? reco::deltaR(*cand, *subjets.at(1)) : -1); - } else { - fts.fill("pfcand_drsubjet1", -1); - fts.fill("pfcand_drsubjet2", -1); - } - } - // using HLT features - else { - pat::PackedCandidate candidate; - math::XYZPoint pv_ass_pos; - // In case input is a packed candidate (evaluate HLT network on offline) - if (packed_cand) { - candidate = *packed_cand; - pv_ass = reco::VertexRef(vtxs_, 0); - pv_ass_pos = pv_ass->position(); - } - // In case input is a reco::PFCandidate - else if (reco_cand) { - // follow what is done in PhysicsTools/PatAlgos/plugins/PATPackedCandidateProducer.cc to minimize differences between HLT and RECO choices of input observables - // no reference to vertex, take closest in dz or position 0 - if (not pv_ass.isNonnull()) { - if (track) { - float z_dist = 99999; - int pv_pos = -1; - for (size_t iv = 0; iv < vtxs_->size(); iv++) { - float dz = std::abs(track->dz(((*vtxs_)[iv]).position())); - if (dz < z_dist) { - z_dist = dz; - pv_pos = iv; - } - } - pv_ass = reco::VertexRef(vtxs_, pv_pos); - } else - pv_ass = reco::VertexRef(vtxs_, 0); - } - pv_ass_pos = pv_ass->position(); - - // create a transient packed candidate - if (track) { - candidate = pat::PackedCandidate(cand->polarP4(), - track->referencePoint(), - track->pt(), - track->eta(), - track->phi(), - cand->pdgId(), - PVRefProd, - pv_ass.key()); - candidate.setAssociationQuality(pat::PackedCandidate::PVAssociationQuality( - btagbtvdeep::vtx_ass_from_pfcand(*reco_cand, pv_ass_quality, pv_ass))); - candidate.setCovarianceVersion(0); - pat::PackedCandidate::LostInnerHits lostHits = pat::PackedCandidate::noLostInnerHits; - int nlost = track->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS); - if (nlost == 0) { - if (track->hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) - lostHits = pat::PackedCandidate::validHitInFirstPixelBarrelLayer; - } else - lostHits = (nlost == 1 ? pat::PackedCandidate::oneLostInnerHit : pat::PackedCandidate::moreLostInnerHits); - candidate.setLostInnerHits(lostHits); - - if (useTrackProperties(reco_cand) or - std::find(whiteListSV.begin(), whiteListSV.end(), icand) != whiteListSV.end() or - std::find(whiteListTk.begin(), whiteListTk.end(), reco_cand->trackRef()) != whiteListTk.end()) { - candidate.setFirstHit(track->hitPattern().getHitPattern(reco::HitPattern::TRACK_HITS, 0)); - if (abs(cand->pdgId()) == 22) - candidate.setTrackProperties(*track, 0, 0); - else { - if (track->hitPattern().numberOfValidPixelHits() > min_valid_pixel_hits_) - candidate.setTrackProperties(*track, 8, 0); - else - candidate.setTrackProperties(*track, 264, 0); - } - } else { - if (candidate.pt() > min_track_pt_property_) { - if (track->hitPattern().numberOfValidPixelHits() > 0) - candidate.setTrackProperties(*track, 520, 0); - else - candidate.setTrackProperties(*track, 776, 0); - } - } - candidate.setTrackHighPurity(reco_cand->trackRef().isNonnull() and - reco_cand->trackRef()->quality(reco::Track::highPurity)); - } else { - candidate = pat::PackedCandidate( - cand->polarP4(), pv_ass_pos, cand->pt(), cand->eta(), cand->phi(), cand->pdgId(), PVRefProd, pv_ass.key()); - candidate.setAssociationQuality( - pat::PackedCandidate::PVAssociationQuality(pat::PackedCandidate::UsedInFitTight)); - } - /// override track - track = candidate.bestTrack(); - } - - TVector3 cand_direction(candP3.x(), candP3.y(), candP3.z()); - - fts.fill("jet_pfcand_pt_log", std::log(candP4.pt())); - fts.fill("jet_pfcand_energy_log", std::log(candP4.energy())); - fts.fill("jet_pfcand_eta", candP4.eta()); - fts.fill("jet_pfcand_deta", jet_direction.Eta() - cand_direction.Eta()); - fts.fill("jet_pfcand_dphi", jet_direction.DeltaPhi(cand_direction)); - fts.fill("jet_pfcand_charge", cand->charge()); - fts.fill("jet_pfcand_etarel", reco::btau::etaRel(jet_dir, candP3)); - fts.fill("jet_pfcand_pperp_ratio", jet_direction.Perp(cand_direction) / cand_direction.Mag()); - fts.fill("jet_pfcand_ppara_ratio", jet_direction.Dot(cand_direction) / cand_direction.Mag()); - fts.fill("jet_pfcand_frompv", candidate.fromPV()); - fts.fill("jet_pfcand_dz", candidate.dz(pv_ass_pos)); - fts.fill("jet_pfcand_dxy", candidate.dxy(pv_ass_pos)); - fts.fill("jet_pfcand_puppiw", puppi_wgt_cache.at(cand.key())); - fts.fill("jet_pfcand_nlostinnerhits", candidate.lostInnerHits()); - fts.fill("jet_pfcand_nhits", candidate.numberOfHits()); - fts.fill("jet_pfcand_npixhits", candidate.numberOfPixelHits()); - fts.fill("jet_pfcand_nstriphits", candidate.stripLayersWithMeasurement()); - fts.fill("pfcand_mask", 1); - - if (track) { - fts.fill("jet_pfcand_dzsig", fabs(candidate.dz(pv_ass_pos)) / candidate.dzError()); - fts.fill("jet_pfcand_dxysig", fabs(candidate.dxy(pv_ass_pos)) / candidate.dxyError()); - fts.fill("jet_pfcand_track_chi2", track->normalizedChi2()); - fts.fill("jet_pfcand_track_qual", track->qualityMask()); - - reco::TransientTrack transientTrack = track_builder_->build(*track); - Measurement1D meas_ip2d = - IPTools::signedTransverseImpactParameter(transientTrack, jet_ref_track_dir, *pv_).second; - Measurement1D meas_ip3d = IPTools::signedImpactParameter3D(transientTrack, jet_ref_track_dir, *pv_).second; - Measurement1D meas_jetdist = IPTools::jetTrackDistance(transientTrack, jet_ref_track_dir, *pv_).second; - Measurement1D meas_decayl = IPTools::signedDecayLength3D(transientTrack, jet_ref_track_dir, *pv_).second; - - fts.fill("jet_pfcand_trackjet_d3d", meas_ip3d.value()); - fts.fill("jet_pfcand_trackjet_d3dsig", fabs(meas_ip3d.significance())); - fts.fill("jet_pfcand_trackjet_dist", -meas_jetdist.value()); - fts.fill("jet_pfcand_trackjet_decayL", meas_decayl.value()); - } else { - fts.fill("jet_pfcand_dzsig", 0); - fts.fill("jet_pfcand_dxysig", 0); - fts.fill("jet_pfcand_track_chi2", 0); - fts.fill("jet_pfcand_track_qual", 0); - fts.fill("jet_pfcand_trackjet_d3d", 0); - fts.fill("jet_pfcand_trackjet_d3dsig", 0); - fts.fill("jet_pfcand_trackjet_dist", 0); - fts.fill("jet_pfcand_trackjet_decayL", 0); - } + } else { + // candidate track + const reco::Track *track = nullptr; + if (packed_cand) + track = packed_cand->bestTrack(); + else if (reco_cand and useTrackProperties(reco_cand)) + track = reco_cand->bestTrack(); + + // reco-vertex association + int pv_ass_quality = 0; + reco::VertexRef pv_ass; + float vtx_ass = 0; + if (reco_cand) { + if (use_pvasq_value_map_) { + pv_ass_quality = (*pvasq_value_map_)[cand]; + pv_ass = (*pvas_)[cand]; + vtx_ass = vtx_ass_from_pfcand(*reco_cand, pv_ass_quality, pv_ass); + } else + throw edm::Exception(edm::errors::InvalidReference) << "Vertex association missing"; + } + + // Building offline features + if (not use_hlt_features_) { + // in case of packed candidate + if (packed_cand) { + float hcal_fraction = 0.; + if (packed_cand->pdgId() == 1 or packed_cand->pdgId() == 130) + hcal_fraction = packed_cand->hcalFraction(); + else if (packed_cand->isIsolatedChargedHadron()) + hcal_fraction = packed_cand->rawHcalFraction(); + + fts.fill("pfcand_hcalFrac", hcal_fraction); + fts.fill("pfcand_VTX_ass", packed_cand->pvAssociationQuality()); + fts.fill("pfcand_lostInnerHits", packed_cand->lostInnerHits()); + fts.fill("pfcand_quality", track ? track->qualityMask() : 0); + fts.fill("pfcand_charge", packed_cand->charge()); + fts.fill("pfcand_isEl", std::abs(packed_cand->pdgId()) == 11); + fts.fill("pfcand_isMu", std::abs(packed_cand->pdgId()) == 13); + fts.fill("pfcand_isChargedHad", std::abs(packed_cand->pdgId()) == 211); + fts.fill("pfcand_isGamma", std::abs(packed_cand->pdgId()) == 22); + fts.fill("pfcand_isNeutralHad", std::abs(packed_cand->pdgId()) == 130); + fts.fill("pfcand_dz", ip_sign * packed_cand->dz()); + fts.fill("pfcand_dxy", ip_sign * packed_cand->dxy()); + fts.fill("pfcand_dzsig", track ? ip_sign * packed_cand->dz() / packed_cand->dzError() : 0); + fts.fill("pfcand_dxysig", track ? ip_sign * packed_cand->dxy() / packed_cand->dxyError() : 0); + + } + // in the case of reco candidate + else if (reco_cand) { + fts.fill("pfcand_hcalFrac", reco_cand->hcalEnergy() / (reco_cand->ecalEnergy() + reco_cand->hcalEnergy())); + fts.fill("pfcand_VTX_ass", vtx_ass); + fts.fill("pfcand_lostInnerHits", useTrackProperties(reco_cand) ? lost_inner_hits_from_pfcand(*reco_cand) : 0); + fts.fill("pfcand_quality", useTrackProperties(reco_cand) ? quality_from_pfcand(*reco_cand) : 0); + fts.fill("pfcand_charge", reco_cand->charge()); + fts.fill("pfcand_isEl", std::abs(reco_cand->pdgId()) == 11); + fts.fill("pfcand_isMu", std::abs(reco_cand->pdgId()) == 13); + fts.fill("pfcand_isChargedHad", std::abs(reco_cand->pdgId()) == 211); + fts.fill("pfcand_isGamma", std::abs(reco_cand->pdgId()) == 22); + fts.fill("pfcand_isNeutralHad", std::abs(reco_cand->pdgId()) == 130); + fts.fill("pfcand_dz", track ? ip_sign * track->dz(pv_->position()) : 0); + fts.fill("pfcand_dzsig", track ? ip_sign * track->dz(pv_->position()) / track->dzError() : 0); + fts.fill("pfcand_dxy", track ? ip_sign * track->dxy(pv_->position()) : 0); + fts.fill("pfcand_dxysig", track ? ip_sign * track->dxy(pv_->position()) / track->dxyError() : 0); + } + + // generic candidate observables + fts.fill("pfcand_puppiw", puppi_wgt_cache.at(cand.key())); + fts.fill("pfcand_phirel", reco::deltaPhi(candP4, jet)); + fts.fill("pfcand_etarel", etasign * (candP4.eta() - jet.eta())); + fts.fill("pfcand_deltaR", reco::deltaR(candP4, jet)); + fts.fill("pfcand_abseta", std::abs(candP4.eta())); + + fts.fill("pfcand_ptrel_log", std::log(candP4.pt() / jet.pt())); + fts.fill("pfcand_ptrel", candP4.pt() / jet.pt()); + fts.fill("pfcand_erel_log", std::log(candP4.energy() / jet.energy())); + fts.fill("pfcand_erel", candP4.energy() / jet.energy()); + fts.fill("pfcand_pt_log", std::log(candP4.pt())); + + fts.fill("pfcand_mask", 1); + fts.fill("pfcand_pt_log_nopuppi", std::log(cand->pt())); + fts.fill("pfcand_e_log_nopuppi", std::log(cand->energy())); + + float drminpfcandsv = btagbtvdeep::mindrsvpfcand(*svs_, &(*cand), std::numeric_limits::infinity()); + fts.fill("pfcand_drminsv", drminpfcandsv); + + if (track) { + auto cov = [&](unsigned i, unsigned j) { return track->covariance(i, j); }; + fts.fill("pfcand_dptdpt", cov(0, 0)); + fts.fill("pfcand_detadeta", cov(1, 1)); + fts.fill("pfcand_dphidphi", cov(2, 2)); + fts.fill("pfcand_dxydxy", cov(3, 3)); + fts.fill("pfcand_dzdz", cov(4, 4)); + fts.fill("pfcand_dxydz", cov(3, 4)); + fts.fill("pfcand_dphidxy", cov(2, 3)); + fts.fill("pfcand_dlambdadz", cov(1, 4)); + + fts.fill("pfcand_normchi2", std::floor(track->normalizedChi2())); + + (*trackinfo).buildTrackInfo(&(*cand), jet_dir, jet_ref_track_dir, *pv_); + fts.fill("pfcand_btagEtaRel", (*trackinfo).getTrackEtaRel()); + fts.fill("pfcand_btagPtRatio", (*trackinfo).getTrackPtRatio()); + fts.fill("pfcand_btagPParRatio", (*trackinfo).getTrackPParRatio()); + fts.fill("pfcand_btagSip2dVal", ip_sign * (*trackinfo).getTrackSip2dVal()); + fts.fill("pfcand_btagSip2dSig", ip_sign * (*trackinfo).getTrackSip2dSig()); + fts.fill("pfcand_btagSip3dVal", ip_sign * (*trackinfo).getTrackSip3dVal()); + fts.fill("pfcand_btagSip3dSig", ip_sign * (*trackinfo).getTrackSip3dSig()); + fts.fill("pfcand_btagJetDistVal", (*trackinfo).getTrackJetDistVal()); + } else { + fts.fill("pfcand_normchi2", 999); + fts.fill("pfcand_dptdpt", 0); + fts.fill("pfcand_detadeta", 0); + fts.fill("pfcand_dphidphi", 0); + fts.fill("pfcand_dxydxy", 0); + fts.fill("pfcand_dzdz", 0); + fts.fill("pfcand_dxydz", 0); + fts.fill("pfcand_dphidxy", 0); + fts.fill("pfcand_dlambdadz", 0); + fts.fill("pfcand_btagEtaRel", 0); + fts.fill("pfcand_btagPtRatio", 0); + fts.fill("pfcand_btagPParRatio", 0); + fts.fill("pfcand_btagSip2dVal", 0); + fts.fill("pfcand_btagSip2dSig", 0); + fts.fill("pfcand_btagSip3dVal", 0); + fts.fill("pfcand_btagSip3dSig", 0); + fts.fill("pfcand_btagJetDistVal", 0); + } + + // subjets only if the incomming jets is a PAT one + const auto *patJet = dynamic_cast(&jet); + if (patJet and patJet->nSubjetCollections() > 0) { + auto subjets = patJet->subjets(); + std::sort(subjets.begin(), subjets.end(), [](const edm::Ptr &p1, const edm::Ptr &p2) { + return p1->pt() > p2->pt(); + }); + fts.fill("pfcand_drsubjet1", !subjets.empty() ? reco::deltaR(*cand, *subjets.at(0)) : -1); + fts.fill("pfcand_drsubjet2", subjets.size() > 1 ? reco::deltaR(*cand, *subjets.at(1)) : -1); + } else { + fts.fill("pfcand_drsubjet1", -1); + fts.fill("pfcand_drsubjet2", -1); + } + } + // using HLT features + else { + pat::PackedCandidate candidate; + math::XYZPoint pv_ass_pos; + // In case input is a packed candidate (evaluate HLT network on offline) + if (packed_cand) { + candidate = *packed_cand; + pv_ass = reco::VertexRef(vtxs_, 0); + pv_ass_pos = pv_ass->position(); + } + // In case input is a reco::PFCandidate + else if (reco_cand) { + // follow what is done in PhysicsTools/PatAlgos/plugins/PATPackedCandidateProducer.cc to minimize differences between HLT and RECO choices of input observables + // no reference to vertex, take closest in dz or position 0 + if (not pv_ass.isNonnull()) { + if (track) { + float z_dist = 99999; + int pv_pos = -1; + for (size_t iv = 0; iv < vtxs_->size(); iv++) { + float dz = std::abs(track->dz(((*vtxs_)[iv]).position())); + if (dz < z_dist) { + z_dist = dz; + pv_pos = iv; + } + } + pv_ass = reco::VertexRef(vtxs_, pv_pos); + } else + pv_ass = reco::VertexRef(vtxs_, 0); + } + pv_ass_pos = pv_ass->position(); + + // create a transient packed candidate + if (track) { + candidate = pat::PackedCandidate(cand->polarP4(), + track->referencePoint(), + track->pt(), + track->eta(), + track->phi(), + cand->pdgId(), + (*PVRefProd), + pv_ass.key()); + candidate.setAssociationQuality(pat::PackedCandidate::PVAssociationQuality( + btagbtvdeep::vtx_ass_from_pfcand(*reco_cand, pv_ass_quality, pv_ass))); + candidate.setCovarianceVersion(0); + pat::PackedCandidate::LostInnerHits lostHits = pat::PackedCandidate::noLostInnerHits; + int nlost = track->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS); + if (nlost == 0) { + if (track->hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector::PixelBarrel, 1)) + lostHits = pat::PackedCandidate::validHitInFirstPixelBarrelLayer; + } else + lostHits = (nlost == 1 ? pat::PackedCandidate::oneLostInnerHit : pat::PackedCandidate::moreLostInnerHits); + candidate.setLostInnerHits(lostHits); + + if (useTrackProperties(reco_cand) or + std::find(whiteListSV.begin(), whiteListSV.end(), icand) != whiteListSV.end() or + std::find(whiteListTk.begin(), whiteListTk.end(), reco_cand->trackRef()) != whiteListTk.end()) { + candidate.setFirstHit(track->hitPattern().getHitPattern(reco::HitPattern::TRACK_HITS, 0)); + if (abs(cand->pdgId()) == 22) + candidate.setTrackProperties(*track, 0, 0); + else { + if (track->hitPattern().numberOfValidPixelHits() > min_valid_pixel_hits_) + candidate.setTrackProperties(*track, 8, 0); + else + candidate.setTrackProperties(*track, 264, 0); + } + } else { + if (candidate.pt() > min_track_pt_property_) { + if (track->hitPattern().numberOfValidPixelHits() > 0) + candidate.setTrackProperties(*track, 520, 0); + else + candidate.setTrackProperties(*track, 776, 0); + } + } + candidate.setTrackHighPurity(reco_cand->trackRef().isNonnull() and + reco_cand->trackRef()->quality(reco::Track::highPurity)); + } else { + candidate = pat::PackedCandidate( + cand->polarP4(), pv_ass_pos, cand->pt(), cand->eta(), cand->phi(), cand->pdgId(), (*PVRefProd), pv_ass.key()); + candidate.setAssociationQuality( + pat::PackedCandidate::PVAssociationQuality(pat::PackedCandidate::UsedInFitTight)); + } + /// override track + track = candidate.bestTrack(); + } + + TVector3 cand_direction(candP3.x(), candP3.y(), candP3.z()); + + fts.fill("jet_pfcand_pt_log", std::log(candP4.pt())); + fts.fill("jet_pfcand_energy_log", std::log(candP4.energy())); + fts.fill("jet_pfcand_eta", candP4.eta()); + fts.fill("jet_pfcand_deta", jet_direction.Eta() - cand_direction.Eta()); + fts.fill("jet_pfcand_dphi", jet_direction.DeltaPhi(cand_direction)); + fts.fill("jet_pfcand_charge", cand->charge()); + fts.fill("jet_pfcand_etarel", reco::btau::etaRel(jet_dir, candP3)); + fts.fill("jet_pfcand_pperp_ratio", jet_direction.Perp(cand_direction) / cand_direction.Mag()); + fts.fill("jet_pfcand_ppara_ratio", jet_direction.Dot(cand_direction) / cand_direction.Mag()); + fts.fill("jet_pfcand_frompv", candidate.fromPV()); + fts.fill("jet_pfcand_dz", candidate.dz(pv_ass_pos)); + fts.fill("jet_pfcand_dxy", candidate.dxy(pv_ass_pos)); + fts.fill("jet_pfcand_puppiw", puppi_wgt_cache.at(cand.key())); + fts.fill("jet_pfcand_nlostinnerhits", candidate.lostInnerHits()); + fts.fill("jet_pfcand_nhits", candidate.numberOfHits()); + fts.fill("jet_pfcand_npixhits", candidate.numberOfPixelHits()); + fts.fill("jet_pfcand_nstriphits", candidate.stripLayersWithMeasurement()); + fts.fill("pfcand_mask", 1); + + if (track) { + fts.fill("jet_pfcand_dzsig", fabs(candidate.dz(pv_ass_pos)) / candidate.dzError()); + fts.fill("jet_pfcand_dxysig", fabs(candidate.dxy(pv_ass_pos)) / candidate.dxyError()); + fts.fill("jet_pfcand_track_chi2", track->normalizedChi2()); + fts.fill("jet_pfcand_track_qual", track->qualityMask()); + + reco::TransientTrack transientTrack = track_builder_->build(*track); + Measurement1D meas_ip2d = + IPTools::signedTransverseImpactParameter(transientTrack, jet_ref_track_dir, *pv_).second; + Measurement1D meas_ip3d = IPTools::signedImpactParameter3D(transientTrack, jet_ref_track_dir, *pv_).second; + Measurement1D meas_jetdist = IPTools::jetTrackDistance(transientTrack, jet_ref_track_dir, *pv_).second; + Measurement1D meas_decayl = IPTools::signedDecayLength3D(transientTrack, jet_ref_track_dir, *pv_).second; + + fts.fill("jet_pfcand_trackjet_d3d", meas_ip3d.value()); + fts.fill("jet_pfcand_trackjet_d3dsig", fabs(meas_ip3d.significance())); + fts.fill("jet_pfcand_trackjet_dist", -meas_jetdist.value()); + fts.fill("jet_pfcand_trackjet_decayL", meas_decayl.value()); + } else { + fts.fill("jet_pfcand_dzsig", 0); + fts.fill("jet_pfcand_dxysig", 0); + fts.fill("jet_pfcand_track_chi2", 0); + fts.fill("jet_pfcand_track_qual", 0); + fts.fill("jet_pfcand_trackjet_d3d", 0); + fts.fill("jet_pfcand_trackjet_d3dsig", 0); + fts.fill("jet_pfcand_trackjet_dist", 0); + fts.fill("jet_pfcand_trackjet_decayL", 0); + } + } } icand++; } diff --git a/RecoBTag/ONNXRuntime/plugins/BoostedJetONNXValueMapProducer.cc b/RecoBTag/ONNXRuntime/plugins/BoostedJetONNXValueMapProducer.cc new file mode 100644 index 0000000000000..53bce37911dce --- /dev/null +++ b/RecoBTag/ONNXRuntime/plugins/BoostedJetONNXValueMapProducer.cc @@ -0,0 +1,237 @@ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/Framework/interface/makeRefToBaseProdFrom.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/BTauReco/interface/JetTag.h" + +#include "DataFormats/BTauReco/interface/DeepBoostedJetTagInfo.h" + +#include "PhysicsTools/ONNXRuntime/interface/ONNXRuntime.h" + +#include "RecoBTag/FeatureTools/interface/deep_helpers.h" + +#include +#include +#include +#include +#include + +using namespace cms::Ort; +using namespace btagbtvdeep; + +class BoostedJetONNXValueMapProducer : public edm::stream::EDProducer> { +public: + explicit BoostedJetONNXValueMapProducer(const edm::ParameterSet &, const ONNXRuntime *); + ~BoostedJetONNXValueMapProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions &); + + static std::unique_ptr initializeGlobalCache(const edm::ParameterSet &); + static void globalEndJob(const ONNXRuntime *); + +private: + typedef std::vector TagInfoCollection; + typedef reco::JetTagCollection JetTagCollection; + + void produce(edm::Event &, const edm::EventSetup &) override; + + void make_inputs(const reco::DeepBoostedJetTagInfo &taginfo); + + const edm::EDGetTokenT src_; + std::vector flav_names_; // names of the output scores + edm::EDGetTokenT> jet_token_; + std::vector input_names_; // names of each input group - the ordering is important! + std::vector> input_shapes_; // shapes of each input group (-1 for dynamic axis) + std::vector input_sizes_; // total length of each input vector + std::unordered_map prep_info_map_; // preprocessing info for each input group + + FloatArrays data_; + + bool debug_ = false; +}; + +BoostedJetONNXValueMapProducer::BoostedJetONNXValueMapProducer(const edm::ParameterSet &iConfig, const ONNXRuntime *cache) + : src_(consumes(iConfig.getParameter("src"))), + flav_names_(iConfig.getParameter>("flav_names")), + jet_token_(consumes>(iConfig.getParameter("jets"))), + debug_(iConfig.getUntrackedParameter("debugMode", false)) { + ParticleNetConstructor(iConfig, true, input_names_, prep_info_map_, input_shapes_, input_sizes_, &data_); + + if (debug_) { + LogDebug("BoostedJetONNXValueMapProducer") << ":" << std::endl; + for (unsigned i = 0; i < input_names_.size(); ++i) { + const auto &group_name = input_names_.at(i); + if (!input_shapes_.empty()) { + LogDebug("BoostedJetONNXValueMapProducer") << group_name << "\nshapes: "; + for (const auto &x : input_shapes_.at(i)) { + LogDebug("BoostedJetONNXValueMapProducer") << x << ", "; + } + } + LogDebug("BoostedJetONNXValueMapProducer") << "\nvariables: "; + for (const auto &x : prep_info_map_.at(group_name).var_names) { + LogDebug("BoostedJetONNXValueMapProducer") << x << ", "; + } + LogDebug("BoostedJetONNXValueMapProducer") << "\n"; + } + LogDebug("BoostedJetONNXValueMapProducer") << "flav_names: "; + for (const auto &flav_name : flav_names_) { + LogDebug("BoostedJetONNXValueMapProducer") << flav_name << ", "; + } + LogDebug("BoostedJetONNXValueMapProducer") << "\n"; + } + + // get output names from flav_names + for (const auto &flav_name : flav_names_) { + produces>(flav_name); + } +} + +BoostedJetONNXValueMapProducer::~BoostedJetONNXValueMapProducer() {} + +void BoostedJetONNXValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + // pfDeepBoostedJetTags + edm::ParameterSetDescription desc; + desc.add("src", edm::InputTag("pfDeepBoostedJetTagInfos")); + desc.add("preprocess_json", ""); + desc.add("jets", edm::InputTag("ak8PFJetsPuppi")); + // `preprocessParams` is deprecated -- use the preprocessing json file instead + edm::ParameterSetDescription preprocessParams; + preprocessParams.setAllowAnything(); + preprocessParams.setComment("`preprocessParams` is deprecated, please use `preprocess_json` instead."); + desc.addOptional("preprocessParams", preprocessParams); + desc.add("model_path", + edm::FileInPath("RecoBTag/Combined/data/DeepBoostedJet/V02/full/resnet.onnx")); + desc.add>("flav_names", + std::vector{ + "probTbcq", + "probTbqq", + "probTbc", + "probTbq", + "probWcq", + "probWqq", + "probZbb", + "probZcc", + "probZqq", + "probHbb", + "probHcc", + "probHqqqq", + "probQCDbb", + "probQCDcc", + "probQCDb", + "probQCDc", + "probQCDothers", + }); + desc.addOptionalUntracked("debugMode", false); + + descriptions.addWithDefaultLabel(desc); +} + +std::unique_ptr BoostedJetONNXValueMapProducer::initializeGlobalCache(const edm::ParameterSet &iConfig) { + return std::make_unique(iConfig.getParameter("model_path").fullPath()); +} + +void BoostedJetONNXValueMapProducer::globalEndJob(const ONNXRuntime *cache) {} + +void BoostedJetONNXValueMapProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { + edm::Handle tag_infos; + iEvent.getByToken(src_, tag_infos); + auto jets = iEvent.getHandle(jet_token_); + + // initialize output collection + std::vector> output_scores(flav_names_.size(), std::vector(tag_infos->size(), -1.0)); + + for (unsigned jet_n = 0; jet_n < tag_infos->size(); ++jet_n) { + const auto &taginfo = (*tag_infos)[jet_n]; + std::vector outputs(flav_names_.size(), 0); // init as all zeros + + if (!taginfo.features().empty()) { + // convert inputs + make_inputs(taginfo); + // run prediction and get outputs + outputs = globalCache()->run(input_names_, data_, input_shapes_)[0]; + assert(outputs.size() == flav_names_.size()); + } + + for (std::size_t flav_n = 0; flav_n < flav_names_.size(); flav_n++) { + output_scores[flav_n][jet_n] = outputs[flav_n]; + } + } + + if (debug_) { + LogDebug("produce") << ":" << std::endl + << "=== " << iEvent.id().run() << ":" << iEvent.id().luminosityBlock() << ":" + << iEvent.id().event() << " ===" << std::endl; + for (unsigned jet_n = 0; jet_n < tag_infos->size(); ++jet_n) { + const auto &jet_ref = tag_infos->at(jet_n).jet(); + LogDebug("produce") << " - Jet #" << jet_n << ", pt=" << jet_ref->pt() << ", eta=" << jet_ref->eta() + << ", phi=" << jet_ref->phi() << std::endl; + for (std::size_t flav_n = 0; flav_n < flav_names_.size(); ++flav_n) { + LogDebug("produce") << " " << flav_names_.at(flav_n) << " = " << output_scores[flav_n][jet_n] + << std::endl; + } + } + } + + // put into the event + for (size_t k = 0; k < output_scores.size(); k++) { + std::unique_ptr> VM(new edm::ValueMap()); + edm::ValueMap::Filler filler(*VM); + filler.insert(jets, output_scores.at(k).begin(), output_scores.at(k).end()); + filler.fill(); + iEvent.put(std::move(VM), flav_names_[k]); + } +} + +void BoostedJetONNXValueMapProducer::make_inputs(const reco::DeepBoostedJetTagInfo &taginfo) { + for (unsigned igroup = 0; igroup < input_names_.size(); ++igroup) { + const auto &group_name = input_names_[igroup]; + const auto &prep_params = prep_info_map_.at(group_name); + auto &group_values = data_[igroup]; + group_values.resize(input_sizes_[igroup]); + // first reset group_values to 0 + std::fill(group_values.begin(), group_values.end(), 0); + unsigned curr_pos = 0; + // transform/pad + for (unsigned i = 0; i < prep_params.var_names.size(); ++i) { + const auto &varname = prep_params.var_names[i]; + const auto &raw_value = taginfo.features().get(varname); + const auto &info = prep_params.info(varname); + int insize = center_norm_pad(raw_value, + info.center, + info.norm_factor, + prep_params.min_length, + prep_params.max_length, + group_values, + curr_pos, + info.pad, + info.replace_inf_value, + info.lower_bound, + info.upper_bound); + curr_pos += insize; + if (i == 0 && (!input_shapes_.empty())) { + input_shapes_[igroup][2] = insize; + } + + if (debug_) { + LogDebug("make_inputs") << ":" << std::endl + << " -- var=" << varname << ", center=" << info.center << ", scale=" << info.norm_factor + << ", replace=" << info.replace_inf_value << ", pad=" << info.pad << std::endl; + for (unsigned i = curr_pos - insize; i < curr_pos; i++) { + LogDebug("make_inputs") << group_values[i] << ","; + } + LogDebug("make_inputs") << std::endl; + } + } + group_values.resize(curr_pos); + } +} + +//define this as a plug-in +DEFINE_FWK_MODULE(BoostedJetONNXValueMapProducer);