Skip to content

Commit

Permalink
Merge pull request #31271 from lathomas/HFShowerShape_master
Browse files Browse the repository at this point in the history
Shower shape variables for HF jets
  • Loading branch information
cmsbuild authored Sep 10, 2020
2 parents 67425f0 + efc2ace commit d8780f6
Show file tree
Hide file tree
Showing 6 changed files with 290 additions and 2 deletions.
27 changes: 27 additions & 0 deletions PhysicsTools/NanoAOD/python/jets_cff.py
Original file line number Diff line number Diff line change
Expand Up @@ -237,6 +237,10 @@
puId = Var("userInt('pileupJetId:fullId')",int,doc="Pilup ID flags with 80X (2016) training"),
jetId = Var("userInt('tightId')*2+4*userInt('tightIdLepVeto')",int,doc="Jet ID flags bit1 is loose (always false in 2017 since it does not exist), bit2 is tight, bit3 is tightLepVeto"),
qgl = Var("userFloat('qgl')",float,doc="Quark vs Gluon likelihood discriminator",precision=10),
hfsigmaEtaEta = Var("userFloat('hfJetShowerShape:sigmaEtaEta')",float,doc="sigmaEtaEta for HF jets (noise discriminating variable)",precision=10),
hfsigmaPhiPhi = Var("userFloat('hfJetShowerShape:sigmaPhiPhi')",float,doc="sigmaPhiPhi for HF jets (noise discriminating variable)",precision=10),
hfcentralEtaStripSize = Var("userInt('hfJetShowerShape:centralEtaStripSize')", int, doc="eta size of the central tower strip in HF (noise discriminating variable) "),
hfadjacentEtaStripsSize = Var("userInt('hfJetShowerShape:adjacentEtaStripsSize')", int, doc="eta size of the strips next to the central tower strip in HF (noise discriminating variable) "),
nConstituents = Var("numberOfDaughters()",int,doc="Number of particles in the jet"),
rawFactor = Var("1.-jecFactor('Uncorrected')",float,doc="1 - Factor to get back to raw pT",precision=6),
chHEF = Var("chargedHadronEnergyFraction()", float, doc="charged Hadron Energy Fraction", precision= 6),
Expand Down Expand Up @@ -671,6 +675,7 @@
from RecoJets.JetProducers.QGTagger_cfi import QGTagger
qgtagger=QGTagger.clone(srcJets="updatedJets",srcVertexCollection="offlineSlimmedPrimaryVertices")


from RecoJets.JetProducers.PileupJetID_cfi import pileupJetId, _chsalgos_94x, _chsalgos_102x
pileupJetId94X=pileupJetId.clone(jets="updatedJets",algos = cms.VPSet(_chsalgos_94x),inputIsCorrected=True,applyJec=False,vertexes="offlineSlimmedPrimaryVertices")
pileupJetId102X=pileupJetId.clone(jets="updatedJets",algos = cms.VPSet(_chsalgos_102x),inputIsCorrected=True,applyJec=False,vertexes="offlineSlimmedPrimaryVertices")
Expand All @@ -685,6 +690,28 @@
for modifier in run2_miniAOD_80XLegacy, run2_nanoAOD_94X2016:
modifier.toReplaceWith(jetSequence, _jetSequence_2016)

