From 703a7cdd4118dc05c6571e4216a35155583998aa Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Mon, 13 Feb 2023 00:29:49 +0900 Subject: [PATCH 1/7] Fix PUID Variable bug and implement PUPPI Weight using a Value Map --- PhysicsTools/NanoAOD/python/custom_jme_cff.py | 2 +- .../JetProducers/interface/PileupJetIdAlgo.h | 2 +- .../plugins/PileupJetIdProducer.cc | 28 ++- .../plugins/PileupJetIdProducer.h | 6 +- .../JetProducers/python/PileupJetID_cfi.py | 2 +- RecoJets/JetProducers/src/PileupJetIdAlgo.cc | 209 +++++++++++------- 6 files changed, 164 insertions(+), 85 deletions(-) diff --git a/PhysicsTools/NanoAOD/python/custom_jme_cff.py b/PhysicsTools/NanoAOD/python/custom_jme_cff.py index 7c4494c6c9c02..ab920cfaab0f0 100644 --- a/PhysicsTools/NanoAOD/python/custom_jme_cff.py +++ b/PhysicsTools/NanoAOD/python/custom_jme_cff.py @@ -274,7 +274,7 @@ def AddPileUpJetIDVars(proc, jetName="", jetSrc="", jetTableName="", jetTaskName vertexes = "offlineSlimmedPrimaryVertices", inputIsCorrected = True, applyJec = False, - usePuppi = True if "PUPPI" in jetName.upper() else False + srcConstituentWeights = "packedPFCandidatespuppi" if "PUPPI" in jetName.upper() else "" ) ) getattr(proc,jetTaskName).add(getattr(proc, puJetIdVarsCalculator)) diff --git a/RecoJets/JetProducers/interface/PileupJetIdAlgo.h b/RecoJets/JetProducers/interface/PileupJetIdAlgo.h index 3969257b30c4b..fa46d1b6b9b86 100644 --- a/RecoJets/JetProducers/interface/PileupJetIdAlgo.h +++ b/RecoJets/JetProducers/interface/PileupJetIdAlgo.h @@ -29,7 +29,7 @@ class PileupJetIdAlgo { ~PileupJetIdAlgo(); PileupJetIdentifier computeIdVariables( - const reco::Jet* jet, float jec, const reco::Vertex*, const reco::VertexCollection&, double rho, bool usePuppi); + const reco::Jet* jet, float jec, const reco::Vertex*, const reco::VertexCollection&, double rho, edm::ValueMap& constituentWeights , bool applyConstituentWeight); void set(const PileupJetIdentifier&); float getMVAval(const std::vector&, const std::unique_ptr&); diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc index 2edfeb24729d7..e295394285549 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -27,7 +27,7 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig) applyJec_(iConfig.getParameter("applyJec")), jec_(iConfig.getParameter("jec")), residualsFromTxt_(iConfig.getParameter("residualsFromTxt")), - usePuppi_(iConfig.getParameter("usePuppi")) { + applyConstituentWeight_(false){ if (residualsFromTxt_) { residualsTxt_ = iConfig.getParameter("residualsTxt"); } @@ -40,6 +40,12 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig) if (!runMvas_) { assert(algos.size() == 1); } + + edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); + if (!srcConstituentWeights.label().empty()){ + applyConstituentWeight_ = true; + } + } // ------------------------------------------------------------------------------------------ @@ -62,6 +68,17 @@ PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRFo consumes>(iConfig.getParameter("jetids")); input_rho_token_ = consumes(iConfig.getParameter("rho")); parameters_token_ = esConsumes(edm::ESInputTag("", globalCache->jec())); + + +// edm::InputTag src_ = iConfig.getParameter("src"); +// edm::InputTag srcWeights = iConfig.getParameter("srcWeights"); +// if (iConfig.getParameter("src").label() == iConfig.getParameter("srcConstituentWeights").label()) +// edm::LogWarning("PileupJetIdAlgo") +// << "Particle and weights collection have the same label. You may be applying the same weights twice. \n"; + edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); + if (!srcConstituentWeights.label().empty()){ + input_constituent_weights_token_ = consumes>(srcConstituentWeights); + } } // ------------------------------------------------------------------------------------------ @@ -80,6 +97,12 @@ void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe iEvent.getByToken(input_jet_token_, jetHandle); const View& jets = *jetHandle; + // Constituent weight (e.g PUPPI) Value Map + edm::ValueMap constituentWeights; + if (!input_constituent_weights_token_.isUninitialized()) { + constituentWeights = iEvent.get(input_constituent_weights_token_); + } + // input variables Handle> vmap; if (!gc->produceJetIds()) { @@ -167,7 +190,8 @@ void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe PileupJetIdentifier puIdentifier; if (gc->produceJetIds()) { // Compute the input variables - puIdentifier = ialgo->computeIdVariables(theJet, jec, &(*vtx), *vertexes, rho, gc->usePuppi()); + ////////////////////////////// added PUPPI weight Value Map + puIdentifier = ialgo->computeIdVariables(theJet, jec, &(*vtx), *vertexes, rho, constituentWeights, gc->applyConstituentWeight()); ids.push_back(puIdentifier); } else { // Or read it from the value map diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.h b/RecoJets/JetProducers/plugins/PileupJetIdProducer.h index 479179f175114..5bf91438b9c5c 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.h +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.h @@ -61,7 +61,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_; } + bool applyConstituentWeight() const { return applyConstituentWeight_; } private: std::vector vAlgoGBRForestsAndConstants_; @@ -73,7 +73,7 @@ class GBRForestsAndConstants { std::string jec_; bool residualsFromTxt_; edm::FileInPath residualsTxt_; - bool usePuppi_; + bool applyConstituentWeight_; }; class PileupJetIdProducer : public edm::stream::EDProducer> { @@ -99,6 +99,8 @@ class PileupJetIdProducer : public edm::stream::EDProducer jecCor_; std::vector jetCorPars_; + edm::ValueMap constituentWeights_; + edm::EDGetTokenT> input_constituent_weights_token_; edm::EDGetTokenT> input_jet_token_; edm::EDGetTokenT input_vertex_token_; edm::EDGetTokenT> input_vm_pujetid_token_; diff --git a/RecoJets/JetProducers/python/PileupJetID_cfi.py b/RecoJets/JetProducers/python/PileupJetID_cfi.py index 1fdafecaca58c..f326da1846e60 100644 --- a/RecoJets/JetProducers/python/PileupJetID_cfi.py +++ b/RecoJets/JetProducers/python/PileupJetID_cfi.py @@ -32,7 +32,7 @@ applyJec = cms.bool(True), inputIsCorrected = cms.bool(False), residualsFromTxt = cms.bool(False), - usePuppi = cms.bool(False), + srcConstituentWeights = cms.InputTag(""), # residualsTxt = cms.FileInPath("RecoJets/JetProducers/data/download.url") # must be an existing file ) diff --git a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc index b31a33f6d8b9f..b39bd1fbb0849 100644 --- a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc +++ b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc @@ -263,7 +263,8 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, const reco::Vertex* vtx, const reco::VertexCollection& allvtx, double rho, - bool usePuppi) { + edm::ValueMap& constituentWeights, + bool applyConstituentWeight) { // initialize all variables to 0 resetVariables(); @@ -280,6 +281,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, const reco::Candidate *lLead = nullptr, *lSecond = nullptr, *lLeadNeut = nullptr, *lLeadEm = nullptr, *lLeadCh = nullptr, *lTrail = nullptr; + std::vector frac, fracCh, fracEm, fracNeut; float cones[] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7}; size_t ncones = sizeof(cones) / sizeof(float); @@ -316,6 +318,9 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, 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; + float sumW2(0.0); + float sum_deta(0.0), sum_dphi(0.0); + float ave_deta(0.0), ave_dphi(0.0); setPtEtaPhi( *jet, internalId_.jetPt_, internalId_.jetEta_, internalId_.jetPhi_); // use corrected pt for jet kinematics internalId_.jetM_ = jet->mass(); @@ -323,10 +328,16 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.rho_ = rho; float dRmin(1000); + float LeadcandWeight = 1.0; + float SecondcandWeight = 1.0; + float LeadNeutcandWeight = 1.0; + float LeadEmcandWeight = 1.0; + float LeadChcandWeight = 1.0; + float TrailcandWeight = 1.0; + for (unsigned i = 0; i < jet->numberOfSourceCandidatePtrs(); ++i) { reco::CandidatePtr pfJetConstituent = jet->sourceCandidatePtr(i); - const reco::Candidate* icand = pfJetConstituent.get(); const pat::PackedCandidate* lPack = dynamic_cast(icand); const reco::PFCandidate* lPF = dynamic_cast(icand); @@ -334,10 +345,12 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, if (lPack == nullptr) { isPacked = false; } - float candPuppiWeight = 1.0; - if (usePuppi && isPacked) - candPuppiWeight = lPack->puppiWeight(); - float candPt = (icand->pt()) * candPuppiWeight; + + float candWeight = 1.0; + if (applyConstituentWeight){ // PUPPI Jet weight should be pulled up from valuemap, not packed candidate + candWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + float candPt = (icand->pt()) * candWeight; float candPtFrac = candPt / jetPt; float candDr = reco::deltaR(*icand, *jet); float candDeta = icand->eta() - jet->eta(); @@ -348,12 +361,19 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, if (candDr < dRmin) dRmin = candDr; - // // all particles - if (lLead == nullptr || candPt > lLead->pt()) { + // // all particles; PUPPI weights multiplied to leading and subleading constituent if it is for PUPPI + if (lLead == nullptr || candPt > (lLead->pt())*LeadcandWeight) { lSecond = lLead; + SecondcandWeight = LeadcandWeight; lLead = icand; - } else if ((lSecond == nullptr || candPt > lSecond->pt()) && (candPt < lLead->pt())) { + if (applyConstituentWeight) { + LeadcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + } else if ((lSecond == nullptr || candPt > (lSecond->pt())*SecondcandWeight) && (candPt < (lLead->pt())*LeadcandWeight)) { lSecond = icand; + if (applyConstituentWeight) { + SecondcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } } // // average shapes @@ -373,11 +393,15 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, *coneFracs[icone] += candPt; } - // neutrals + // neutrals Neutral hadrons if (abs(icand->pdgId()) == 130) { - if (lLeadNeut == nullptr || candPt > lLeadNeut->pt()) { + if (lLeadNeut == nullptr || candPt > (lLeadNeut->pt())*LeadNeutcandWeight) { lLeadNeut = icand; + if (applyConstituentWeight) { + LeadNeutcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } } + internalId_.dRMeanNeut_ += candPtDr; fracNeut.push_back(candPtFrac); if (icone < ncones) { @@ -385,12 +409,16 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } internalId_.ptDNe_ += candPt * candPt; sumPtNe += candPt; - multNeut += candPuppiWeight; + multNeut += candWeight; } - // EM candidated + + // EM candidated photon if (icand->pdgId() == 22) { - if (lLeadEm == nullptr || candPt > lLeadEm->pt()) { + if (lLeadEm == nullptr || candPt > (lLeadEm->pt())*LeadEmcandWeight) { lLeadEm = icand; + if (applyConstituentWeight) { + LeadEmcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } } internalId_.dRMeanEm_ += candPtDr; fracEm.push_back(candPtFrac); @@ -399,16 +427,19 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } internalId_.ptDNe_ += candPt * candPt; sumPtNe += candPt; - multNeut += candPuppiWeight; + multNeut += candWeight; } + // hadrons and EM in HF if ((abs(icand->pdgId()) == 1) || (abs(icand->pdgId()) == 2)) - multNeut += candPuppiWeight; + multNeut += candWeight; // Charged particles if (icand->charge() != 0) { - if (lLeadCh == nullptr || candPt > lLeadCh->pt()) { + if (lLeadCh == nullptr || candPt > (lLeadCh->pt())*LeadChcandWeight) { lLeadCh = icand; - + if (applyConstituentWeight) { + LeadChcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } const reco::Track* pfTrk = icand->bestTrack(); if (lPF && std::abs(icand->pdgId()) == 13 && pfTrk == nullptr) { reco::MuonRef lmuRef = lPF->muonRef(); @@ -491,7 +522,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, double dZ0 = 9999.; double dZ_tmp = 9999.; for (unsigned vtx_i = 0; vtx_i < allvtx.size(); vtx_i++) { - const auto& iv = allvtx[vtx_i]; + auto iv = allvtx[vtx_i]; if (iv.isFake()) continue; @@ -525,27 +556,79 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } // trailing candidate - if (lTrail == nullptr || candPt < lTrail->pt()) { + if (lTrail == nullptr || candPt < (lTrail->pt())*TrailcandWeight) { lTrail = icand; + if (applyConstituentWeight) { + TrailcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + } + + // average for pull variable + + float weight2 = candPt * candPt; + sumW2 += weight2; + float deta = icand->eta() - jet->eta(); + float dphi = reco::deltaPhi(*icand, *jet); + sum_deta += deta * weight2; + sum_dphi += dphi * weight2; + if (sumW2 > 0) { + ave_deta = sum_deta / sumW2; + ave_dphi = sum_dphi / sumW2; } } // // Finalize all variables - assert(!(lLead == nullptr)); + // Most of Below values are not used for puID variable generation at the moment, except lLeadCh Pt for JetRchg, so I assign that zero if there is no charged constituent. - if (lSecond == nullptr) { - lSecond = lTrail; + + assert(!(lLead == nullptr)); + internalId_.leadPt_ = lLead->pt()*LeadcandWeight; + internalId_.leadEta_ = lLead->eta(); + internalId_.leadPhi_ = lLead->phi(); + + if (lSecond != nullptr){ + internalId_.secondPt_ = lSecond->pt()*SecondcandWeight; + internalId_.secondEta_ = lSecond->eta(); + internalId_.secondPhi_ = lSecond->phi(); + } else { + internalId_.secondPt_ = 0.0; + internalId_.secondEta_ = large_val; + internalId_.secondPhi_ = large_val; } - if (lLeadNeut == nullptr) { - lLeadNeut = lTrail; + + if (lLeadNeut != nullptr){ + internalId_.leadNeutPt_ = lLeadNeut->pt()*LeadNeutcandWeight; + internalId_.leadNeutEta_ = lLeadNeut->eta(); + internalId_.leadNeutPhi_ = lLeadNeut->phi(); + } else { + internalId_.leadNeutPt_ = 0.0; + internalId_.leadNeutEta_ = large_val; + internalId_.leadNeutPhi_ = large_val; } - if (lLeadEm == nullptr) { - lLeadEm = lTrail; + + if (lLeadEm != nullptr){ + internalId_.leadEmPt_ = lLeadEm->pt()*LeadEmcandWeight; + internalId_.leadEmEta_ = lLeadEm->eta(); + internalId_.leadEmPhi_ = lLeadEm->phi(); + } else { + internalId_.leadEmPt_ = 0.0; + internalId_.leadEmEta_ = large_val; + internalId_.leadEmPhi_ = large_val; } - if (lLeadCh == nullptr) { - lLeadCh = lTrail; + + if (lLeadCh != nullptr){ + internalId_.leadChPt_ = lLeadCh->pt()*LeadChcandWeight; + internalId_.leadChEta_ = lLeadCh->eta(); + internalId_.leadChPhi_ = lLeadCh->phi(); + } else { + internalId_.leadChPt_ = 0.0; + internalId_.leadChEta_ = large_val; + internalId_.leadChPhi_ = large_val; } + + + if (patjet != nullptr) { // to enable running on MiniAOD slimmedJets internalId_.nCharged_ = patjet->chargedMultiplicity(); internalId_.nNeutrals_ = patjet->neutralMultiplicity(); @@ -553,7 +636,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.neuEMfrac_ = patjet->neutralEmEnergy() / jet->energy(); internalId_.chgHadrfrac_ = patjet->chargedHadronEnergy() / jet->energy(); internalId_.neuHadrfrac_ = patjet->neutralHadronEnergy() / jet->energy(); - if (usePuppi) + if (applyConstituentWeight) internalId_.nNeutrals_ = multNeut; } else { internalId_.nCharged_ = pfjet->chargedMultiplicity(); @@ -562,54 +645,24 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.neuEMfrac_ = pfjet->neutralEmEnergy() / jet->energy(); internalId_.chgHadrfrac_ = pfjet->chargedHadronEnergy() / jet->energy(); internalId_.neuHadrfrac_ = pfjet->neutralHadronEnergy() / jet->energy(); + if (applyConstituentWeight) + internalId_.nNeutrals_ = multNeut; } internalId_.nParticles_ = jet->nConstituents(); ///////////////////////pull variable/////////////////////////////////// - float sumW2(0.0); - float sum_deta(0.0), sum_dphi(0.0); - float ave_deta(0.0), ave_dphi(0.0); - for (size_t j = 0; j < jet->numberOfDaughters(); j++) { - const auto& part = jet->daughterPtr(j); - if (!(part.isAvailable() && part.isNonnull())) { - continue; - } - - float partPuppiWeight = 1.0; - if (usePuppi) { - const pat::PackedCandidate* partpack = dynamic_cast(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(); - float dphi = reco::deltaPhi(*part, *jet); - sum_deta += deta * weight2; - sum_dphi += dphi * weight2; - if (sumW2 > 0) { - ave_deta = sum_deta / sumW2; - ave_dphi = sum_dphi / sumW2; - } - } - float ddetaR_sum(0.0), ddphiR_sum(0.0), pull_tmp(0.0); - for (size_t i = 0; i < jet->numberOfDaughters(); i++) { - const auto& part = jet->daughterPtr(i); - if (!(part.isAvailable() && part.isNonnull())) { - continue; - } + for (unsigned k = 0; k < jet->numberOfSourceCandidatePtrs(); k++) { + reco::CandidatePtr temp_pfJetConstituent = jet->sourceCandidatePtr(k); +// reco::CandidatePtr temp_weightpfJetConstituent = jet->sourceCandidatePtr(k); + const reco::Candidate* part = temp_pfJetConstituent.get(); + + float candWeight = 1.0; - float partPuppiWeight = 1.0; - if (usePuppi) { - const pat::PackedCandidate* partpack = dynamic_cast(part.get()); - if (partpack != nullptr) - partPuppiWeight = partpack->puppiWeight(); - } + if (applyConstituentWeight) + candWeight = constituentWeights[jet->sourceCandidatePtr(k)]; - float weight = partPuppiWeight * (part->pt()) * partPuppiWeight * (part->pt()); + float weight = candWeight * (part->pt()) * candWeight * (part->pt()); float deta = part->eta() - jet->eta(); float dphi = reco::deltaPhi(*part, *jet); float ddeta, ddphi, ddR; @@ -627,11 +680,6 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.pull_ = pull_tmp; /////////////////////////////////////////////////////////////////////// - setPtEtaPhi(*lLead, internalId_.leadPt_, internalId_.leadEta_, internalId_.leadPhi_); - setPtEtaPhi(*lSecond, internalId_.secondPt_, internalId_.secondEta_, internalId_.secondPhi_); - setPtEtaPhi(*lLeadNeut, internalId_.leadNeutPt_, internalId_.leadNeutEta_, internalId_.leadNeutPhi_); - setPtEtaPhi(*lLeadEm, internalId_.leadEmPt_, internalId_.leadEmEta_, internalId_.leadEmPhi_); - setPtEtaPhi(*lLeadCh, internalId_.leadChPt_, internalId_.leadChEta_, internalId_.leadChPhi_); std::sort(frac.begin(), frac.end(), std::greater()); std::sort(fracCh.begin(), fracCh.end(), std::greater()); @@ -698,8 +746,13 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.sumChPt_ = sumPtCh; internalId_.sumNePt_ = sumPtNe; - internalId_.jetR_ = lLead->pt() / sumPt; - internalId_.jetRchg_ = lLeadCh->pt() / sumPt; + internalId_.jetR_ = (lLead->pt())*LeadcandWeight / sumPt; + if (lLeadCh != nullptr){ + internalId_.jetRchg_ = (lLeadCh->pt())*LeadChcandWeight / sumPt; + } else { + internalId_.jetRchg_ = 0; + } + internalId_.dRMatch_ = dRmin; if (sumTkPt != 0.) { From 7174e335d5a1c7081c9e8ad7628518ec282bd89b Mon Sep 17 00:00:00 2001 From: Nurfikri Norjoharuddeen Date: Mon, 13 Feb 2023 07:22:01 +0100 Subject: [PATCH 2/7] Make photon protection for existing puppi weights an option. Default is switched off. --- CommonTools/PileupAlgos/plugins/PuppiProducer.cc | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc index ff386726f54c4..515fd2b1e125b 100644 --- a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc +++ b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc @@ -80,6 +80,7 @@ class PuppiProducer : public edm::stream::EDProducer<> { uint fNumOfPUVtxsForCharged; double fDZCutForChargedFromPUVtxs; bool fUseExistingWeights; + bool fApplyPhotonProtectionForExistingWeights; bool fClonePackedCands; int fVtxNdofCut; double fVtxZCut; @@ -104,6 +105,7 @@ PuppiProducer::PuppiProducer(const edm::ParameterSet& iConfig) { fNumOfPUVtxsForCharged = iConfig.getParameter("NumOfPUVtxsForCharged"); fDZCutForChargedFromPUVtxs = iConfig.getParameter("DeltaZCutForChargedFromPUVtxs"); fUseExistingWeights = iConfig.getParameter("useExistingWeights"); + fApplyPhotonProtectionForExistingWeights = iConfig.getParameter("applyPhotonProtectionForExistingWeights"); fClonePackedCands = iConfig.getParameter("clonePackedCands"); fVtxNdofCut = iConfig.getParameter("vtxNdofCut"); fVtxZCut = iConfig.getParameter("vtxZCut"); @@ -363,9 +365,9 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { curpupweight = lPack->puppiWeight(); } } - // Protect high pT photons (important for gamma to hadronic recoil balance) - if ((fPtMaxPhotons > 0) && (lPack->pdgId() == 22) && (std::abs(lPack->eta()) < fEtaMaxPhotons) && - (lPack->pt() > fPtMaxPhotons)) + // Optional: Protect high pT photons (important for gamma to hadronic recoil balance) for existing weights. + if (fApplyPhotonProtectionForExistingWeights && (fPtMaxPhotons > 0) && (lPack->pdgId() == 22) && + (std::abs(lPack->eta()) < fEtaMaxPhotons) && (lPack->pt() > fPtMaxPhotons)) curpupweight = 1; lWeights.push_back(curpupweight); lPackCtr++; @@ -511,6 +513,7 @@ void PuppiProducer::fillDescriptions(edm::ConfigurationDescriptions& description desc.add("NumOfPUVtxsForCharged", 0); desc.add("DeltaZCutForChargedFromPUVtxs", 0.2); desc.add("useExistingWeights", false); + desc.add("applyPhotonProtectionForExistingWeights", false); desc.add("clonePackedCands", false); desc.add("vtxNdofCut", 4); desc.add("vtxZCut", 24); From f3a673dcd37b27e654b0593c64e5893308437fe0 Mon Sep 17 00:00:00 2001 From: jshin96 <87310162+jshin96@users.noreply.github.com> Date: Tue, 14 Feb 2023 22:24:31 +0900 Subject: [PATCH 3/7] Update PileupJetIdProducer.cc --- RecoJets/JetProducers/plugins/PileupJetIdProducer.cc | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc index e295394285549..76bc490be9d7e 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -70,11 +70,7 @@ PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRFo parameters_token_ = esConsumes(edm::ESInputTag("", globalCache->jec())); -// edm::InputTag src_ = iConfig.getParameter("src"); -// edm::InputTag srcWeights = iConfig.getParameter("srcWeights"); -// if (iConfig.getParameter("src").label() == iConfig.getParameter("srcConstituentWeights").label()) -// edm::LogWarning("PileupJetIdAlgo") -// << "Particle and weights collection have the same label. You may be applying the same weights twice. \n"; + edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); if (!srcConstituentWeights.label().empty()){ input_constituent_weights_token_ = consumes>(srcConstituentWeights); From ec4a43ed71b8be653468fb8ef003fb31adf4b183 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Wed, 15 Feb 2023 15:26:38 +0900 Subject: [PATCH 4/7] apply code-checks patch --- .../JetProducers/interface/PileupJetIdAlgo.h | 9 ++- .../plugins/PileupJetIdProducer.cc | 17 ++--- RecoJets/JetProducers/src/PileupJetIdAlgo.cc | 65 +++++++++---------- 3 files changed, 43 insertions(+), 48 deletions(-) diff --git a/RecoJets/JetProducers/interface/PileupJetIdAlgo.h b/RecoJets/JetProducers/interface/PileupJetIdAlgo.h index fa46d1b6b9b86..aecda1c1159d7 100644 --- a/RecoJets/JetProducers/interface/PileupJetIdAlgo.h +++ b/RecoJets/JetProducers/interface/PileupJetIdAlgo.h @@ -28,8 +28,13 @@ class PileupJetIdAlgo { PileupJetIdAlgo(AlgoGBRForestsAndConstants const* cache); ~PileupJetIdAlgo(); - PileupJetIdentifier computeIdVariables( - const reco::Jet* jet, float jec, const reco::Vertex*, const reco::VertexCollection&, double rho, edm::ValueMap& constituentWeights , bool applyConstituentWeight); + PileupJetIdentifier computeIdVariables(const reco::Jet* jet, + float jec, + const reco::Vertex*, + const reco::VertexCollection&, + double rho, + edm::ValueMap& constituentWeights, + bool applyConstituentWeight); void set(const PileupJetIdentifier&); float getMVAval(const std::vector&, const std::unique_ptr&); diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc index e295394285549..fff45fccb91c9 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -27,7 +27,7 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig) applyJec_(iConfig.getParameter("applyJec")), jec_(iConfig.getParameter("jec")), residualsFromTxt_(iConfig.getParameter("residualsFromTxt")), - applyConstituentWeight_(false){ + applyConstituentWeight_(false) { if (residualsFromTxt_) { residualsTxt_ = iConfig.getParameter("residualsTxt"); } @@ -42,10 +42,9 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig) } edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); - if (!srcConstituentWeights.label().empty()){ + if (!srcConstituentWeights.label().empty()) { applyConstituentWeight_ = true; } - } // ------------------------------------------------------------------------------------------ @@ -70,11 +69,6 @@ PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRFo parameters_token_ = esConsumes(edm::ESInputTag("", globalCache->jec())); -// edm::InputTag src_ = iConfig.getParameter("src"); -// edm::InputTag srcWeights = iConfig.getParameter("srcWeights"); -// if (iConfig.getParameter("src").label() == iConfig.getParameter("srcConstituentWeights").label()) -// edm::LogWarning("PileupJetIdAlgo") -// << "Particle and weights collection have the same label. You may be applying the same weights twice. \n"; edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); if (!srcConstituentWeights.label().empty()){ input_constituent_weights_token_ = consumes>(srcConstituentWeights); @@ -97,12 +91,12 @@ void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe iEvent.getByToken(input_jet_token_, jetHandle); const View& jets = *jetHandle; - // Constituent weight (e.g PUPPI) Value Map + // Constituent weight (e.g PUPPI) Value Map edm::ValueMap constituentWeights; if (!input_constituent_weights_token_.isUninitialized()) { constituentWeights = iEvent.get(input_constituent_weights_token_); } - + // input variables Handle> vmap; if (!gc->produceJetIds()) { @@ -191,7 +185,8 @@ void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe if (gc->produceJetIds()) { // Compute the input variables ////////////////////////////// added PUPPI weight Value Map - puIdentifier = ialgo->computeIdVariables(theJet, jec, &(*vtx), *vertexes, rho, constituentWeights, gc->applyConstituentWeight()); + puIdentifier = ialgo->computeIdVariables( + theJet, jec, &(*vtx), *vertexes, rho, constituentWeights, gc->applyConstituentWeight()); ids.push_back(puIdentifier); } else { // Or read it from the value map diff --git a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc index b39bd1fbb0849..a18aea1fc11f6 100644 --- a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc +++ b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc @@ -281,7 +281,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, const reco::Candidate *lLead = nullptr, *lSecond = nullptr, *lLeadNeut = nullptr, *lLeadEm = nullptr, *lLeadCh = nullptr, *lTrail = nullptr; - + std::vector frac, fracCh, fracEm, fracNeut; float cones[] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7}; size_t ncones = sizeof(cones) / sizeof(float); @@ -335,7 +335,6 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, float LeadChcandWeight = 1.0; float TrailcandWeight = 1.0; - for (unsigned i = 0; i < jet->numberOfSourceCandidatePtrs(); ++i) { reco::CandidatePtr pfJetConstituent = jet->sourceCandidatePtr(i); const reco::Candidate* icand = pfJetConstituent.get(); @@ -347,7 +346,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } float candWeight = 1.0; - if (applyConstituentWeight){ // PUPPI Jet weight should be pulled up from valuemap, not packed candidate + if (applyConstituentWeight) { // PUPPI Jet weight should be pulled up from valuemap, not packed candidate candWeight = constituentWeights[jet->sourceCandidatePtr(i)]; } float candPt = (icand->pt()) * candWeight; @@ -362,14 +361,15 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, dRmin = candDr; // // all particles; PUPPI weights multiplied to leading and subleading constituent if it is for PUPPI - if (lLead == nullptr || candPt > (lLead->pt())*LeadcandWeight) { + if (lLead == nullptr || candPt > (lLead->pt()) * LeadcandWeight) { lSecond = lLead; SecondcandWeight = LeadcandWeight; lLead = icand; if (applyConstituentWeight) { LeadcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; } - } else if ((lSecond == nullptr || candPt > (lSecond->pt())*SecondcandWeight) && (candPt < (lLead->pt())*LeadcandWeight)) { + } else if ((lSecond == nullptr || candPt > (lSecond->pt()) * SecondcandWeight) && + (candPt < (lLead->pt()) * LeadcandWeight)) { lSecond = icand; if (applyConstituentWeight) { SecondcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; @@ -395,13 +395,13 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, // neutrals Neutral hadrons if (abs(icand->pdgId()) == 130) { - if (lLeadNeut == nullptr || candPt > (lLeadNeut->pt())*LeadNeutcandWeight) { + if (lLeadNeut == nullptr || candPt > (lLeadNeut->pt()) * LeadNeutcandWeight) { lLeadNeut = icand; if (applyConstituentWeight) { LeadNeutcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; } } - + internalId_.dRMeanNeut_ += candPtDr; fracNeut.push_back(candPtFrac); if (icone < ncones) { @@ -411,10 +411,10 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, sumPtNe += candPt; multNeut += candWeight; } - + // EM candidated photon if (icand->pdgId() == 22) { - if (lLeadEm == nullptr || candPt > (lLeadEm->pt())*LeadEmcandWeight) { + if (lLeadEm == nullptr || candPt > (lLeadEm->pt()) * LeadEmcandWeight) { lLeadEm = icand; if (applyConstituentWeight) { LeadEmcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; @@ -429,13 +429,13 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, sumPtNe += candPt; multNeut += candWeight; } - // hadrons and EM in HF + // hadrons and EM in HF if ((abs(icand->pdgId()) == 1) || (abs(icand->pdgId()) == 2)) multNeut += candWeight; // Charged particles if (icand->charge() != 0) { - if (lLeadCh == nullptr || candPt > (lLeadCh->pt())*LeadChcandWeight) { + if (lLeadCh == nullptr || candPt > (lLeadCh->pt()) * LeadChcandWeight) { lLeadCh = icand; if (applyConstituentWeight) { LeadChcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; @@ -556,7 +556,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } // trailing candidate - if (lTrail == nullptr || candPt < (lTrail->pt())*TrailcandWeight) { + if (lTrail == nullptr || candPt < (lTrail->pt()) * TrailcandWeight) { lTrail = icand; if (applyConstituentWeight) { TrailcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; @@ -580,24 +580,23 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, // // Finalize all variables // Most of Below values are not used for puID variable generation at the moment, except lLeadCh Pt for JetRchg, so I assign that zero if there is no charged constituent. - assert(!(lLead == nullptr)); - internalId_.leadPt_ = lLead->pt()*LeadcandWeight; + internalId_.leadPt_ = lLead->pt() * LeadcandWeight; internalId_.leadEta_ = lLead->eta(); internalId_.leadPhi_ = lLead->phi(); - if (lSecond != nullptr){ - internalId_.secondPt_ = lSecond->pt()*SecondcandWeight; + if (lSecond != nullptr) { + internalId_.secondPt_ = lSecond->pt() * SecondcandWeight; internalId_.secondEta_ = lSecond->eta(); internalId_.secondPhi_ = lSecond->phi(); - } else { + } else { internalId_.secondPt_ = 0.0; internalId_.secondEta_ = large_val; internalId_.secondPhi_ = large_val; } - if (lLeadNeut != nullptr){ - internalId_.leadNeutPt_ = lLeadNeut->pt()*LeadNeutcandWeight; + if (lLeadNeut != nullptr) { + internalId_.leadNeutPt_ = lLeadNeut->pt() * LeadNeutcandWeight; internalId_.leadNeutEta_ = lLeadNeut->eta(); internalId_.leadNeutPhi_ = lLeadNeut->phi(); } else { @@ -606,8 +605,8 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.leadNeutPhi_ = large_val; } - if (lLeadEm != nullptr){ - internalId_.leadEmPt_ = lLeadEm->pt()*LeadEmcandWeight; + if (lLeadEm != nullptr) { + internalId_.leadEmPt_ = lLeadEm->pt() * LeadEmcandWeight; internalId_.leadEmEta_ = lLeadEm->eta(); internalId_.leadEmPhi_ = lLeadEm->phi(); } else { @@ -616,8 +615,8 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.leadEmPhi_ = large_val; } - if (lLeadCh != nullptr){ - internalId_.leadChPt_ = lLeadCh->pt()*LeadChcandWeight; + if (lLeadCh != nullptr) { + internalId_.leadChPt_ = lLeadCh->pt() * LeadChcandWeight; internalId_.leadChEta_ = lLeadCh->eta(); internalId_.leadChPhi_ = lLeadCh->phi(); } else { @@ -626,9 +625,6 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.leadChPhi_ = large_val; } - - - if (patjet != nullptr) { // to enable running on MiniAOD slimmedJets internalId_.nCharged_ = patjet->chargedMultiplicity(); internalId_.nNeutrals_ = patjet->neutralMultiplicity(); @@ -654,10 +650,10 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, float ddetaR_sum(0.0), ddphiR_sum(0.0), pull_tmp(0.0); for (unsigned k = 0; k < jet->numberOfSourceCandidatePtrs(); k++) { reco::CandidatePtr temp_pfJetConstituent = jet->sourceCandidatePtr(k); -// reco::CandidatePtr temp_weightpfJetConstituent = jet->sourceCandidatePtr(k); + // reco::CandidatePtr temp_weightpfJetConstituent = jet->sourceCandidatePtr(k); const reco::Candidate* part = temp_pfJetConstituent.get(); - - float candWeight = 1.0; + + float candWeight = 1.0; if (applyConstituentWeight) candWeight = constituentWeights[jet->sourceCandidatePtr(k)]; @@ -680,7 +676,6 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.pull_ = pull_tmp; /////////////////////////////////////////////////////////////////////// - std::sort(frac.begin(), frac.end(), std::greater()); std::sort(fracCh.begin(), fracCh.end(), std::greater()); std::sort(fracEm.begin(), fracEm.end(), std::greater()); @@ -746,11 +741,11 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.sumChPt_ = sumPtCh; internalId_.sumNePt_ = sumPtNe; - internalId_.jetR_ = (lLead->pt())*LeadcandWeight / sumPt; - if (lLeadCh != nullptr){ - internalId_.jetRchg_ = (lLeadCh->pt())*LeadChcandWeight / sumPt; - } else { - internalId_.jetRchg_ = 0; + internalId_.jetR_ = (lLead->pt()) * LeadcandWeight / sumPt; + if (lLeadCh != nullptr) { + internalId_.jetRchg_ = (lLeadCh->pt()) * LeadChcandWeight / sumPt; + } else { + internalId_.jetRchg_ = 0; } internalId_.dRMatch_ = dRmin; From 9c6f631b27d20d1aa04555cbfa37008885bd0fe7 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Wed, 15 Feb 2023 15:47:05 +0900 Subject: [PATCH 5/7] apply code-checks patch --- RecoJets/JetProducers/plugins/PileupJetIdProducer.cc | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc index 9470ecdbc8acb..1d76b9b4d65fe 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -40,13 +40,9 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig) if (!runMvas_) { assert(algos.size() == 1); } -// edm::InputTag src_ = iConfig.getParameter("src"); -// edm::InputTag srcWeights = iConfig.getParameter("srcWeights"); -// if (iConfig.getParameter("src").label() == iConfig.getParameter("srcConstituentWeights").label()) -// edm::LogWarning("PileupJetIdAlgo") -// << "Particle and weights collection have the same label. You may be applying the same weights twice. \n"; + edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); - if (!srcConstituentWeights.label().empty()){ + if (!srcConstituentWeights.label().empty()) { applyConstituentWeight_ = true; } } @@ -78,7 +74,7 @@ PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRFo >>>>>>> f3a673dcd37b27e654b0593c64e5893308437fe0 edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); - if (!srcConstituentWeights.label().empty()){ + if (!srcConstituentWeights.label().empty()) { input_constituent_weights_token_ = consumes>(srcConstituentWeights); } } From 7600128c05dc3450181d7f9fe3d30d885f48082e Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Wed, 15 Feb 2023 16:35:46 +0900 Subject: [PATCH 6/7] marker error removed --- RecoJets/JetProducers/plugins/PileupJetIdProducer.cc | 4 ---- 1 file changed, 4 deletions(-) diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc index 1d76b9b4d65fe..2e67c354622e9 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -69,10 +69,6 @@ PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRFo parameters_token_ = esConsumes(edm::ESInputTag("", globalCache->jec())); -<<<<<<< HEAD -======= - ->>>>>>> f3a673dcd37b27e654b0593c64e5893308437fe0 edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); if (!srcConstituentWeights.label().empty()) { input_constituent_weights_token_ = consumes>(srcConstituentWeights); From 283d569a9919bd9dd664f4d2006b6761b00fa987 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Wed, 15 Feb 2023 16:55:08 +0900 Subject: [PATCH 7/7] apply code-checks patch second time --- RecoJets/JetProducers/plugins/PileupJetIdProducer.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc index 2e67c354622e9..45fecafc022d4 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -68,7 +68,6 @@ PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRFo input_rho_token_ = consumes(iConfig.getParameter("rho")); parameters_token_ = esConsumes(edm::ESInputTag("", globalCache->jec())); - edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); if (!srcConstituentWeights.label().empty()) { input_constituent_weights_token_ = consumes>(srcConstituentWeights);