From 526c244e5c530130d3c6d0e5794e02087e345eb3 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Mon, 13 Feb 2023 00:29:49 +0900 Subject: [PATCH 01/12] 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 | 207 +++++++++++------- 6 files changed, 163 insertions(+), 84 deletions(-) diff --git a/PhysicsTools/NanoAOD/python/custom_jme_cff.py b/PhysicsTools/NanoAOD/python/custom_jme_cff.py index c6bf82d5dab86..112e979bfc8fe 100644 --- a/PhysicsTools/NanoAOD/python/custom_jme_cff.py +++ b/PhysicsTools/NanoAOD/python/custom_jme_cff.py @@ -272,7 +272,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 a5994a7b67fb4..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(); @@ -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 9a95c252a64b647cfdf82619ee56a5121d243287 Mon Sep 17 00:00:00 2001 From: Nurfikri Norjoharuddeen Date: Mon, 13 Feb 2023 07:22:01 +0100 Subject: [PATCH 02/12] 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 214f9278579a92c42f720d8f78597546d8d56600 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Wed, 15 Feb 2023 15:26:38 +0900 Subject: [PATCH 03/12] 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 b9a69ca655f96d7bee4f144e067cb6f5932b4f41 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Wed, 15 Feb 2023 15:47:05 +0900 Subject: [PATCH 04/12] apply code-checks patch --- RecoJets/JetProducers/plugins/PileupJetIdProducer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc index fff45fccb91c9..2e67c354622e9 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -70,7 +70,7 @@ PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRFo edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); - if (!srcConstituentWeights.label().empty()){ + if (!srcConstituentWeights.label().empty()) { input_constituent_weights_token_ = consumes>(srcConstituentWeights); } } From 853463d7575633f828ffb6d415206767c3d90720 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Wed, 15 Feb 2023 16:55:08 +0900 Subject: [PATCH 05/12] 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); From 2b4106bb474dda48d391658cb680d3b7eb01285b Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Thu, 6 Apr 2023 12:03:05 +0900 Subject: [PATCH 06/12] Added applybuggy variable to optionalize using bug fixed version --- .../PileupAlgos/plugins/PuppiProducer.cc | 14 +- PhysicsTools/NanoAOD/python/custom_jme_cff.py | 4 +- .../JetProducers/interface/PileupJetIdAlgo.h | 2 +- .../plugins/PileupJetIdProducer.cc | 21 +- .../plugins/PileupJetIdProducer.h | 6 + .../JetProducers/python/PileupJetID_cfi.py | 2 + RecoJets/JetProducers/src/PileupJetIdAlgo.cc | 422 +++++++++++++----- 7 files changed, 357 insertions(+), 114 deletions(-) diff --git a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc index ff386726f54c4..e4f45595bc723 100644 --- a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc +++ b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc @@ -80,7 +80,9 @@ class PuppiProducer : public edm::stream::EDProducer<> { uint fNumOfPUVtxsForCharged; double fDZCutForChargedFromPUVtxs; bool fUseExistingWeights; + bool fApplyPhotonProtectionForExistingWeights; bool fClonePackedCands; + bool fapplybuggy; int fVtxNdofCut; double fVtxZCut; bool fUsePUProxyValue; @@ -104,6 +106,8 @@ PuppiProducer::PuppiProducer(const edm::ParameterSet& iConfig) { fNumOfPUVtxsForCharged = iConfig.getParameter("NumOfPUVtxsForCharged"); fDZCutForChargedFromPUVtxs = iConfig.getParameter("DeltaZCutForChargedFromPUVtxs"); fUseExistingWeights = iConfig.getParameter("useExistingWeights"); + fApplyPhotonProtectionForExistingWeights = iConfig.getParameter("applyPhotonProtectionForExistingWeights"); + fapplybuggy = iConfig.getParameter("applybuggy"); fClonePackedCands = iConfig.getParameter("clonePackedCands"); fVtxNdofCut = iConfig.getParameter("vtxNdofCut"); fVtxZCut = iConfig.getParameter("vtxZCut"); @@ -365,7 +369,13 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { } // 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)) + (lPack->pt() > fPtMaxPhotons) && fapplybuggy) + curpupweight = 1; + + + // 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) && fapplybuggy) curpupweight = 1; lWeights.push_back(curpupweight); lPackCtr++; @@ -511,6 +521,8 @@ 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("applybuggy",false); desc.add("clonePackedCands", false); desc.add("vtxNdofCut", 4); desc.add("vtxZCut", 24); diff --git a/PhysicsTools/NanoAOD/python/custom_jme_cff.py b/PhysicsTools/NanoAOD/python/custom_jme_cff.py index c6bf82d5dab86..f4bdfbd91a540 100644 --- a/PhysicsTools/NanoAOD/python/custom_jme_cff.py +++ b/PhysicsTools/NanoAOD/python/custom_jme_cff.py @@ -272,7 +272,9 @@ def AddPileUpJetIDVars(proc, jetName="", jetSrc="", jetTableName="", jetTaskName vertexes = "offlineSlimmedPrimaryVertices", inputIsCorrected = True, applyJec = False, - usePuppi = True if "PUPPI" in jetName.upper() else False + usePuppi = True if "PUPPI" in jetName.upper() else False, + srcConstituentWeights = "packedPFCandidatespuppi" if "PUPPI" in jetName.upper() else "", + applybuggy = False ) ) getattr(proc,jetTaskName).add(getattr(proc, puJetIdVarsCalculator)) diff --git a/RecoJets/JetProducers/interface/PileupJetIdAlgo.h b/RecoJets/JetProducers/interface/PileupJetIdAlgo.h index 3969257b30c4b..06293d0a453fc 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, bool usePuppi, edm::ValueMap& constituentWeights, bool applyConstituentWeight, bool applybuggy); 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..55ad8ffaf2703 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -27,7 +27,9 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig) applyJec_(iConfig.getParameter("applyJec")), jec_(iConfig.getParameter("jec")), residualsFromTxt_(iConfig.getParameter("residualsFromTxt")), - usePuppi_(iConfig.getParameter("usePuppi")) { + usePuppi_(iConfig.getParameter("usePuppi")), + applyConstituentWeight_(false), + applybuggy_(iConfig.getParameter("applybuggy")) { if (residualsFromTxt_) { residualsTxt_ = iConfig.getParameter("residualsTxt"); } @@ -40,6 +42,11 @@ 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 +69,10 @@ 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 srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); + if (!srcConstituentWeights.label().empty()){ + input_constituent_weights_token_ = consumes>(srcConstituentWeights); + } } // ------------------------------------------------------------------------------------------ @@ -79,6 +90,11 @@ void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe Handle> jetHandle; iEvent.getByToken(input_jet_token_, jetHandle); const View& jets = *jetHandle; + edm::ValueMap constituentWeights; + if (!input_constituent_weights_token_.isUninitialized()) { + constituentWeights = iEvent.get(input_constituent_weights_token_); + } + // input variables Handle> vmap; @@ -167,7 +183,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()); + puIdentifier = ialgo->computeIdVariables(theJet, jec, &(*vtx), *vertexes, rho, gc->usePuppi(), constituentWeights, gc->applyConstituentWeight(), gc->applybuggy()); + 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..b18595a6e2fff 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.h +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.h @@ -62,6 +62,8 @@ class GBRForestsAndConstants { bool residualsFromTxt() const { return residualsFromTxt_; } edm::FileInPath const& residualsTxt() const { return residualsTxt_; } bool usePuppi() const { return usePuppi_; } + bool applyConstituentWeight() const { return applyConstituentWeight_; } + bool applybuggy() const {return applybuggy_; } private: std::vector vAlgoGBRForestsAndConstants_; @@ -74,6 +76,8 @@ class GBRForestsAndConstants { bool residualsFromTxt_; edm::FileInPath residualsTxt_; bool usePuppi_; + bool applyConstituentWeight_; + bool applybuggy_; }; class PileupJetIdProducer : public edm::stream::EDProducer> { @@ -99,6 +103,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..84537d724275e 100644 --- a/RecoJets/JetProducers/python/PileupJetID_cfi.py +++ b/RecoJets/JetProducers/python/PileupJetID_cfi.py @@ -33,6 +33,8 @@ inputIsCorrected = cms.bool(False), residualsFromTxt = cms.bool(False), usePuppi = cms.bool(False), + srcConstituentWeights = cms.InputTag(""), + applybuggy = cms.bool(False) # 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 a5994a7b67fb4..20a098980b779 100644 --- a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc +++ b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc @@ -263,7 +263,11 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, const reco::Vertex* vtx, const reco::VertexCollection& allvtx, double rho, - bool usePuppi) { + bool usePuppi, + edm::ValueMap& constituentWeights, + bool applyConstituentWeight, + bool applybuggy) { + // initialize all variables to 0 resetVariables(); @@ -315,6 +319,10 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, 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 sumW2(0.0); + float sum_deta(0.0), sum_dphi(0.0); + float ave_deta(0.0), ave_dphi(0.0); float multNeut = 0.0; setPtEtaPhi( *jet, internalId_.jetPt_, internalId_.jetEta_, internalId_.jetPhi_); // use corrected pt for jet kinematics @@ -323,8 +331,20 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.rho_ = rho; float dRmin(1000); - - for (unsigned i = 0; i < jet->numberOfSourceCandidatePtrs(); ++i) { + float LeadcandWeight = 1.0; + float SecondcandWeight = 1.0; + float LeadNeutcandWeight = 1.0; + float LeadEmcandWeight = 1.0; + float LeadChcandWeight = 1.0; + float TrailcandWeight = 1.0; + unsigned nCandPtrs(0); + if (applybuggy) { + nCandPtrs = jet->numberOfDaughters(); + } + else { + nCandPtrs = jet->numberOfSourceCandidatePtrs(); + } + for (unsigned i = 0; i < nCandPtrs; ++i) { reco::CandidatePtr pfJetConstituent = jet->sourceCandidatePtr(i); const reco::Candidate* icand = pfJetConstituent.get(); @@ -335,9 +355,18 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, isPacked = false; } float candPuppiWeight = 1.0; - if (usePuppi && isPacked) + float candWeight = 1.0; + if (applybuggy && usePuppi && isPacked) { candPuppiWeight = lPack->puppiWeight(); - float candPt = (icand->pt()) * candPuppiWeight; + } + else if (!applybuggy && applyConstituentWeight) { + candWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + float candPt = 0.0; + if (applybuggy) + candPt = (icand->pt()) * candPuppiWeight; + else + candPt = (icand->pt()) * candWeight; float candPtFrac = candPt / jetPt; float candDr = reco::deltaR(*icand, *jet); float candDeta = icand->eta() - jet->eta(); @@ -349,13 +378,29 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, dRmin = candDr; // // all particles - if (lLead == nullptr || candPt > lLead->pt()) { - lSecond = lLead; - lLead = icand; - } else if ((lSecond == nullptr || candPt > lSecond->pt()) && (candPt < lLead->pt())) { - lSecond = icand; + if (applybuggy) { + if (lLead == nullptr || candPt > lLead->pt()) { + lSecond = lLead; + lLead = icand; + } else if ((lSecond == nullptr || candPt > lSecond->pt()) && (candPt < lLead->pt())) { + lSecond = icand; + } + } + else { + 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)) { + lSecond = icand; + if (applyConstituentWeight) { + SecondcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + } } - // // average shapes internalId_.dRMean_ += candPtDr; internalId_.dR2Mean_ += candPtDr * candPtDr; @@ -375,8 +420,18 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, // neutrals if (abs(icand->pdgId()) == 130) { - if (lLeadNeut == nullptr || candPt > lLeadNeut->pt()) { - lLeadNeut = icand; + if (applybuggy) { + if (lLeadNeut == nullptr || candPt > lLeadNeut->pt()) { + lLeadNeut = icand; + } + } + else { + if (lLeadNeut == nullptr || candPt > (lLeadNeut->pt())*LeadNeutcandWeight) { + lLeadNeut = icand; + if (applyConstituentWeight) { + LeadNeutcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + } } internalId_.dRMeanNeut_ += candPtDr; fracNeut.push_back(candPtFrac); @@ -385,55 +440,109 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } internalId_.ptDNe_ += candPt * candPt; sumPtNe += candPt; - multNeut += candPuppiWeight; + if (applybuggy) + multNeut += candPuppiWeight; + else + multNeut += candWeight; } // EM candidated if (icand->pdgId() == 22) { - if (lLeadEm == nullptr || candPt > lLeadEm->pt()) { - lLeadEm = icand; + if (applybuggy) { + if (lLeadEm == nullptr || candPt > lLeadEm->pt()) { + lLeadEm = icand; + } + internalId_.dRMeanEm_ += candPtDr; + fracEm.push_back(candPtFrac); + if (icone < ncones) { + *coneEmFracs[icone] += candPt; + } + internalId_.ptDNe_ += candPt * candPt; + sumPtNe += candPt; + multNeut += candPuppiWeight; } - internalId_.dRMeanEm_ += candPtDr; - fracEm.push_back(candPtFrac); - if (icone < ncones) { - *coneEmFracs[icone] += candPt; + else { + if (lLeadEm == nullptr || candPt > (lLeadEm->pt())*LeadEmcandWeight) { + lLeadEm = icand; + if (applyConstituentWeight) { + LeadEmcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + } + internalId_.dRMeanEm_ += candPtDr; + fracEm.push_back(candPtFrac); + if (icone < ncones) { + *coneEmFracs[icone] += candPt; + } + internalId_.ptDNe_ += candPt * candPt; + sumPtNe += candPt; + multNeut += candWeight; } - internalId_.ptDNe_ += candPt * candPt; - sumPtNe += candPt; - multNeut += candPuppiWeight; } - if ((abs(icand->pdgId()) == 1) || (abs(icand->pdgId()) == 2)) - multNeut += candPuppiWeight; - + if (abs(icand->pdgId()) == 1 || abs(icand->pdgId()) == 2) { + if (applybuggy) + multNeut += candPuppiWeight; + else + multNeut += candWeight; + } // Charged particles if (icand->charge() != 0) { - if (lLeadCh == nullptr || candPt > lLeadCh->pt()) { - lLeadCh = icand; - - const reco::Track* pfTrk = icand->bestTrack(); - if (lPF && std::abs(icand->pdgId()) == 13 && pfTrk == nullptr) { - reco::MuonRef lmuRef = lPF->muonRef(); - if (lmuRef.isNonnull()) { - const reco::Muon& lmu = *lmuRef.get(); - pfTrk = lmu.bestTrack(); - edm::LogWarning("BadMuon") - << "Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack "; + if ((applybuggy) && (lLeadCh == nullptr || candPt > lLeadCh->pt())) { + lLeadCh = icand; + + const reco::Track* pfTrk = icand->bestTrack(); + if (lPF && std::abs(icand->pdgId()) == 13 && pfTrk == nullptr) { + reco::MuonRef lmuRef = lPF->muonRef(); + if (lmuRef.isNonnull()) { + const reco::Muon& lmu = *lmuRef.get(); + pfTrk = lmu.bestTrack(); + edm::LogWarning("BadMuon") + << "Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack "; + } + } + if (pfTrk == nullptr) { //protection against empty pointers for the miniAOD case + //To handle the electron case + if (isPacked) { + internalId_.d0_ = std::abs(lPack->dxy(vtx->position())); + internalId_.dZ_ = std::abs(lPack->dz(vtx->position())); + } else if (lPF != nullptr) { + pfTrk = (lPF->trackRef().get() == nullptr) ? lPF->gsfTrackRef().get() : lPF->trackRef().get(); + internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position())); + internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position())); + } + } else { + internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position())); + internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position())); } } - if (pfTrk == nullptr) { //protection against empty pointers for the miniAOD case - //To handle the electron case - if (isPacked) { - internalId_.d0_ = std::abs(lPack->dxy(vtx->position())); - internalId_.dZ_ = std::abs(lPack->dz(vtx->position())); - } else if (lPF != nullptr) { - pfTrk = (lPF->trackRef().get() == nullptr) ? lPF->gsfTrackRef().get() : lPF->trackRef().get(); + else if ((!applybuggy) && (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(); + if (lmuRef.isNonnull()) { + const reco::Muon& lmu = *lmuRef.get(); + pfTrk = lmu.bestTrack(); + edm::LogWarning("BadMuon") + << "Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack "; + } + } + if (pfTrk == nullptr) { //protection against empty pointers for the miniAOD case + //To handle the electron case + if (isPacked) { + internalId_.d0_ = std::abs(lPack->dxy(vtx->position())); + internalId_.dZ_ = std::abs(lPack->dz(vtx->position())); + } else if (lPF != nullptr) { + pfTrk = (lPF->trackRef().get() == nullptr) ? lPF->gsfTrackRef().get() : lPF->trackRef().get(); + internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position())); + internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position())); + } + } else { internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position())); internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position())); } - } else { - internalId_.d0_ = std::abs(pfTrk->dxy(vtx->position())); - internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position())); } - } internalId_.dRMeanCh_ += candPtDr; internalId_.ptDCh_ += candPt * candPt; fracCh.push_back(candPtFrac); @@ -525,27 +634,93 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } // trailing candidate - if (lTrail == nullptr || candPt < lTrail->pt()) { + if ((applybuggy) && (lTrail == nullptr || candPt < lTrail->pt())) { lTrail = icand; } + else if ((!applybuggy) && (lTrail == nullptr || candPt < (lTrail->pt())*TrailcandWeight)){ + lTrail = icand; + if (applyConstituentWeight) { + TrailcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + } + + // average for pull variavble + if (!applybuggy) { + 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)); - - if (lSecond == nullptr) { - lSecond = lTrail; - } - if (lLeadNeut == nullptr) { - lLeadNeut = lTrail; - } - if (lLeadEm == nullptr) { - lLeadEm = lTrail; + if (applybuggy) { + if (lSecond == nullptr) { + lSecond = lTrail; + } + if (lLeadNeut == nullptr) { + lLeadNeut = lTrail; + } + if (lLeadEm == nullptr) { + lLeadEm = lTrail; + } + if (lLeadCh == nullptr) { + lLeadCh = lTrail; + } } - if (lLeadCh == nullptr) { - lLeadCh = lTrail; + else { + 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){ + 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){ + 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){ + 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 +728,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 (usePuppi || applyConstituentWeight) internalId_.nNeutrals_ = multNeut; } else { internalId_.nCharged_ = pfjet->chargedMultiplicity(); @@ -566,52 +741,72 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, 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(); - } + if (applybuggy) { +// 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 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++) { + + for (size_t i = 0; i < nCandPtrs; i++) { const auto& part = jet->daughterPtr(i); - if (!(part.isAvailable() && part.isNonnull())) { - continue; - } + const reco::CandidatePtr temp_pfJetConsituent = jet->sourceCandidatePtr(i); + const reco::Candidate* fix_part = temp_pfJetConsituent.get(); + float deta(0.0); + float dphi(0.0); + float weight(0.0); + if (applybuggy) { + 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 partPuppiWeight = 1.0; + if (usePuppi) { + const pat::PackedCandidate* partpack = dynamic_cast(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); + weight = partPuppiWeight * (part->pt()) * partPuppiWeight * (part->pt()); + deta = part->eta() - jet->eta(); + dphi = reco::deltaPhi(*part, *jet); + } + else { + float candWeight = 1.0; + if (applyConstituentWeight) { + candWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + weight = candWeight * (part->pt()) * candWeight * (part->pt()); + deta = fix_part->eta() - jet->eta(); + dphi = reco::deltaPhi(*fix_part, *jet); + } float ddeta, ddphi, ddR; ddeta = deta - ave_deta; ddphi = dphi - ave_dphi; @@ -626,13 +821,13 @@ 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_); - + if (applybuggy) { + 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()); std::sort(fracEm.begin(), fracEm.end(), std::greater()); @@ -697,9 +892,18 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.sumPt_ = sumPt; internalId_.sumChPt_ = sumPtCh; internalId_.sumNePt_ = sumPtNe; - - internalId_.jetR_ = lLead->pt() / sumPt; - internalId_.jetRchg_ = lLeadCh->pt() / sumPt; + if (applybuggy) { + internalId_.jetR_ = lLead->pt() / sumPt; + internalId_.jetRchg_ = lLeadCh->pt() / sumPt; + } + else { + 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 27cfac46894ccdd77d8cc017511cda44f1974a69 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Fri, 14 Apr 2023 00:17:59 +0900 Subject: [PATCH 07/12] change name applybuggy to useBugFix and fApplyPhotonProtectionForExistingWeights flag double counting --- .../PileupAlgos/plugins/PuppiProducer.cc | 10 ++--- PhysicsTools/NanoAOD/python/custom_jme_cff.py | 2 +- .../JetProducers/interface/PileupJetIdAlgo.h | 2 +- .../plugins/PileupJetIdProducer.cc | 4 +- .../plugins/PileupJetIdProducer.h | 4 +- .../JetProducers/python/PileupJetID_cfi.py | 2 +- RecoJets/JetProducers/src/PileupJetIdAlgo.cc | 40 +++++++++---------- 7 files changed, 32 insertions(+), 32 deletions(-) diff --git a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc index e4f45595bc723..2241a2cf5fba4 100644 --- a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc +++ b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc @@ -82,7 +82,7 @@ class PuppiProducer : public edm::stream::EDProducer<> { bool fUseExistingWeights; bool fApplyPhotonProtectionForExistingWeights; bool fClonePackedCands; - bool fapplybuggy; + bool fuseBugFix; int fVtxNdofCut; double fVtxZCut; bool fUsePUProxyValue; @@ -107,7 +107,7 @@ PuppiProducer::PuppiProducer(const edm::ParameterSet& iConfig) { fDZCutForChargedFromPUVtxs = iConfig.getParameter("DeltaZCutForChargedFromPUVtxs"); fUseExistingWeights = iConfig.getParameter("useExistingWeights"); fApplyPhotonProtectionForExistingWeights = iConfig.getParameter("applyPhotonProtectionForExistingWeights"); - fapplybuggy = iConfig.getParameter("applybuggy"); + fuseBugFix = iConfig.getParameter("useBugFix"); fClonePackedCands = iConfig.getParameter("clonePackedCands"); fVtxNdofCut = iConfig.getParameter("vtxNdofCut"); fVtxZCut = iConfig.getParameter("vtxZCut"); @@ -369,13 +369,13 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { } // 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) && fapplybuggy) + (lPack->pt() > fPtMaxPhotons) && fuseBugFix) curpupweight = 1; // 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) && fapplybuggy) + (std::abs(lPack->eta()) < fEtaMaxPhotons) && (lPack->pt() > fPtMaxPhotons)) curpupweight = 1; lWeights.push_back(curpupweight); lPackCtr++; @@ -522,7 +522,7 @@ void PuppiProducer::fillDescriptions(edm::ConfigurationDescriptions& description desc.add("DeltaZCutForChargedFromPUVtxs", 0.2); desc.add("useExistingWeights", false); desc.add("applyPhotonProtectionForExistingWeights", false); - desc.add("applybuggy",false); + desc.add("useBugFix",false); desc.add("clonePackedCands", false); desc.add("vtxNdofCut", 4); desc.add("vtxZCut", 24); diff --git a/PhysicsTools/NanoAOD/python/custom_jme_cff.py b/PhysicsTools/NanoAOD/python/custom_jme_cff.py index f4bdfbd91a540..e2b75f40fed9b 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 applyJec = False, usePuppi = True if "PUPPI" in jetName.upper() else False, srcConstituentWeights = "packedPFCandidatespuppi" if "PUPPI" in jetName.upper() else "", - applybuggy = False + useBugFix = False ) ) getattr(proc,jetTaskName).add(getattr(proc, puJetIdVarsCalculator)) diff --git a/RecoJets/JetProducers/interface/PileupJetIdAlgo.h b/RecoJets/JetProducers/interface/PileupJetIdAlgo.h index 06293d0a453fc..d6349eeb61a69 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, edm::ValueMap& constituentWeights, bool applyConstituentWeight, bool applybuggy); + const reco::Jet* jet, float jec, const reco::Vertex*, const reco::VertexCollection&, double rho, bool usePuppi, edm::ValueMap& constituentWeights, bool applyConstituentWeight, bool useBugFix); 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 55ad8ffaf2703..b810b09283bc9 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -29,7 +29,7 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig) residualsFromTxt_(iConfig.getParameter("residualsFromTxt")), usePuppi_(iConfig.getParameter("usePuppi")), applyConstituentWeight_(false), - applybuggy_(iConfig.getParameter("applybuggy")) { + useBugFix_(iConfig.getParameter("useBugFix")) { if (residualsFromTxt_) { residualsTxt_ = iConfig.getParameter("residualsTxt"); } @@ -183,7 +183,7 @@ 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(), constituentWeights, gc->applyConstituentWeight(), gc->applybuggy()); + puIdentifier = ialgo->computeIdVariables(theJet, jec, &(*vtx), *vertexes, rho, gc->usePuppi(), constituentWeights, gc->applyConstituentWeight(), gc->useBugFix()); ids.push_back(puIdentifier); } else { diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.h b/RecoJets/JetProducers/plugins/PileupJetIdProducer.h index b18595a6e2fff..5f630a09a4855 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.h +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.h @@ -63,7 +63,7 @@ class GBRForestsAndConstants { edm::FileInPath const& residualsTxt() const { return residualsTxt_; } bool usePuppi() const { return usePuppi_; } bool applyConstituentWeight() const { return applyConstituentWeight_; } - bool applybuggy() const {return applybuggy_; } + bool useBugFix() const {return useBugFix_; } private: std::vector vAlgoGBRForestsAndConstants_; @@ -77,7 +77,7 @@ class GBRForestsAndConstants { edm::FileInPath residualsTxt_; bool usePuppi_; bool applyConstituentWeight_; - bool applybuggy_; + bool useBugFix_; }; class PileupJetIdProducer : public edm::stream::EDProducer> { diff --git a/RecoJets/JetProducers/python/PileupJetID_cfi.py b/RecoJets/JetProducers/python/PileupJetID_cfi.py index 84537d724275e..d09fdb6f7d54f 100644 --- a/RecoJets/JetProducers/python/PileupJetID_cfi.py +++ b/RecoJets/JetProducers/python/PileupJetID_cfi.py @@ -34,7 +34,7 @@ residualsFromTxt = cms.bool(False), usePuppi = cms.bool(False), srcConstituentWeights = cms.InputTag(""), - applybuggy = cms.bool(False) + useBugFix = cms.bool(False) # 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 20a098980b779..714115fdcee3d 100644 --- a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc +++ b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc @@ -266,7 +266,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, bool usePuppi, edm::ValueMap& constituentWeights, bool applyConstituentWeight, - bool applybuggy) { + bool useBugFix) { // initialize all variables to 0 resetVariables(); @@ -338,7 +338,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, float LeadChcandWeight = 1.0; float TrailcandWeight = 1.0; unsigned nCandPtrs(0); - if (applybuggy) { + if (!useBugFix) { nCandPtrs = jet->numberOfDaughters(); } else { @@ -356,14 +356,14 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } float candPuppiWeight = 1.0; float candWeight = 1.0; - if (applybuggy && usePuppi && isPacked) { + if (!useBugFix && usePuppi && isPacked) { candPuppiWeight = lPack->puppiWeight(); } - else if (!applybuggy && applyConstituentWeight) { + else if (useBugFix && applyConstituentWeight) { candWeight = constituentWeights[jet->sourceCandidatePtr(i)]; } float candPt = 0.0; - if (applybuggy) + if (!useBugFix) candPt = (icand->pt()) * candPuppiWeight; else candPt = (icand->pt()) * candWeight; @@ -378,7 +378,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, dRmin = candDr; // // all particles - if (applybuggy) { + if (!useBugFix) { if (lLead == nullptr || candPt > lLead->pt()) { lSecond = lLead; lLead = icand; @@ -420,7 +420,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, // neutrals if (abs(icand->pdgId()) == 130) { - if (applybuggy) { + if (!useBugFix) { if (lLeadNeut == nullptr || candPt > lLeadNeut->pt()) { lLeadNeut = icand; } @@ -440,14 +440,14 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } internalId_.ptDNe_ += candPt * candPt; sumPtNe += candPt; - if (applybuggy) + if (!useBugFix) multNeut += candPuppiWeight; else multNeut += candWeight; } // EM candidated if (icand->pdgId() == 22) { - if (applybuggy) { + if (!useBugFix) { if (lLeadEm == nullptr || candPt > lLeadEm->pt()) { lLeadEm = icand; } @@ -478,14 +478,14 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } } if (abs(icand->pdgId()) == 1 || abs(icand->pdgId()) == 2) { - if (applybuggy) + if (!useBugFix) multNeut += candPuppiWeight; else multNeut += candWeight; } // Charged particles if (icand->charge() != 0) { - if ((applybuggy) && (lLeadCh == nullptr || candPt > lLeadCh->pt())) { + if ((!useBugFix) && (lLeadCh == nullptr || candPt > lLeadCh->pt())) { lLeadCh = icand; const reco::Track* pfTrk = icand->bestTrack(); @@ -513,7 +513,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.dZ_ = std::abs(pfTrk->dz(vtx->position())); } } - else if ((!applybuggy) && (lLeadCh == nullptr || candPt > (lLeadCh->pt()) * LeadChcandWeight)) { + else if ((useBugFix) && (lLeadCh == nullptr || candPt > (lLeadCh->pt()) * LeadChcandWeight)) { lLeadCh = icand; if (applyConstituentWeight) { LeadChcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; @@ -634,10 +634,10 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } // trailing candidate - if ((applybuggy) && (lTrail == nullptr || candPt < lTrail->pt())) { + if ((!useBugFix) && (lTrail == nullptr || candPt < lTrail->pt())) { lTrail = icand; } - else if ((!applybuggy) && (lTrail == nullptr || candPt < (lTrail->pt())*TrailcandWeight)){ + else if ((useBugFix) && (lTrail == nullptr || candPt < (lTrail->pt())*TrailcandWeight)){ lTrail = icand; if (applyConstituentWeight) { TrailcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; @@ -645,7 +645,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } // average for pull variavble - if (!applybuggy) { + if (useBugFix) { float weight2 = candPt * candPt; sumW2 += weight2; float deta = icand->eta() - jet->eta(); @@ -662,7 +662,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, // // Finalize all variables assert(!(lLead == nullptr)); - if (applybuggy) { + if (!useBugFix) { if (lSecond == nullptr) { lSecond = lTrail; } @@ -742,7 +742,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, ///////////////////////pull variable/////////////////////////////////// - if (applybuggy) { + if (!useBugFix) { // float sumW2(0.0); // float sum_deta(0.0), sum_dphi(0.0); // float ave_deta(0.0), ave_dphi(0.0); @@ -782,7 +782,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, float deta(0.0); float dphi(0.0); float weight(0.0); - if (applybuggy) { + if (!useBugFix) { if (!(part.isAvailable() && part.isNonnull())) { continue; } @@ -821,7 +821,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } internalId_.pull_ = pull_tmp; /////////////////////////////////////////////////////////////////////// - if (applybuggy) { + if (!useBugFix) { setPtEtaPhi(*lLead, internalId_.leadPt_, internalId_.leadEta_, internalId_.leadPhi_); setPtEtaPhi(*lSecond, internalId_.secondPt_, internalId_.secondEta_, internalId_.secondPhi_); setPtEtaPhi(*lLeadNeut, internalId_.leadNeutPt_, internalId_.leadNeutEta_, internalId_.leadNeutPhi_); @@ -892,7 +892,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, internalId_.sumPt_ = sumPt; internalId_.sumChPt_ = sumPtCh; internalId_.sumNePt_ = sumPtNe; - if (applybuggy) { + if (!useBugFix) { internalId_.jetR_ = lLead->pt() / sumPt; internalId_.jetRchg_ = lLeadCh->pt() / sumPt; } From 4ba3efa1a37f3a5e1b49032b5363ac0ffa18f099 Mon Sep 17 00:00:00 2001 From: jshin96 <87310162+jshin96@users.noreply.github.com> Date: Mon, 24 Apr 2023 09:52:42 +0900 Subject: [PATCH 08/12] Update PileupJetIdProducer.cc deleted extra lines that were duplicate --- RecoJets/JetProducers/plugins/PileupJetIdProducer.cc | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc index b14fd9a594dc8..e000d9eb2306a 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -91,12 +91,7 @@ void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe Handle> jetHandle; iEvent.getByToken(input_jet_token_, jetHandle); const View& jets = *jetHandle; - edm::ValueMap constituentWeights; - if (!input_constituent_weights_token_.isUninitialized()) { - constituentWeights = iEvent.get(input_constituent_weights_token_); - } - - + // Constituent weight (e.g PUPPI) Value Map edm::ValueMap constituentWeights; if (!input_constituent_weights_token_.isUninitialized()) { From 1bc0651e66171cfca52062992022692b2953e0bc Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Mon, 8 May 2023 11:55:58 +0900 Subject: [PATCH 09/12] Fix bug and made bug fixing optional --- .../PileupAlgos/plugins/PuppiProducer.cc | 1 + .../JetProducers/interface/PileupJetIdAlgo.h | 9 ++--- .../plugins/PileupJetIdProducer.cc | 6 ++-- .../plugins/PileupJetIdProducer.h | 1 - .../JetProducers/python/PileupJetID_cfi.py | 1 - RecoJets/JetProducers/src/PileupJetIdAlgo.cc | 34 ++++--------------- 6 files changed, 11 insertions(+), 41 deletions(-) diff --git a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc index 68979e62c149c..2241a2cf5fba4 100644 --- a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc +++ b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc @@ -372,6 +372,7 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { (lPack->pt() > fPtMaxPhotons) && fuseBugFix) curpupweight = 1; + // 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)) diff --git a/RecoJets/JetProducers/interface/PileupJetIdAlgo.h b/RecoJets/JetProducers/interface/PileupJetIdAlgo.h index aecda1c1159d7..d6349eeb61a69 100644 --- a/RecoJets/JetProducers/interface/PileupJetIdAlgo.h +++ b/RecoJets/JetProducers/interface/PileupJetIdAlgo.h @@ -28,13 +28,8 @@ 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, bool usePuppi, edm::ValueMap& constituentWeights, bool applyConstituentWeight, bool useBugFix); 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 e000d9eb2306a..b810b09283bc9 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.cc @@ -30,7 +30,6 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig) usePuppi_(iConfig.getParameter("usePuppi")), applyConstituentWeight_(false), useBugFix_(iConfig.getParameter("useBugFix")) { - if (residualsFromTxt_) { residualsTxt_ = iConfig.getParameter("residualsTxt"); } @@ -43,11 +42,11 @@ GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig) if (!runMvas_) { assert(algos.size() == 1); } - edm::InputTag srcConstituentWeights = iConfig.getParameter("srcConstituentWeights"); if (!srcConstituentWeights.label().empty()){ applyConstituentWeight_ = true; } + } // ------------------------------------------------------------------------------------------ @@ -91,13 +90,12 @@ void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe Handle> jetHandle; 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()) { diff --git a/RecoJets/JetProducers/plugins/PileupJetIdProducer.h b/RecoJets/JetProducers/plugins/PileupJetIdProducer.h index de3c054795bfb..5f630a09a4855 100644 --- a/RecoJets/JetProducers/plugins/PileupJetIdProducer.h +++ b/RecoJets/JetProducers/plugins/PileupJetIdProducer.h @@ -65,7 +65,6 @@ class GBRForestsAndConstants { bool applyConstituentWeight() const { return applyConstituentWeight_; } bool useBugFix() const {return useBugFix_; } - private: std::vector vAlgoGBRForestsAndConstants_; diff --git a/RecoJets/JetProducers/python/PileupJetID_cfi.py b/RecoJets/JetProducers/python/PileupJetID_cfi.py index e99088f5fa59d..d09fdb6f7d54f 100644 --- a/RecoJets/JetProducers/python/PileupJetID_cfi.py +++ b/RecoJets/JetProducers/python/PileupJetID_cfi.py @@ -35,7 +35,6 @@ usePuppi = cms.bool(False), srcConstituentWeights = cms.InputTag(""), useBugFix = cms.bool(False) - # 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 b4f4c79e70334..714115fdcee3d 100644 --- a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc +++ b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc @@ -267,6 +267,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, edm::ValueMap& constituentWeights, bool applyConstituentWeight, bool useBugFix) { + // initialize all variables to 0 resetVariables(); @@ -283,7 +284,6 @@ 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); @@ -324,9 +324,6 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, float sum_deta(0.0), sum_dphi(0.0); float ave_deta(0.0), ave_dphi(0.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(); @@ -349,6 +346,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } for (unsigned i = 0; i < nCandPtrs; ++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); @@ -378,6 +376,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, if (candDr < dRmin) dRmin = candDr; + // // all particles if (!useBugFix) { if (lLead == nullptr || candPt > lLead->pt()) { @@ -418,7 +417,8 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, if (icone < ncones) { *coneFracs[icone] += candPt; } - // neutrals Neutral hadrons + + // neutrals if (abs(icand->pdgId()) == 130) { if (!useBugFix) { if (lLeadNeut == nullptr || candPt > lLeadNeut->pt()) { @@ -433,7 +433,6 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } } } - internalId_.dRMeanNeut_ += candPtDr; fracNeut.push_back(candPtFrac); if (icone < ncones) { @@ -446,8 +445,7 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, else multNeut += candWeight; } - - // EM candidated photon + // EM candidated if (icand->pdgId() == 22) { if (!useBugFix) { if (lLeadEm == nullptr || candPt > lLeadEm->pt()) { @@ -638,22 +636,6 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, // trailing candidate if ((!useBugFix) && (lTrail == nullptr || candPt < lTrail->pt())) { 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; } else if ((useBugFix) && (lTrail == nullptr || candPt < (lTrail->pt())*TrailcandWeight)){ lTrail = icand; @@ -679,8 +661,6 @@ 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)); if (!useBugFix) { if (lSecond == nullptr) { @@ -757,8 +737,6 @@ 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(); From 3b952b1c2fc43ee0ed5f1ba8be64d0ce4e312673 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Fri, 12 May 2023 23:17:46 +0900 Subject: [PATCH 10/12] changed useBugFix to True --- PhysicsTools/NanoAOD/python/custom_jme_cff.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PhysicsTools/NanoAOD/python/custom_jme_cff.py b/PhysicsTools/NanoAOD/python/custom_jme_cff.py index e2b75f40fed9b..6224b1a085177 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 applyJec = False, usePuppi = True if "PUPPI" in jetName.upper() else False, srcConstituentWeights = "packedPFCandidatespuppi" if "PUPPI" in jetName.upper() else "", - useBugFix = False + useBugFix = True ) ) getattr(proc,jetTaskName).add(getattr(proc, puJetIdVarsCalculator)) From 50cea296015910453429fe71fe43d016bb63ffc5 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Sat, 3 Jun 2023 00:36:37 +0900 Subject: [PATCH 11/12] move average calculation out of the loop and remove redundant useBugFix in PuppiProducer --- CommonTools/PileupAlgos/plugins/PuppiProducer.cc | 7 ------- RecoJets/JetProducers/src/PileupJetIdAlgo.cc | 6 +++--- 2 files changed, 3 insertions(+), 10 deletions(-) diff --git a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc index 2241a2cf5fba4..573997a9bd90e 100644 --- a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc +++ b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc @@ -82,7 +82,6 @@ class PuppiProducer : public edm::stream::EDProducer<> { bool fUseExistingWeights; bool fApplyPhotonProtectionForExistingWeights; bool fClonePackedCands; - bool fuseBugFix; int fVtxNdofCut; double fVtxZCut; bool fUsePUProxyValue; @@ -107,7 +106,6 @@ PuppiProducer::PuppiProducer(const edm::ParameterSet& iConfig) { fDZCutForChargedFromPUVtxs = iConfig.getParameter("DeltaZCutForChargedFromPUVtxs"); fUseExistingWeights = iConfig.getParameter("useExistingWeights"); fApplyPhotonProtectionForExistingWeights = iConfig.getParameter("applyPhotonProtectionForExistingWeights"); - fuseBugFix = iConfig.getParameter("useBugFix"); fClonePackedCands = iConfig.getParameter("clonePackedCands"); fVtxNdofCut = iConfig.getParameter("vtxNdofCut"); fVtxZCut = iConfig.getParameter("vtxZCut"); @@ -367,11 +365,6 @@ 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) && fuseBugFix) - curpupweight = 1; - // Optional: Protect high pT photons (important for gamma to hadronic recoil balance) for existing weights. if (fApplyPhotonProtectionForExistingWeights && (fPtMaxPhotons > 0) && (lPack->pdgId() == 22) && diff --git a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc index 714115fdcee3d..91a4dcd099267 100644 --- a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc +++ b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc @@ -767,10 +767,10 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, 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; } + 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); From 2d370bf37ecafa7872a4532fba87ecdbda2cbbb8 Mon Sep 17 00:00:00 2001 From: Jihoon Shin Date: Sat, 3 Jun 2023 00:46:33 +0900 Subject: [PATCH 12/12] move average calculation out of the loop --- RecoJets/JetProducers/src/PileupJetIdAlgo.cc | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc index 91a4dcd099267..a9c9cbe7e5534 100644 --- a/RecoJets/JetProducers/src/PileupJetIdAlgo.cc +++ b/RecoJets/JetProducers/src/PileupJetIdAlgo.cc @@ -652,10 +652,6 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, 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; - } } } @@ -767,12 +763,12 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, 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; } } + 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 < nCandPtrs; i++) {