#HF shower shape recomputation
#Only run if needed (i.e. if default MINIAOD info is missing or outdated because of new JECs...)
from RecoJets.JetProducers.hfJetShowerShape_cfi import hfJetShowerShape
hfJetShowerShapeforNanoAOD = hfJetShowerShape.clone(jets="updatedJets",vertices="offlineSlimmedPrimaryVertices")
for modifier in run2_miniAOD_80XLegacy, run2_nanoAOD_94X2016, run2_nanoAOD_94XMiniAODv1, run2_nanoAOD_94XMiniAODv2, run2_nanoAOD_102Xv1, run2_nanoAOD_106Xv1:
modifier.toModify(updatedJetsWithUserData.userFloats,
hfsigmaEtaEta = cms.InputTag('hfJetShowerShapeforNanoAOD:sigmaEtaEta'),
hfsigmaPhiPhi = cms.InputTag('hfJetShowerShapeforNanoAOD:sigmaPhiPhi'),
)
modifier.toModify(updatedJetsWithUserData.userInts,
hfcentralEtaStripSize = cms.InputTag('hfJetShowerShapeforNanoAOD:centralEtaStripSize'),
hfadjacentEtaStripsSize = cms.InputTag('hfJetShowerShapeforNanoAOD:adjacentEtaStripsSize'),
)
modifier.toModify( jetTable.variables, hfsigmaEtaEta = Var("userFloat('hfsigmaEtaEta')",float,doc="sigmaEtaEta for HF jets (noise discriminating variable)",precision=10))
modifier.toModify( jetTable.variables, hfsigmaPhiPhi = Var("userFloat('hfsigmaPhiPhi')",float,doc="sigmaPhiPhi for HF jets (noise discriminating variable)",precision=10))
modifier.toModify( jetTable.variables, hfcentralEtaStripSize = Var("userInt('hfcentralEtaStripSize')", int, doc="eta size of the central tower strip in HF (noise discriminating variable) "))
modifier.toModify( jetTable.variables, hfadjacentEtaStripsSize = Var("userInt('hfadjacentEtaStripsSize')", int, doc="eta size of the strips next to the central tower strip in HF (noise discriminating variable) "))
_jetSequence_rerunHFshowershape = jetSequence.copy()
_jetSequence_rerunHFshowershape.insert(_jetSequence_rerunHFshowershape.index(updatedJetsWithUserData), hfJetShowerShapeforNanoAOD)
modifier.toReplaceWith(jetSequence, _jetSequence_rerunHFshowershape)


#after lepton collections have been run
jetLepSequence = cms.Sequence(lepInJetVars)

Expand Down
4 changes: 4 additions & 0 deletions PhysicsTools/NanoAOD/python/nanoDQM_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -373,6 +373,10 @@
Plot1D('puId', 'puId', 8, -0.5, 7.5, 'Pilup ID flags'),
Plot1D('puIdDisc', 'puIdDisc', 20, -1, 1, 'Pilup ID discriminant with 102X (2018) training'),
Plot1D('qgl', 'qgl', 20, 0, 1, 'Quark vs Gluon likelihood discriminator'),
Plot1D('hfsigmaEtaEta', 'hfsigmaEtaEta', 20, 0, 0.2, 'sigmaEtaEta for HF jets (noise discriminating variable)'),
Plot1D('hfsigmaPhiPhi', 'hfsigmaPhiPhi', 20, 0, 0.2, 'sigmaPhiPhi for HF jets (noise discriminating variable)'),
Plot1D('hfcentralEtaStripSize', 'hfcentralEtaStripSize', 10, 0, 10, 'eta size of the central tower strip in HF (noise discriminating variable)'),
Plot1D('hfadjacentEtaStripsSize', 'hfadjacentEtaStripsSize', 10, 0, 10, 'eta size of the strips next to the central tower strip in HF (noise discriminating variable)'),
Plot1D('rawFactor', 'rawFactor', 20, -0.5, 0.5, '1 - Factor to get back to raw pT'),
)
),
Expand Down
11 changes: 9 additions & 2 deletions PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py
Original file line number Diff line number Diff line change
Expand Up @@ -239,7 +239,14 @@ def miniAOD_customizeCommon(process):
process.load('RecoJets.JetProducers.QGTagger_cfi')
task.add(process.QGTaggerTask)

process.patJets.userData.userFloats.src += [ cms.InputTag('QGTagger:qgLikelihood'), ]
process.patJets.userData.userFloats.src += [ 'QGTagger:qgLikelihood', ]

#HF jet shower shape
process.load('RecoJets.JetProducers.hfJetShowerShape_cfi')
task.add(process.hfJetShowerShape)

process.patJets.userData.userFloats.src += [ 'hfJetShowerShape:sigmaEtaEta', 'hfJetShowerShape:sigmaPhiPhi']
process.patJets.userData.userInts.src += [ 'hfJetShowerShape:centralEtaStripSize', 'hfJetShowerShape:adjacentEtaStripsSize']

