diff --git a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc index ff386726f54c4..573997a9bd90e 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,10 @@ 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 +514,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("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 c6bf82d5dab86..6224b1a085177 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 "", + useBugFix = True ) ) getattr(proc,jetTaskName).add(getattr(proc, puJetIdVarsCalculator)) diff --git a/RecoJets/JetProducers/interface/PileupJetIdAlgo.h b/RecoJets/JetProducers/interface/PileupJetIdAlgo.h index 3969257b30c4b..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); + 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 2edfeb24729d7..b810b09283bc9 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), + useBugFix_(iConfig.getParameter("useBugFix")) { 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->useBugFix()); + 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..5f630a09a4855 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 useBugFix() const {return useBugFix_; } private: std::vector vAlgoGBRForestsAndConstants_; @@ -74,6 +76,8 @@ class GBRForestsAndConstants { bool residualsFromTxt_; edm::FileInPath residualsTxt_; bool usePuppi_; + bool applyConstituentWeight_; + bool useBugFix_; }; 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..d09fdb6f7d54f 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(""), + 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 a5994a7b67fb4..a9c9cbe7e5534 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 useBugFix) { + // 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 (!useBugFix) { + 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 (!useBugFix && usePuppi && isPacked) { candPuppiWeight = lPack->puppiWeight(); - float candPt = (icand->pt()) * candPuppiWeight; + } + else if (useBugFix && applyConstituentWeight) { + candWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + float candPt = 0.0; + if (!useBugFix) + 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 (!useBugFix) { + 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 (!useBugFix) { + 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 (!useBugFix) + multNeut += candPuppiWeight; + else + multNeut += candWeight; } // EM candidated if (icand->pdgId() == 22) { - if (lLeadEm == nullptr || candPt > lLeadEm->pt()) { - lLeadEm = icand; + if (!useBugFix) { + 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 (!useBugFix) + 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 ((!useBugFix) && (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 ((useBugFix) && (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,89 @@ PileupJetIdentifier PileupJetIdAlgo::computeIdVariables(const reco::Jet* jet, } // trailing candidate - if (lTrail == nullptr || candPt < lTrail->pt()) { + if ((!useBugFix) && (lTrail == nullptr || candPt < lTrail->pt())) { lTrail = icand; } + else if ((useBugFix) && (lTrail == nullptr || candPt < (lTrail->pt())*TrailcandWeight)){ + lTrail = icand; + if (applyConstituentWeight) { + TrailcandWeight = constituentWeights[jet->sourceCandidatePtr(i)]; + } + } + + // average for pull variavble + if (useBugFix) { + 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; + } } + // // Finalize all variables assert(!(lLead == nullptr)); - - if (lSecond == nullptr) { - lSecond = lTrail; - } - if (lLeadNeut == nullptr) { - lLeadNeut = lTrail; - } - if (lLeadEm == nullptr) { - lLeadEm = lTrail; + if (!useBugFix) { + 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 +724,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 +737,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 (!useBugFix) { +// 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 (!useBugFix) { + 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 +817,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 (!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_); + 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 +888,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 (!useBugFix) { + 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.) {