Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Puppi pileupjetID #24287

Merged
merged 8 commits into from Aug 21, 2018
Merged

Puppi pileupjetID #24287

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
d019c48
update pileupjetId to include puppi weights
knash Aug 14, 2018
66a2709
few small things
knash Aug 14, 2018
3bf728f
add puppi weights for pull variables
knash Aug 14, 2018
68d0925
final changes
knash Aug 14, 2018
39d96a0
actual final changes
knash Aug 14, 2018
071f9f8
few whitespace issues
knash Aug 14, 2018
1afcc6b
few whitespace issues
knash Aug 14, 2018
e5fdda6
calculate puppi mult by hand
knash Aug 15, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion RecoJets/JetProducers/interface/PileupJetIdAlgo.h
View file
Open in desktop
Expand Up @@ -34,7 +34,7 @@ class PileupJetIdAlgo {
~PileupJetIdAlgo();

PileupJetIdentifier computeIdVariables(const reco::Jet * jet,
float jec, const reco::Vertex *, const reco::VertexCollection &, double rho);
float jec, const reco::Vertex *, const reco::VertexCollection &, double rho, bool usePuppi);

void set(const PileupJetIdentifier &);
float getMVAval(const std::vector<std::string> &, const std::unique_ptr<const GBRForest> &);
Expand Down
5 changes: 3 additions & 2 deletions RecoJets/JetProducers/plugins/PileupJetIdProducer.cc
View file
Open in desktop
Expand Up @@ -24,7 +24,8 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig)
inputIsCorrected_(iConfig.getParameter<bool>("inputIsCorrected")),
applyJec_(iConfig.getParameter<bool>("applyJec")),
jec_(iConfig.getParameter<std::string>("jec")),
residualsFromTxt_(iConfig.getParameter<bool>("residualsFromTxt")) {
residualsFromTxt_(iConfig.getParameter<bool>("residualsFromTxt")),
usePuppi_(iConfig.getParameter<bool>("usePuppi")) {

if (residualsFromTxt_) {
residualsTxt_ = iConfig.getParameter<edm::FileInPath>("residualsTxt");
Expand Down Expand Up @@ -167,7 +168,7 @@ PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
PileupJetIdentifier puIdentifier;
if( gc->produceJetIds() ) {
// Compute the input variables
puIdentifier = ialgo->computeIdVariables(theJet, jec, &(*vtx), vertexes, rho);
puIdentifier = ialgo->computeIdVariables(theJet, jec, &(*vtx), vertexes, rho,gc->usePuppi());
ids.push_back( puIdentifier );
} else {
// Or read it from the value map
Expand Down
2 changes: 2 additions & 0 deletions RecoJets/JetProducers/plugins/PileupJetIdProducer.h
View file
Open in desktop
Expand Up @@ -63,6 +63,7 @@ class GBRForestsAndConstants {
std::string const& jec() const { return jec_; }
bool residualsFromTxt() const { return residualsFromTxt_; }
edm::FileInPath const& residualsTxt() const { return residualsTxt_; }
bool usePuppi() const { return usePuppi_; }

private:

Expand All @@ -75,6 +76,7 @@ class GBRForestsAndConstants {
std::string jec_;
bool residualsFromTxt_;
edm::FileInPath residualsTxt_;
bool usePuppi_;
};

class PileupJetIdProducer : public edm::stream::EDProducer<edm::GlobalCache<GBRForestsAndConstants>> {
Expand Down
1 change: 1 addition & 0 deletions RecoJets/JetProducers/python/PileupJetID_cfi.py
View file
Open in desktop
Expand Up @@ -27,6 +27,7 @@
applyJec = cms.bool(True),
inputIsCorrected = cms.bool(False),
residualsFromTxt = cms.bool(False),
usePuppi = cms.bool(False),
# residualsTxt = cms.FileInPath("RecoJets/JetProducers/data/download.url") # must be an existing file
)

Expand Down
35 changes: 28 additions & 7 deletions RecoJets/JetProducers/src/PileupJetIdAlgo.cc
View file
Open in desktop
Expand Up @@ -265,7 +265,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeMva()

// ------------------------------------------------------------------------------------------
PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, float jec, const reco::Vertex * vtx,
const reco::VertexCollection & allvtx, double rho)
const reco::VertexCollection & allvtx, double rho, bool usePuppi)
{

static std::atomic<int> printWarning{10};
Expand Down Expand Up @@ -299,13 +299,14 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, f
TMatrixDSym covMatrix(2); covMatrix = 0.;
float jetPt = jet->pt() / jec; // use uncorrected pt for shape variables
float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
float multNeut = 0.0;
setPtEtaPhi(*jet,internalId_.jetPt_,internalId_.jetEta_,internalId_.jetPhi_); // use corrected pt for jet kinematics
internalId_.jetM_ = jet->mass();
internalId_.nvtx_ = allvtx.size();
internalId_.rho_ = rho;

float dRmin(1000);

for ( unsigned i = 0; i < jet->numberOfSourceCandidatePtrs(); ++i ) {
reco::CandidatePtr pfJetConstituent = jet->sourceCandidatePtr(i);

Expand All @@ -316,8 +317,9 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, f
if (lPack == nullptr){
isPacked = false;
}

float candPt = icand->pt();
float candPuppiWeight = 1.0;
if (usePuppi && isPacked) candPuppiWeight = lPack->puppiWeight();
float candPt = (icand->pt())*candPuppiWeight;
float candPtFrac = candPt/jetPt;
float candDr = reco::deltaR(*icand,*jet);
float candDeta = icand->eta() - jet->eta();
Expand Down Expand Up @@ -358,6 +360,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, f
if( icone < ncones ) { *coneNeutFracs[icone] += candPt; }
internalId_.ptDNe_ += candPt*candPt;
sumPtNe += candPt;
multNeut += candPuppiWeight;
}
// EM candidated
if( icand->pdgId() == 22 ) {
Expand All @@ -367,7 +370,10 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, f
if( icone < ncones ) { *coneEmFracs[icone] += candPt; }
internalId_.ptDNe_ += candPt*candPt;
sumPtNe += candPt;
multNeut += candPuppiWeight;
}
if((abs(icand->pdgId()) == 1) || (abs(icand->pdgId()) == 2)) multNeut += candPuppiWeight;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the purpose of this line? PdgId 1,2 seem not to be used in PF? https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideParticleFlow#Particle_types

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are HF had/gamma. The multNeut should reproduce patjet->neutralMultiplicity() exactly in the case that usePuppi is False.

So we have events like this (setting usePuppi to False)

Adding 130
Adding 130
Adding 22
Adding 22
Adding 22
Adding 22
Adding 22
Adding 22
multNeut = 8 ; patjet->neutralMultiplicity() = 8

and events like this

Adding 1
Adding 1
Adding 1
Adding 1
Adding 1
multNeut = 5 ; patjet->neutralMultiplicity() = 5

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok makes sense!


// Charged particles
if( icand->charge() != 0 ) {
if (lLeadCh == nullptr || candPt > lLeadCh->pt()) {
Expand Down Expand Up @@ -483,6 +489,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, f
}
}
}

// trailing candidate
if( lTrail == nullptr || candPt < lTrail->pt() ) {
lTrail = icand;
Expand All @@ -505,6 +512,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, f
internalId_.neuEMfrac_ = patjet->neutralEmEnergy() /jet->energy();
internalId_.chgHadrfrac_ = patjet->chargedHadronEnergy()/jet->energy();
internalId_.neuHadrfrac_ = patjet->neutralHadronEnergy()/jet->energy();
if (usePuppi) internalId_.nNeutrals_ = multNeut;
} else {
internalId_.nCharged_ = pfjet->chargedMultiplicity();
internalId_.nNeutrals_ = pfjet->neutralMultiplicity();
Expand All @@ -525,7 +533,13 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, f
continue;
}

float weight = part->pt();
float partPuppiWeight=1.0;
if (usePuppi){
const pat::PackedCandidate* partpack = dynamic_cast<const pat::PackedCandidate *>( part.get() );
if (partpack!=nullptr) partPuppiWeight = partpack->puppiWeight();
}

float weight = (part->pt())*partPuppiWeight;
float weight2 = weight * weight;
sumW2 += weight2;
float deta = part->eta() - jet->eta();
Expand All @@ -544,7 +558,14 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, f
if (!(part.isAvailable() && part.isNonnull()) ){
continue;
}
float weight =part->pt()*part->pt();

float partPuppiWeight=1.0;
if (usePuppi){
const pat::PackedCandidate* partpack = dynamic_cast<const pat::PackedCandidate *>( part.get() );
if (partpack!=nullptr) partPuppiWeight = partpack->puppiWeight();
}

float weight = partPuppiWeight*(part->pt())*partPuppiWeight*(part->pt());
float deta = part->eta() - jet->eta();
float dphi = reco::deltaPhi(*part, *jet);
float ddeta, ddphi, ddR;
Expand Down Expand Up @@ -597,7 +618,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet * jet, f
internalId_.dRMeanEm_ /= jetPt;
internalId_.dRMeanCh_ /= jetPt;
internalId_.dR2Mean_ /= sumPt2;

for(size_t ic=0; ic<ncones; ++ic){
*coneFracs[ic] /= jetPt;
*coneEmFracs[ic] /= jetPt;
Expand Down
17 changes: 9 additions & 8 deletions RecoJets/JetProducers/test/testJetTools_onMiniAOD_cfg.py
View file
Open in desktop
Expand Up @@ -75,15 +75,15 @@
#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#ECF

process.load('RecoJets.JetProducers.ECF_cfi')
patAlgosToolsTask.add(process.ECF)
process.ECFAK8 = process.ECF.clone()
patAlgosToolsTask.add(process.ECFAK8)
process.ECFAK8.cone = cms.double(0.8)
process.ECFAK8.src = cms.InputTag("slimmedJetsAK8")
process.load('RecoJets.JetProducers.ECF_cff')
patAlgosToolsTask.add(process.ecf)
process.ecfAK8 = process.ecf.clone()
patAlgosToolsTask.add(process.ecfAK8)
#process.ecfAK8.cone = cms.double(0.8)
process.ecfAK8.src = cms.InputTag("slimmedJetsAK8")

patJetsAK8.userData.userFloats.src += ['ECFAK8:ecf1','ECFAK8:ecf2','ECFAK8:ecf3']
process.out.outputCommands += ['keep *_ECFAK8_*_*']
patJetsAK8.userData.userFloats.src += ['ecfAK8:ecf1','ecfAK8:ecf2','ecfAK8:ecf3']
process.out.outputCommands += ['keep *_ecfAK8_*_*']

#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#QJetsAdder
Expand Down Expand Up @@ -172,3 +172,4 @@
process.out.fileName = 'testJetTools.root'
# ##
# process.options.wantSummary = False ## (to suppress the long output at the end of the job)