Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

improving performance of PuppiProducer #31164

Merged
merged 6 commits into from Aug 25, 2020
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 0 additions & 1 deletion CommonTools/PileupAlgos/BuildFile.xml
@@ -1,7 +1,6 @@
<use name="FWCore/ParameterSet"/>
<use name="root"/>
<use name="rootcore"/>
<use name="fastjet"/>
<export>
<lib name="1"/>
</export>
64 changes: 33 additions & 31 deletions CommonTools/PileupAlgos/src/PuppiContainer.cc
Expand Up @@ -7,8 +7,6 @@
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/isFinite.h"

#include "fastjet/PseudoJet.hh"

using namespace std;

PuppiContainer::PuppiContainer(const edm::ParameterSet &iConfig) {
Expand Down Expand Up @@ -46,28 +44,33 @@ void PuppiContainer::initialize(const std::vector<RecoObj> &iRecoObjects) {
fPFParticlesForVar.reserve(iRecoObjects.size());
fPFParticlesForVarChargedPV.reserve(iRecoObjects.size());
for (auto const &rParticle : *fRecoParticles) {
fastjet::PseudoJet curPseudoJet;
PuppiCandidate pCand;
pCand.id = rParticle.id;
if (edm::isFinite(rParticle.rapidity)) {
curPseudoJet.reset_PtYPhiM(rParticle.pt, rParticle.rapidity, rParticle.phi, rParticle.m); //hacked
pCand.pt = rParticle.pt;
pCand.eta = rParticle.eta;
pCand.rapidity = rParticle.rapidity;
pCand.phi = rParticle.phi;
pCand.m = rParticle.m;
} else {
curPseudoJet.reset_PtYPhiM(0, 99., 0, 0); //skipping may have been a better choice
pCand.pt = 0.;
pCand.eta = 99.;
pCand.rapidity = 99.;
pCand.phi = 0.;
pCand.m = 0.;
}
PuppiCandidate pCand;
pCand.pt = curPseudoJet.pt();
pCand.eta = curPseudoJet.eta();
pCand.rapidity = curPseudoJet.rap();
pCand.phi = curPseudoJet.phi();
pCand.m = curPseudoJet.m();
pCand.id = rParticle.id;

fPFParticles.push_back(pCand);
//Take Charged particles associated to PV

// skip candidates to be ignored in the computation
// of PUPPI's alphas (e.g. electrons and muons if puppiNoLep=True)
if (std::abs(rParticle.id) == 3)
continue;

fPFParticlesForVar.push_back(pCand);
// charged candidates assigned to LV
if (std::abs(rParticle.id) == 1)
fPFParticlesForVarChargedPV.push_back(pCand);
//if(rParticle.id == 3) _chargedNoPV.push_back(pCand);
// if(fNPV < rParticle.vtxId) fNPV = rParticle.vtxId;
}
}
Expand All @@ -88,17 +91,17 @@ double PuppiContainer::var_within_R(int iId,
if (iId == -1)
return 1.;

double const r2(R * R);
double var(0.);
double const r2 = R * R;
double var = 0.;

for (auto const &cand : particles) {
if (std::abs(cand.rapidity - centre.rapidity) < R) {
auto const dr2y(reco::deltaR2(cand.rapidity, cand.phi, centre.rapidity, centre.phi));
auto const dr2y = reco::deltaR2(cand.rapidity, cand.phi, centre.rapidity, centre.phi);
if (dr2y < r2) {
auto const dr2(reco::deltaR2(cand.eta, cand.phi, centre.eta, centre.phi));
auto const dr2 = reco::deltaR2(cand.eta, cand.phi, centre.eta, centre.phi);
if (dr2 < 0.0001)
continue;
auto const pt(cand.pt);
auto const pt = cand.pt;
if (iId == 5)
var += (pt * pt / dr2);
else if (iId == 4)
Expand Down Expand Up @@ -139,17 +142,17 @@ void PuppiContainer::getRMSAvg(int iOpt,
int pAlgo = fPuppiAlgo[pPupId].algoId(iOpt);
bool pCharged = fPuppiAlgo[pPupId].isCharged(iOpt);
double pCone = fPuppiAlgo[pPupId].coneSize(iOpt);
//Compute the Puppi Metric
// compute the Puppi metric:
// - calculate goodVar only for candidates that (1) will not be assigned a predefined weight (e.g 0, 1),
// or (2) are required for computations inside puppi-algos (see call to PuppiAlgo::add below)
double pVal = -1;
// calculate goodVar only for candidates that (1) will not be assigned a predefined weight (e.g 0, 1),
// or (2) are required for computations inside puppi-algos (see call to PuppiAlgo::add below)
if (((not(fApplyCHS and ((iConstits[i0].id == 1) or (iConstits[i0].id == 2)))) and (iConstits[i0].id != 3)) or
((std::abs(iConstits[i0].eta) < fPuppiAlgo[pPupId].etaMaxExtrap()) and
((iConstits[i0].id == 1) or (iConstits[i0].id == 2)))) {
bool const getsDefaultWgtIfApplyCHS = iConstits[i0].id == 1 or iConstits[i0].id == 2;
if (not((fApplyCHS and getsDefaultWgtIfApplyCHS) or iConstits[i0].id == 3) or
(std::abs(iConstits[i0].eta) < fPuppiAlgo[pPupId].etaMaxExtrap() and getsDefaultWgtIfApplyCHS)) {
pVal = goodVar(iConstits[i0], pCharged ? iChargedParticles : iParticles, pAlgo, pCone);
}
fVals.push_back(pVal);
//if(std::isnan(pVal) || std::isinf(pVal)) cerr << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt << " -- " << iConstits[i0].eta << endl;

if (!edm::isFinite(pVal)) {
LogDebug("NotFound") << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt << " -- " << iConstits[i0].eta
<< endl;
Expand All @@ -159,11 +162,10 @@ void PuppiContainer::getRMSAvg(int iOpt,
// code added by Nhan: now instead for every algorithm give it all the particles
for (int i1 = 0; i1 < fNAlgos; i1++) {
// skip cands outside of algo's etaMaxExtrap, as they would anyway be ignored inside PuppiAlgo::add (see end of the block)
if (not((std::abs(iConstits[i0].eta) < fPuppiAlgo[i1].etaMaxExtrap()) and
((iConstits[i0].id == 1) or (iConstits[i0].id == 2))))
if (not(std::abs(iConstits[i0].eta) < fPuppiAlgo[i1].etaMaxExtrap() and getsDefaultWgtIfApplyCHS))
continue;

auto curVal(pVal);
auto curVal = pVal;
// recompute goodVar if algo has changed
if (i1 != pPupId) {
pAlgo = fPuppiAlgo[i1].algoId(iOpt);
Expand Down Expand Up @@ -305,8 +307,8 @@ std::vector<double> const &PuppiContainer::puppiWeights() {
if (pWeight * fPFParticles[i0].pt < fPuppiAlgo[pPupId].neutralPt(fNPV) && rParticle.id == 0)
pWeight = 0; //threshold cut on the neutral Pt
// Protect high pT photons (important for gamma to hadronic recoil balance)
if ((fPtMaxPhotons > 0) && (rParticle.pdgId == 22) && (std::abs(fPFParticles[i0].eta) < fEtaMaxPhotons) &&
(fPFParticles[i0].pt > fPtMaxPhotons))
if (fPtMaxPhotons > 0 && rParticle.pdgId == 22 && std::abs(fPFParticles[i0].eta) < fEtaMaxPhotons &&
fPFParticles[i0].pt > fPtMaxPhotons)
pWeight = 1.;
// Protect high pT neutrals
else if ((fPtMaxNeutrals > 0) && (rParticle.id == 0))
Expand Down