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

[10_6_X] improving performance of PuppiProducer #31290

Merged
merged 2 commits into from Sep 4, 2020
Merged
Show file tree
Hide file tree
Changes from 1 commit
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
2 changes: 2 additions & 0 deletions CommonTools/PileupAlgos/interface/PuppiAlgo.h
Expand Up @@ -35,6 +35,8 @@ class PuppiAlgo{
inline double rms() const {return cur_RMS;}
inline double median() const {return cur_Med;}

inline double etaMaxExtrap() const { return fEtaMaxExtrap; }

private:
unsigned int fNAlgos;
std::vector<double> fEtaMax;
Expand Down
32 changes: 14 additions & 18 deletions CommonTools/PileupAlgos/interface/PuppiCandidate.h
@@ -1,22 +1,18 @@
#ifndef CommonTools_PileupAlgos_PuppiCandidate
#define CommonTools_PileupAlgos_PuppiCandidate
#ifndef CommonTools_PileupAlgos_PuppiCandidate_h
#define CommonTools_PileupAlgos_PuppiCandidate_h

#include "fastjet/PseudoJet.hh"

// Extension of fastjet::PseudoJet that caches eta and some other info
// Puppi uses register to decide between NHs, PV CHs, and PU CHs.
class PuppiCandidate : public fastjet::PseudoJet {
public:
using fastjet::PseudoJet::PseudoJet;
double pseudorapidity() const { _ensure_valid_eta(); return _eta; }
double eta() const { return pseudorapidity(); }
void _ensure_valid_eta() const { if(_eta==fastjet::pseudojet_invalid_rap) _eta = fastjet::PseudoJet::pseudorapidity(); }
void set_info(int puppi_register) { _puppi_register = puppi_register; }
inline int puppi_register() const { return _puppi_register; }
private:
// variable names in fastjet style
mutable double _eta = fastjet::pseudojet_invalid_rap;
int _puppi_register;
struct PuppiCandidate {
double pt{0};
double eta{0};
double phi{0};
double m{0};
double rapidity{0};
double px{0};
double py{0};
double pz{0};
double e{0};
int id{0};
int puppi_register{0};
};

#endif
14 changes: 11 additions & 3 deletions CommonTools/PileupAlgos/plugins/PuppiProducer.cc
Expand Up @@ -236,8 +236,16 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
else{ curpupweight = lPack->puppiWeight(); }
}
lWeights.push_back(curpupweight);
PuppiCandidate curjet( curpupweight*lPack->px(), curpupweight*lPack->py(), curpupweight*lPack->pz(), curpupweight*lPack->energy());
curjet.set_user_index(lPackCtr);
PuppiCandidate curjet;
curjet.px = curpupweight*lPack->px();
curjet.py = curpupweight*lPack->py();
curjet.pz = curpupweight*lPack->pz();
curjet.e = curpupweight*lPack->energy();
curjet.pt = curpupweight*lPack->pt();
curjet.eta = lPack->eta();
curjet.rapidity = lPack->rapidity();
curjet.phi = lPack->phi();
curjet.m = curpupweight*lPack->mass();
lCandidates.push_back(curjet);
lPackCtr++;
}
Expand Down Expand Up @@ -285,7 +293,7 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
int iPuppiMatched = fUseExistingWeights ? val : fPuppiContainer->recoToPup()[val];
if ( iPuppiMatched >= 0 ) {
auto const& puppiMatched = lCandidates[iPuppiMatched];
pVec.SetPxPyPzE(puppiMatched.px(),puppiMatched.py(),puppiMatched.pz(),puppiMatched.E());
pVec.SetPxPyPzE(puppiMatched.px,puppiMatched.py,puppiMatched.pz,puppiMatched.e);
if(fClonePackedCands && (!fUseExistingWeights)) {
if(fPuppiForLeptons)
pCand->setPuppiWeight(pCand->puppiWeight(),lWeights[val]);
Expand Down
9 changes: 4 additions & 5 deletions CommonTools/PileupAlgos/src/PuppiAlgo.cc
Expand Up @@ -87,12 +87,12 @@ void PuppiAlgo::fixAlgoEtaBin(int i_eta) {
}

void PuppiAlgo::add(const PuppiCandidate &iParticle,const double &iVal,const unsigned int iAlgo) {
if(iParticle.pt() < fRMSPtMin[iAlgo]) return;
if(iParticle.pt < fRMSPtMin[iAlgo]) return;
// Change from SRR : Previously used fastjet::PseudoJet::user_index to decide the particle type.
// In CMSSW we use the user_index to specify the index in the input collection, so I invented
// a new mechanism using the fastjet UserInfo functionality. Of course, it's still just an integer
// but that interface could be changed (or augmented) if desired / needed.
int puppi_register = iParticle.puppi_register();
int puppi_register = iParticle.puppi_register;
if ( puppi_register == std::numeric_limits<int>::lowest() ) {
throw cms::Exception("PuppiRegisterNotSet") << "The puppi register is not set. This must be set before use.\n";
}
Expand All @@ -107,13 +107,13 @@ void PuppiAlgo::add(const PuppiCandidate &iParticle,const double &iVal,const uns

// added by Nhan -- for all eta regions, compute mean/RMS from the central charged PU
//std::cout << "std::abs(puppi_register) = " << std::abs(puppi_register) << std::endl;
if ((std::abs(iParticle.eta()) < fEtaMaxExtrap) && (std::abs(puppi_register) >= 3)){
if ((std::abs(iParticle.eta) < fEtaMaxExtrap) && (std::abs(puppi_register) >= 3)){
fPups.push_back(iVal);
// fPupsPV.push_back(iVal);
fNCount[iAlgo]++;
}
// for the low PU case, correction. for checking that the PU-only median will be below the PV particles
if(std::abs(iParticle.eta()) < fEtaMaxExtrap && (std::abs(puppi_register) >=1 && std::abs(puppi_register) <=2)) fPupsPV.push_back(iVal);
if(std::abs(iParticle.eta) < fEtaMaxExtrap && (std::abs(puppi_register) >=1 && std::abs(puppi_register) <=2)) fPupsPV.push_back(iVal);

}

Expand Down Expand Up @@ -162,7 +162,6 @@ void PuppiAlgo::computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac) {

if(fAdjust[iAlgo]){
//Adjust the p-value to correspond to the median
std::sort(fPupsPV.begin(),fPupsPV.end());
int lNPV = 0;
for(unsigned int i0 = 0; i0 < fPupsPV.size(); i0++) if(fPupsPV[i0] <= lMed ) lNPV++;
double lAdjust = double(lNPV)/double(lNPV+0.5*fNCount[iAlgo]);
Expand Down
1 change: 0 additions & 1 deletion CommonTools/PileupAlgos/src/PuppiCandidate.cc

This file was deleted.

167 changes: 91 additions & 76 deletions CommonTools/PileupAlgos/src/PuppiContainer.cc
Expand Up @@ -7,6 +7,8 @@
#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 @@ -35,36 +37,44 @@ void PuppiContainer::initialize(const std::vector<RecoObj> &iRecoObjects) {
fRawAlphas.resize(0);
fAlphaMed .resize(0);
fAlphaRMS .resize(0);
//fChargedNoPV.resize(0);
//Link to the RecoObjects
fPVFrac = 0.;
fNPV = 1.;
//Link to the RecoObjects
fRecoParticles = &iRecoObjects;
fPFParticles.reserve(iRecoObjects.size());
fChargedPV.reserve(iRecoObjects.size());
fRecoToPup.clear();
fRecoToPup.reserve(fRecoParticles->size());
for (auto const& rParticle : *fRecoParticles){
PuppiCandidate curPseudoJet;
// float nom = sqrt((rParticle.m)*(rParticle.m) + (rParticle.pt)*(rParticle.pt)*(cosh(rParticle.eta))*(cosh(rParticle.eta))) + (rParticle.pt)*sinh(rParticle.eta);//hacked
// float denom = sqrt((rParticle.m)*(rParticle.m) + (rParticle.pt)*(rParticle.pt));//hacked
// float rapidity = log(nom/denom);//hacked
// usage of fastjet::PseudoJet is only needed to enforce
// the no-change policy for backports (no numerical differences in outputs)
fastjet::PseudoJet curPseudoJet;
if (edm::isFinite(rParticle.rapidity)){
curPseudoJet.reset_PtYPhiM(rParticle.pt,rParticle.rapidity,rParticle.phi,rParticle.m);//hacked
} else {
curPseudoJet.reset_PtYPhiM(0, 99., 0, 0);//skipping may have been a better choice
}
//curPseudoJet.reset_PtYPhiM(rParticle.pt,rParticle.eta,rParticle.phi,rParticle.m);
int puppi_register = 0;
if(rParticle.id == 0 or rParticle.charge == 0) puppi_register = 0; // zero is neutral hadron
if(rParticle.id == 1 and rParticle.charge != 0) puppi_register = rParticle.charge; // from PV use the
if(rParticle.id == 2 and rParticle.charge != 0) puppi_register = rParticle.charge+5; // from NPV use the charge as key +5 as key
curPseudoJet.set_info( puppi_register );
PuppiCandidate pCand;
pCand.id = rParticle.id;
pCand.puppi_register = puppi_register;
pCand.pt = curPseudoJet.pt();
pCand.rapidity = curPseudoJet.rap();
pCand.eta = curPseudoJet.eta();
pCand.phi = curPseudoJet.phi();
pCand.m = curPseudoJet.m();
pCand.px = curPseudoJet.px();
pCand.py = curPseudoJet.py();
pCand.pz = curPseudoJet.pz();
pCand.e = curPseudoJet.E();
// fill vector of pseudojets for internal references
fPFParticles.push_back(curPseudoJet);
fPFParticles.push_back(pCand);
//Take Charged particles associated to PV
if(std::abs(rParticle.id) == 1) fChargedPV.push_back(curPseudoJet);
if(std::abs(rParticle.id) == 1) fChargedPV.push_back(pCand);
if(std::abs(rParticle.id) >= 1 ) fPVFrac+=1.;
//if(rParticle.id == 3) _chargedNoPV.push_back(curPseudoJet);
// if(fNPV < rParticle.vtxId) fNPV = rParticle.vtxId;
}
if (fPVFrac != 0) fPVFrac = double(fChargedPV.size())/fPVFrac;
else fPVFrac = 0;
Expand All @@ -76,52 +86,49 @@ double PuppiContainer::goodVar(PuppiCandidate const &iPart,std::vector<PuppiCand
}

double PuppiContainer::var_within_R(int iId, const vector<PuppiCandidate> & particles, const PuppiCandidate& centre, const double R){
if(iId == -1) return 1;
if (iId == -1)
return 1.;

//this is a circle in rapidity-phi
//it would make more sense to have var definition consistent
//fastjet::Selector sel = fastjet::SelectorCircle(R);
//sel.set_reference(centre);
//the original code used Selector infrastructure: it is too heavy here
//logic of SelectorCircle is preserved below
double const r2 = R * R;
double var = 0.;

vector<double > near_dR2s; near_dR2s.reserve(std::min(50UL, particles.size()));
vector<double > near_pts; near_pts.reserve(std::min(50UL, particles.size()));
const double r2 = R*R;
for (auto const& part : particles){
//squared_distance is in (y,phi) coords: rap() has faster access -> check it first
if ( std::abs(part.rap()-centre.rap()) < R && part.squared_distance(centre) < r2 ){
near_dR2s.push_back(reco::deltaR2(part, centre));
near_pts.push_back(part.pt());
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);
if (dr2y < r2) {
auto const dr2 = reco::deltaR2(cand.eta, cand.phi, centre.eta, centre.phi);
if (dr2 < 0.0001)
continue;
auto const pt = cand.pt;
if (iId == 5)
var += (pt * pt / dr2);
else if (iId == 4)
var += pt;
else if (iId == 3)
var += (1. / dr2);
else if (iId == 2)
var += (1. / dr2);
else if (iId == 1)
var += pt;
else if (iId == 0)
var += (pt / dr2);
}
}
}
double var = 0;
//double lSumPt = 0;
//if(iId == 1) for(auto pt : near_pts) lSumPt += pt;
auto nParts = near_dR2s.size();
for(auto i = 0UL; i < nParts; ++i){
auto dr2 = near_dR2s[i];
auto pt = near_pts[i];
if(dr2 < 0.0001) continue;
if(iId == 0) var += (pt/dr2);
else if(iId == 1) var += pt;
else if(iId == 2) var += (1./dr2);
else if(iId == 3) var += (1./dr2);
else if(iId == 4) var += pt;
else if(iId == 5) var += (pt * pt/dr2);
}
if(iId == 1) var += centre.pt(); //Sum in a cone
else if(iId == 0 && var != 0) var = log(var);
else if(iId == 3 && var != 0) var = log(var);
else if(iId == 5 && var != 0) var = log(var);

if ((var != 0.) and ((iId == 0) or (iId == 3) or (iId == 5)))
var = log(var);
else if (iId == 1)
var += centre.pt;

return var;
}

//In fact takes the median not the average
void PuppiContainer::getRMSAvg(int iOpt,std::vector<PuppiCandidate> const &iConstits,std::vector<PuppiCandidate> const &iParticles,std::vector<PuppiCandidate> const &iChargedParticles) {
for(unsigned int i0 = 0; i0 < iConstits.size(); i0++ ) {
double pVal = -1;
//Calculate the Puppi Algo to use
int pPupId = getPuppiId(iConstits[i0].pt(),iConstits[i0].eta());
int pPupId = getPuppiId(iConstits[i0].pt,iConstits[i0].eta);
if(pPupId == -1 || fPuppiAlgo[pPupId].numAlgos() <= iOpt){
fVals.push_back(-1);
continue;
Expand All @@ -130,51 +137,55 @@ void PuppiContainer::getRMSAvg(int iOpt,std::vector<PuppiCandidate> const &iCons
int pAlgo = fPuppiAlgo[pPupId].algoId (iOpt);
bool pCharged = fPuppiAlgo[pPupId].isCharged(iOpt);
double pCone = fPuppiAlgo[pPupId].coneSize (iOpt);
//Compute the Puppi Metric
if(!pCharged) pVal = goodVar(iConstits[i0],iParticles ,pAlgo,pCone);
if( pCharged) pVal = goodVar(iConstits[i0],iChargedParticles,pAlgo,pCone);
// 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;
bool const getsDefaultWgtIfApplyCHS = iConstits[i0].id == 1 or iConstits[i0].id == 2;
if (not(fApplyCHS and getsDefaultWgtIfApplyCHS) or (std::abs(iConstits[i0].eta) < fPuppiAlgo[pPupId].etaMaxExtrap() and iConstits[i0].puppi_register != 0)) {
Copy link
Contributor

Choose a reason for hiding this comment

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

The master version uses and getsDefaultWgtIfApplyCHS in the second conditional instead of and iConstits[i0].puppi_register != 0. Are these equivalent? If yes, a redefinition of getsDefaultWgtIfApplyCHS may be in order.
If no, please elaborate
(in this case, the separate definition of getsDefaultWgtIfApplyCHS seems unnecessary because it is not reused anywhere else).

also a linebreak at/after or would be nice

Copy link
Contributor Author

@missirol missirol Sep 3, 2020

Choose a reason for hiding this comment

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

In practice, I think the condition id == 1 or id == 2 corresponds to register != 0 (looking at how puppi_register is filled here), but in principle the two things can differ (or at least, I wanted to make no assumptions).

The logic behind how these changes are implemented is the following: goodVar needs to be calculated when the final weight of the candidate is non-trivial (this is the first condition, which corresponds to the logical NOT of these two cases), or when the goodVar of a given candidate is needed internally by PuppiAlgo::add, here; in 10_6_X, the latter condition translates to (.. < etaMaxExtrap) and puppi_register != 0, while in 11_X_Y the criterion is (.. < etaMaxExtrap) and (id == 1 or id == 2) (see here).

So, the logic is the same here and in the original PR; the end result looks different just as a consequence of the fact that the Puppi{Container,Algo} implementation is different in the two releases.

I added a commit to remove the temporary variable getsDefaultWgtIfApplyCHS, which was indeed only used once, and took the opportunity to add the linebreak.

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;
LogDebug( "NotFound" ) << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt << " -- " << iConstits[i0].eta << endl;
continue;
}

// // fPuppiAlgo[pPupId].add(iConstits[i0],pVal,iOpt);
//code added by Nhan, now instead for every algorithm give it all the particles
// code added by Nhan: now instead for every algorithm give it all the particles
for(int i1 = 0; i1 < fNAlgos; i1++){
pAlgo = fPuppiAlgo[i1].algoId (iOpt);
pCharged = fPuppiAlgo[i1].isCharged(iOpt);
pCone = fPuppiAlgo[i1].coneSize (iOpt);
double curVal = -1;
if (i1 != pPupId){
if(!pCharged) curVal = goodVar(iConstits[i0],iParticles ,pAlgo,pCone);
if( pCharged) curVal = goodVar(iConstits[i0],iChargedParticles,pAlgo,pCone);
} else {//no need to repeat the computation
curVal = pVal;
// 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].puppi_register != 0))
continue;

auto curVal = pVal;
// recompute goodVar if algo has changed
if (i1 != pPupId) {
pAlgo = fPuppiAlgo[i1].algoId (iOpt);
pCharged = fPuppiAlgo[i1].isCharged(iOpt);
pCone = fPuppiAlgo[i1].coneSize (iOpt);
curVal = goodVar(iConstits[i0],pCharged ? iChargedParticles : iParticles,pAlgo,pCone);
}
//std::cout << "i1 = " << i1 << ", curVal = " << curVal << ", eta = " << iConstits[i0].eta() << ", pupID = " << pPupId << std::endl;

fPuppiAlgo[i1].add(iConstits[i0],curVal,iOpt);
}

}
for(int i0 = 0; i0 < fNAlgos; i0++) fPuppiAlgo[i0].computeMedRMS(iOpt,fPVFrac);
}

//In fact takes the median not the average
void PuppiContainer::getRawAlphas(int iOpt,std::vector<PuppiCandidate> const &iConstits,std::vector<PuppiCandidate> const &iParticles,std::vector<PuppiCandidate> const &iChargedParticles) {
for(int j0 = 0; j0 < fNAlgos; j0++){
for(unsigned int i0 = 0; i0 < iConstits.size(); i0++ ) {
double pVal = -1;
//Get the Puppi Sub Algo (given iteration)
int pAlgo = fPuppiAlgo[j0].algoId (iOpt);
bool pCharged = fPuppiAlgo[j0].isCharged(iOpt);
double pCone = fPuppiAlgo[j0].coneSize (iOpt);
//Compute the Puppi Metric
if(!pCharged) pVal = goodVar(iConstits[i0],iParticles ,pAlgo,pCone);
if( pCharged) pVal = goodVar(iConstits[i0],iChargedParticles,pAlgo,pCone);
double const pVal = goodVar(iConstits[i0],pCharged ? iChargedParticles : iParticles,pAlgo,pCone);
fRawAlphas.push_back(pVal);
if( ! edm::isFinite(pVal)) {
LogDebug( "NotFound" ) << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt() << " -- " << iConstits[i0].eta() << endl;
LogDebug( "NotFound" ) << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt << " -- " << iConstits[i0].eta << endl;
continue;
}
}
Expand Down Expand Up @@ -269,9 +280,9 @@ std::vector<double> const & PuppiContainer::puppiWeights() {
LogDebug("PuppiWeightError") << "====> Weight is nan : " << pWeight << " : pt " << rParticle.pt << " -- eta : " << rParticle.eta << " -- Value" << fVals[i0] << " -- id : " << rParticle.id << " -- NAlgos: " << lNAlgos << std::endl;
}
//Basic Cuts
if(pWeight*fPFParticles[i0].pt() < fPuppiAlgo[pPupId].neutralPt(fNPV) && rParticle.id == 0 ) pWeight = 0; //threshold cut on the neutral Pt
if(pWeight*fPFParticles[i0].pt < fPuppiAlgo[pPupId].neutralPt(fNPV) && rParticle.id == 0 ) pWeight = 0; //threshold cut on the neutral Pt
if((fPtMax>0) && (rParticle.id == 0))
pWeight = std::clamp((fPFParticles[i0].pt() - fPtMaxNeutralsStartSlope) / (fPtMax - fPtMaxNeutralsStartSlope),
pWeight = std::clamp((fPFParticles[i0].pt - fPtMaxNeutralsStartSlope) / (fPtMax - fPtMaxNeutralsStartSlope),
pWeight,
1.);
if(pWeight < fPuppiWeightCut) pWeight = 0; //==> Elminate the low Weight stuff
Expand All @@ -288,11 +299,15 @@ std::vector<double> const & PuppiContainer::puppiWeights() {
// if(std::abs(pWeight) <= 0. ) continue;

//Produce
PuppiCandidate curjet( pWeight*fPFParticles[i0].px(), pWeight*fPFParticles[i0].py(), pWeight*fPFParticles[i0].pz(), pWeight*fPFParticles[i0].e() );
curjet.set_user_index(i0);
PuppiCandidate curjet(fPFParticles[i0]);
curjet.pt *= pWeight;
curjet.m *= pWeight;
curjet.px *= pWeight;
curjet.py *= pWeight;
curjet.pz *= pWeight;
curjet.e *= pWeight;

fPupParticles.push_back(curjet);
}
return fWeights;
}