diff --git a/RecoJets/JetProducers/plugins/BetaStarVarProducer.cc b/RecoJets/JetProducers/plugins/BetaStarVarProducer.cc new file mode 100644 index 0000000000000..fe32971df557e --- /dev/null +++ b/RecoJets/JetProducers/plugins/BetaStarVarProducer.cc @@ -0,0 +1,171 @@ +// -*- C++ -*- +// +// Package: PhysicsTools/NanoAOD +// Class: BetaStarVarProducer +// +/**\class BetaStarVarProducer BetaStarVarProducer.cc PhysicsTools/NanoAOD/plugins/BetaStarVarProducer.cc + + Description: This produces value maps to store CHS-related variables for JERC. + This includes the charged hadrons associated to CHS jets, + and those associated to PU that are within the CHS jet. + + Implementation: + This uses a ValueMap producer functionality, and + loops over the input candidates (usually PackedCandidates) + that are associated to each jet, counting the candidates associated + to the PV and those not. +*/ +// +// Original Author: Sal Rappoccio +// +// + + +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" +#include "DataFormats/VertexReco/interface/Vertex.h" + +#include "DataFormats/Common/interface/View.h" + +#include "PhysicsTools/NanoAOD/interface/MatchingUtils.h" + +template +class BetaStarVarProducer : public edm::global::EDProducer<> { + public: + explicit BetaStarVarProducer(const edm::ParameterSet &iConfig): + srcJet_(consumes>(iConfig.getParameter("srcJet"))), + srcPF_(consumes>(iConfig.getParameter("srcPF"))), + maxDR_( iConfig.getParameter("maxDR") ) + { + produces>("chargedHadronPUEnergyFraction"); + produces>("chargedHadronCHSEnergyFraction"); + } + ~BetaStarVarProducer() override {}; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + private: + void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; + + std::tuple calculateCHSEnergies( edm::Ptr const & jet, edm::View const & cands) const; + + edm::EDGetTokenT> srcJet_; + edm::EDGetTokenT> srcPF_; + double maxDR_; +}; + +template +void +BetaStarVarProducer::produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const +{ + + edm::Handle> srcJet; + iEvent.getByToken(srcJet_, srcJet); + edm::Handle> srcPF; + iEvent.getByToken(srcPF_, srcPF); + + unsigned int nJet = srcJet->size(); + + std::vector chargedHadronPUEnergyFraction(nJet,-1); + std::vector chargedHadronCHSEnergyFraction(nJet,-1); + + for ( unsigned int ij = 0; ij < nJet; ++ij ) { + auto jet = srcJet->ptrAt(ij); + std::tuple vals = calculateCHSEnergies( jet, *srcPF ); + auto chpuf = std::get<0>(vals); + auto chef = std::get<1>(vals); + chargedHadronPUEnergyFraction[ij] = chpuf; + chargedHadronCHSEnergyFraction[ij] = chef; + } + + std::unique_ptr> chargedHadronPUEnergyFractionV(new edm::ValueMap()); + edm::ValueMap::Filler fillerPU(*chargedHadronPUEnergyFractionV); + fillerPU.insert(srcJet,chargedHadronPUEnergyFraction.begin(),chargedHadronPUEnergyFraction.end()); + fillerPU.fill(); + iEvent.put(std::move(chargedHadronPUEnergyFractionV),"chargedHadronPUEnergyFraction"); + + std::unique_ptr> chargedHadronCHSEnergyFractionV(new edm::ValueMap()); + edm::ValueMap::Filler fillerCHE(*chargedHadronCHSEnergyFractionV); + fillerCHE.insert(srcJet,chargedHadronCHSEnergyFraction.begin(),chargedHadronCHSEnergyFraction.end()); + fillerCHE.fill(); + iEvent.put(std::move(chargedHadronCHSEnergyFractionV),"chargedHadronCHSEnergyFraction"); + +} + +template +std::tuple +BetaStarVarProducer::calculateCHSEnergies( edm::Ptr const & ijet, edm::View const & cands ) const { + + auto rawP4 = ijet->correctedP4(0); + std::vector jet2pu; + + // Get all of the PF candidates within a cone of dR to the jet that + // do NOT belong to the primary vertex. + // Store their indices. + for ( unsigned int icand = 0; icand < cands.size(); ++icand ) { + auto cand = cands.ptrAt(icand); + if (cand->fromPV()!=0) continue; + float dR = reco::deltaR(*ijet,*cand); + if (dR used; + auto jetConstituents = ijet->daughterPtrVector(); + for (auto const & ic : jetConstituents ) { + auto icpc = dynamic_cast(ic.get()); + if ( icpc->charge()!=0) { + che += icpc->energy(); + if (icpc->fromPV()==0) { + used.push_back( ic.key() ); + } + } + } + // Loop through the pileup PF candidates within the jet. + for (auto pidx : jet2pu) { + auto const & dtr = cands.ptrAt(pidx); + // We take the candidates that have not appeared before: these were removed by CHS + if (dtr->charge()!=0 and std::find(used.begin(),used.end(),dtr.key() )==used.end()) + pue += dtr->energy(); + } + + // Now get the fractions relative to the raw jet. + auto puf = pue / rawP4.energy(); + auto chf = che / rawP4.energy(); + + return std::tuple(puf,chf); +} + +template +void +BetaStarVarProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("srcJet")->setComment("jet input collection"); + desc.add("srcPF")->setComment("PF candidate input collection"); + desc.add("maxDR")->setComment("Maximum DR to consider for jet-->pf cand association"); + std::string modname ("BetaStarVarProducer"); + descriptions.add(modname,desc); +} + +typedef BetaStarVarProducer BetaStarPackedCandidateVarProducer; + +DEFINE_FWK_MODULE(BetaStarPackedCandidateVarProducer); + diff --git a/RecoJets/JetProducers/python/ak4PFJetsBetaStar_cfi.py b/RecoJets/JetProducers/python/ak4PFJetsBetaStar_cfi.py new file mode 100644 index 0000000000000..8898fb830d74d --- /dev/null +++ b/RecoJets/JetProducers/python/ak4PFJetsBetaStar_cfi.py @@ -0,0 +1,8 @@ +import FWCore.ParameterSet.Config as cms + + +ak4BetaStar = cms.EDProducer("BetaStarPackedCandidateVarProducer", + srcJet = cms.InputTag("slimmedJets"), + srcPF = cms.InputTag("packedPFCandidates"), + maxDR = cms.double(0.4) +)