Skip to content

Commit

Permalink
Merge pull request #31290 from missirol/devel_106X_puppi_speedup5
Browse files Browse the repository at this point in the history 
[10_6_X] improving performance of PuppiProducer
  • Loading branch information
@cmsbuild
cmsbuild committed Sep 4, 2020
2 parents f9e19a5 + 9ed0cf6 commit 1464697
Show file tree
Hide file tree
Showing 6 changed files with 122 additions and 103 deletions.
2 changes: 2 additions & 0 deletions CommonTools/PileupAlgos/interface/PuppiAlgo.h
View file
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
View file
@@ -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
View file
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
View file
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
View file

This file was deleted.

167 changes: 91 additions & 76 deletions CommonTools/PileupAlgos/src/PuppiContainer.cc
View file
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;
if (not(fApplyCHS and (iConstits[i0].id == 1 or iConstits[i0].id == 2))
or (std::abs(iConstits[i0].eta) < fPuppiAlgo[pPupId].etaMaxExtrap() and iConstits[i0].puppi_register != 0)) {
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;
}


0 comments on commit 1464697

Please sign in to comment.