## DeepCSV meta discriminators (simple arithmethic on output probabilities)
process.load('RecoBTag.Combined.deepFlavour_cff')
Expand All @@ -259,7 +266,7 @@ def miniAOD_customizeCommon(process):
valueLabels = cms.vstring('pt','emEnergyFraction'),
lazyParser = cms.bool(True) )
task.add(process.caloJetMap)
process.patJets.userData.userFloats.src += [ cms.InputTag("caloJetMap:pt"), cms.InputTag("caloJetMap:emEnergyFraction") ]
process.patJets.userData.userFloats.src += [ 'caloJetMap:pt', 'caloJetMap:emEnergyFraction' ]

#Muon object modifications
from PhysicsTools.PatAlgos.slimming.muonIsolationsPUPPI_cfi import makeInputForPUPPIIsolationMuon
Expand Down
229 changes: 229 additions & 0 deletions RecoJets/JetProducers/plugins/HFJetShowerShape.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,229 @@
// -*- C++ -*-
//
// Package: HFJetShowerShape/HFJetShowerShape
// Class: HFJetShowerShape
//
/**\class HFJetShowerShape HFJetShowerShape.cc HFJetShowerShape/HFJetShowerShape/plugins/HFJetShowerShape.cc
Description: [one line class summary]
Implementation:
[Notes on implementation]
*/
//
// Original Author: Laurent Thomas
// Created: Tue, 25 Aug 2020 09:17:42 GMT
//
//

// system include files
#include <memory>

// user include files
#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 "FWCore/Utilities/interface/StreamID.h"

#include "DataFormats/JetReco/interface/PFJet.h"
#include "DataFormats/VertexReco/interface/Vertex.h"

#include "DataFormats/Math/interface/deltaPhi.h"

#include <iostream>

//
// class declaration
//

class HFJetShowerShape : public edm::stream::EDProducer<> {
public:
explicit HFJetShowerShape(const edm::ParameterSet&);

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
void produce(edm::Event&, const edm::EventSetup&) override;
template <typename T>
void putInEvent(const std::string&, const edm::Handle<edm::View<reco::Jet>>&, std::vector<T>, edm::Event&) const;

const edm::EDGetTokenT<edm::View<reco::Jet>> jets_token_;
const edm::EDGetTokenT<std::vector<reco::Vertex>> vertices_token_;

//Jet pt/eta thresholds
const double jetPtThreshold_, jetEtaThreshold_;
//HF geometry
const double hfTowerEtaWidth_, hfTowerPhiWidth_;
//Variables for PU subtraction
const double vertexRecoEffcy_, offsetPerPU_, jetReferenceRadius_;
//Pt thresholds for showershape variable calculation
const double stripPtThreshold_, widthPtThreshold_;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
HFJetShowerShape::HFJetShowerShape(const edm::ParameterSet& iConfig)
: jets_token_(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"))),
vertices_token_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("vertices"))),
jetPtThreshold_(iConfig.getParameter<double>("jetPtThreshold")),
jetEtaThreshold_(iConfig.getParameter<double>("jetEtaThreshold")),
hfTowerEtaWidth_(iConfig.getParameter<double>("hfTowerEtaWidth")),
hfTowerPhiWidth_(iConfig.getParameter<double>("hfTowerPhiWidth")),
vertexRecoEffcy_(iConfig.getParameter<double>("vertexRecoEffcy")),
offsetPerPU_(iConfig.getParameter<double>("offsetPerPU")),
jetReferenceRadius_(iConfig.getParameter<double>("jetReferenceRadius")),
stripPtThreshold_(iConfig.getParameter<double>("stripPtThreshold")),
widthPtThreshold_(iConfig.getParameter<double>("widthPtThreshold")) {
produces<edm::ValueMap<float>>("sigmaEtaEta");
produces<edm::ValueMap<float>>("sigmaPhiPhi");
produces<edm::ValueMap<int>>("centralEtaStripSize");
produces<edm::ValueMap<int>>("adjacentEtaStripsSize");
}

//
// member functions
//

// ------------ method called to produce the data ------------
void HFJetShowerShape::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
using namespace edm;

//Jets
auto theJets = iEvent.getHandle(jets_token_);

//Vertices
int nPV = iEvent.get(vertices_token_).size();

//Products
std::vector<float> v_sigmaEtaEta, v_sigmaPhiPhi;
v_sigmaEtaEta.reserve(theJets->size());
v_sigmaPhiPhi.reserve(theJets->size());
std::vector<int> v_size_CentralEtaStrip, v_size_AdjacentEtaStrips;
v_size_CentralEtaStrip.reserve(theJets->size());
v_size_AdjacentEtaStrips.reserve(theJets->size());

//Et offset for HF PF candidates
double puoffset = offsetPerPU_ / (M_PI * jetReferenceRadius_ * jetReferenceRadius_) * nPV / vertexRecoEffcy_ *
(hfTowerEtaWidth_ * hfTowerPhiWidth_);

for (auto const& jet : *theJets) {
double pt_jet = jet.pt();
double eta_jet = jet.eta();

//If central or low pt jets, fill with dummy variables
if (pt_jet <= jetPtThreshold_ || std::abs(eta_jet) <= jetEtaThreshold_) {
v_sigmaEtaEta.push_back(-1.);
v_sigmaPhiPhi.push_back(-1.);
v_size_CentralEtaStrip.push_back(0);
v_size_AdjacentEtaStrips.push_back(0);
} else {
//First loop over PF candidates to compute some global variables needed for shower shape calculations
double sumptPFcands = 0.;

for (unsigned i = 0; i < jet.numberOfSourceCandidatePtrs(); ++i) {
const reco::Candidate* icand = jet.sourceCandidatePtr(i).get();
//Only look at pdgId =1,2 (HF PF cands)
if (std::abs(icand->pdgId()) > 2)
continue;
double pt_PUsub = icand->pt() - puoffset;
if (pt_PUsub < widthPtThreshold_)
continue;
sumptPFcands += pt_PUsub;
}

//Second loop to compute the various shower shape variables
int size_CentralEtaStrip(0.), size_AdjacentEtaStrips(0.);
double sigmaEtaEtaSq(0.), sigmaPhiPhiSq(0.);
double sumweightsPFcands = 0;
for (unsigned i = 0; i < jet.numberOfSourceCandidatePtrs(); ++i) {
const reco::Candidate* icand = jet.sourceCandidatePtr(i).get();
//Only look at pdgId =1,2 (HF PF cands)
if (std::abs(icand->pdgId()) > 2)
continue;

double deta = std::abs(icand->eta() - jet.eta());
double dphi = std::abs(deltaPhi(icand->phi(), jet.phi()));
double pt_PUsub = icand->pt() - puoffset;

//This is simply the size of the central eta strip and the adjacent strips
if (pt_PUsub >= stripPtThreshold_) {
if (dphi <= hfTowerPhiWidth_ * 0.5)
size_CentralEtaStrip++;
else if (dphi <= hfTowerPhiWidth_ * 1.5)
size_AdjacentEtaStrips++;
}

//Now computing sigmaEtaEta/PhiPhi
if (pt_PUsub >= widthPtThreshold_ && sumptPFcands > 0) {
double weight = pt_PUsub / sumptPFcands;
sigmaEtaEtaSq += deta * deta * weight;
sigmaPhiPhiSq += dphi * dphi * weight;
sumweightsPFcands += weight;
}
}

v_size_CentralEtaStrip.push_back(size_CentralEtaStrip);
v_size_AdjacentEtaStrips.push_back(size_AdjacentEtaStrips);

if (sumweightsPFcands > 0 && sigmaEtaEtaSq > 0 && sigmaPhiPhiSq > 0) {
v_sigmaEtaEta.push_back(sqrt(sigmaEtaEtaSq / sumweightsPFcands));
v_sigmaPhiPhi.push_back(sqrt(sigmaPhiPhiSq / sumweightsPFcands));
} else {
v_sigmaEtaEta.push_back(-1.);
v_sigmaPhiPhi.push_back(-1.);
}

} //End loop over jets
}

putInEvent("sigmaEtaEta", theJets, v_sigmaEtaEta, iEvent);
putInEvent("sigmaPhiPhi", theJets, v_sigmaPhiPhi, iEvent);
putInEvent("centralEtaStripSize", theJets, v_size_CentralEtaStrip, iEvent);
putInEvent("adjacentEtaStripsSize", theJets, v_size_AdjacentEtaStrips, iEvent);
}

/// Function to put product into event
template <typename T>
void HFJetShowerShape::putInEvent(const std::string& name,
const edm::Handle<edm::View<reco::Jet>>& jets,
std::vector<T> product,
edm::Event& iEvent) const {
auto out = std::make_unique<edm::ValueMap<T>>();
typename edm::ValueMap<T>::Filler filler(*out);
filler.insert(jets, product.begin(), product.end());
filler.fill();
iEvent.put(std::move(out), name);
}

// ------------ method fills 'descriptions' with the allowed parameters for the module ------------
void HFJetShowerShape::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
desc.add<edm::InputTag>("vertices", edm::InputTag("offlineSlimmedPrimaryVertices"));
desc.add<double>("jetPtThreshold", 25.);
desc.add<double>("jetEtaThreshold", 2.9);
desc.add<double>("hfTowerEtaWidth", 0.175);
desc.add<double>("hfTowerPhiWidth", 0.175);
desc.add<double>("vertexRecoEffcy", 0.7);
desc.add<double>("offsetPerPU", 0.4);
desc.add<double>("jetReferenceRadius", 0.4);
desc.add<double>("stripPtThreshold", 10.);
desc.add<double>("widthPtThreshold", 3.);
descriptions.add("hfJetShowerShape", desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(HFJetShowerShape);
10 changes: 10 additions & 0 deletions RecoJets/JetProducers/test/testJetTools_cfg.py
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,16 @@
patJetsAK4.userData.userFloats.src += ['QGTagger:qgLikelihood']
process.out.outputCommands += ['keep *_QGTagger_*_*']

#HF shower shape variables
process.load('RecoJets.JetProducers.hfJetShowerShape_cfi')
process.hfJetShowerShape.jets = cms.InputTag("ak4PFJetsCHS")
process.hfJetShowerShape.vertices = cms.InputTag("offlinePrimaryVertices")

patJetsAK4.userData.userFloats.src += ['hfJetShowerShape:sigmaEtaEta','hfJetShowerShape:sigmaPhiPhi']
patJetsAK4.userData.userInts.src += ['hfJetShowerShape:centralEtaStripSize','hfJetShowerShape:adjacentEtaStripsSize']
process.out.outputCommands += ['keep *_hfJetShowerShape_*_*']


#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#Njettiness

Expand Down
11 changes: 11 additions & 0 deletions RecoJets/JetProducers/test/testJetTools_onMiniAOD_cfg.py
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,17 @@
patJetsAK4.userData.userFloats.src += ['QGTagger:qgLikelihood']
process.out.outputCommands += ['keep *_QGTagger_*_*']


#HF shower shape variables
process.load('RecoJets.JetProducers.hfJetShowerShape_cfi')
patAlgosToolsTask.add(process.hfJetShowerShape)
process.hfJetShowerShape.jets=cms.InputTag("slimmedJets")
process.hfJetShowerShape.vertices=cms.InputTag("offlineSlimmedPrimaryVertices")
patJetsAK4.userData.userFloats.src += ['hfJetShowerShape:sigmaEtaEta','hfJetShowerShape:sigmaPhiPhi']
patJetsAK4.userData.userInts.src += ['hfJetShowerShape:centralEtaStripSize','hfJetShowerShape:adjacentEtaStripsSize']
process.out.outputCommands += ['keep *_hfJetShowerShape_*_*']


#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#Njettiness

Expand Down

0 comments on commit d8780f6

Please sign in to comment.