Merge pull request #23270 from kpedro88/SpeedupPuppi1020pre3

Speedup puppi
cms-sw · May 31, 2018 · 8b2af09 · 8b2af09
2 parents c247bad + d0482fd
commit 8b2af09
Show file tree

Hide file tree

Showing 7 changed files with 55 additions and 70 deletions.
diff --git a/CommonTools/PileupAlgos/interface/PuppiAlgo.h b/CommonTools/PileupAlgos/interface/PuppiAlgo.h
@@ -3,7 +3,7 @@
 
 #include "FWCore/Framework/interface/Frameworkfwd.h"
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
-#include "fastjet/PseudoJet.hh"
+#include "CommonTools/PileupAlgos/interface/PuppiCandidate.h"
 #include <vector>
 
 class PuppiAlgo{ 
@@ -13,7 +13,7 @@ class PuppiAlgo{
   //Computing Mean and RMS
   void   reset();
   void   fixAlgoEtaBin(int i_eta );
-  void   add(const fastjet::PseudoJet &iParticle,const double &iVal,const unsigned int iAlgo);
+  void   add(const PuppiCandidate &iParticle,const double &iVal,const unsigned int iAlgo);
   void   computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac);
   //Get the Weight
   double compute(std::vector<double> const &iVals,double iChi2) const;

diff --git a/CommonTools/PileupAlgos/interface/PuppiCandidate.h b/CommonTools/PileupAlgos/interface/PuppiCandidate.h
@@ -0,0 +1,22 @@
+#ifndef CommonTools_PileupAlgos_PuppiCandidate
+#define CommonTools_PileupAlgos_PuppiCandidate
+
+#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;
+};
+
+#endif
diff --git a/CommonTools/PileupAlgos/interface/PuppiContainer.h b/CommonTools/PileupAlgos/interface/PuppiContainer.h
@@ -3,44 +3,17 @@
 
 #include "CommonTools/PileupAlgos/interface/PuppiAlgo.h"
 #include "CommonTools/PileupAlgos/interface/RecoObj.h"
-#include "fastjet/internal/base.hh"
-#include "fastjet/PseudoJet.hh"
+#include "CommonTools/PileupAlgos/interface/PuppiCandidate.h"
 
-//FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
-
-
-//......................
 class PuppiContainer{
 public:
-
-
-  // Helper class designed to store Puppi information inside of fastjet pseudojets.
-  // In CMSSW we use the user_index to refer to the index of the input collection, 
-  // but Puppi uses it to decide between NHs, PV CHs, and PU CHs. Instead,
-  // make that a register. 
-  class PuppiUserInfo : public fastjet::PseudoJet::UserInfoBase {
-   public : 
-     PuppiUserInfo( int puppi_register = -1) : puppi_register_(puppi_register) {}
-     ~PuppiUserInfo() override{}
-
-     void set_puppi_register(int i) { puppi_register_ = i; }
-
-     inline int puppi_register() const { return puppi_register_; }
-
-   protected : 
-     int puppi_register_;     /// Used by puppi algorithm to decide neutrals vs PV vs PU
-  };
-
-
-
-
     PuppiContainer(const edm::ParameterSet &iConfig);
     ~PuppiContainer(); 
     void initialize(const std::vector<RecoObj> &iRecoObjects);
     void setNPV(int iNPV){ fNPV = iNPV; }
 
-    std::vector<fastjet::PseudoJet> const & pfParticles() const { return fPFParticles; }    
-    std::vector<fastjet::PseudoJet> const & pvParticles() const { return fChargedPV; }        
+    std::vector<PuppiCandidate> const & pfParticles() const { return fPFParticles; }
+    std::vector<PuppiCandidate> const & pvParticles() const { return fChargedPV; }
     std::vector<double> const & puppiWeights();
     const std::vector<double> & puppiRawAlphas(){ return fRawAlphas; }
     const std::vector<double> & puppiAlphas(){ return fVals; }
@@ -49,21 +22,21 @@ class PuppiContainer{
     const std::vector<double> & puppiAlphasRMS() {return fAlphaRMS;}
 
     int puppiNAlgos(){ return fNAlgos; }
-    std::vector<fastjet::PseudoJet> const & puppiParticles() const { return fPupParticles;}
+    std::vector<PuppiCandidate> const & puppiParticles() const { return fPupParticles;}
 
 protected:
-    double  goodVar      (fastjet::PseudoJet const &iPart,std::vector<fastjet::PseudoJet> const &iParts, int iOpt,const double iRCone);
-    void    getRMSAvg    (int iOpt,std::vector<fastjet::PseudoJet> const &iConstits,std::vector<fastjet::PseudoJet> const &iParticles,std::vector<fastjet::PseudoJet> const &iChargeParticles);
-    void    getRawAlphas    (int iOpt,std::vector<fastjet::PseudoJet> const &iConstits,std::vector<fastjet::PseudoJet> const &iParticles,std::vector<fastjet::PseudoJet> const &iChargeParticles);
+    double  goodVar      (PuppiCandidate const &iPart,std::vector<PuppiCandidate> const &iParts, int iOpt,const double iRCone);
+    void    getRMSAvg    (int iOpt,std::vector<PuppiCandidate> const &iConstits,std::vector<PuppiCandidate> const &iParticles,std::vector<PuppiCandidate> const &iChargeParticles);
+    void    getRawAlphas    (int iOpt,std::vector<PuppiCandidate> const &iConstits,std::vector<PuppiCandidate> const &iParticles,std::vector<PuppiCandidate> const &iChargeParticles);
     double  getChi2FromdZ(double iDZ);
     int     getPuppiId   ( float iPt, float iEta);
-    double  var_within_R (int iId, const std::vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet& centre, const double R);  
+    double  var_within_R (int iId, const std::vector<PuppiCandidate> & particles, const PuppiCandidate& centre, const double R);
 
     bool      fPuppiDiagnostics;
     std::vector<RecoObj>   fRecoParticles;
-    std::vector<fastjet::PseudoJet> fPFParticles;
-    std::vector<fastjet::PseudoJet> fChargedPV;
-    std::vector<fastjet::PseudoJet> fPupParticles;
+    std::vector<PuppiCandidate> fPFParticles;
+    std::vector<PuppiCandidate> fChargedPV;
+    std::vector<PuppiCandidate> fPupParticles;
     std::vector<double>    fWeights;
     std::vector<double>    fVals;
     std::vector<double>    fRawAlphas;

diff --git a/CommonTools/PileupAlgos/plugins/PuppiProducer.cc b/CommonTools/PileupAlgos/plugins/PuppiProducer.cc
@@ -22,9 +22,8 @@
 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
 #include "DataFormats/Common/interface/Association.h"
 //Main File
-#include "fastjet/PseudoJet.hh"
 #include "CommonTools/PileupAlgos/plugins/PuppiProducer.h"
-
+#include "CommonTools/PileupAlgos/interface/PuppiCandidate.h"
 
 
 // ------------------------------------------------------------------------------------------
@@ -184,7 +183,7 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
   fPuppiContainer->setNPV( npv );
 
   std::vector<double> lWeights;
-  std::vector<fastjet::PseudoJet> lCandidates;
+  std::vector<PuppiCandidate> lCandidates;
   if (!fUseExistingWeights){
     //Compute the weights and get the particles
     lWeights = fPuppiContainer->puppiWeights();
@@ -205,7 +204,7 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
         else{ curpupweight = lPack->puppiWeight();  }
       }
       lWeights.push_back(curpupweight);
-      fastjet::PseudoJet curjet( curpupweight*lPack->px(), curpupweight*lPack->py(), curpupweight*lPack->pz(), curpupweight*lPack->energy());
+      PuppiCandidate curjet( curpupweight*lPack->px(), curpupweight*lPack->py(), curpupweight*lPack->pz(), curpupweight*lPack->energy());
       curjet.set_user_index(lPackCtr);
       lCandidates.push_back(curjet);
       lPackCtr++;
@@ -252,7 +251,7 @@ void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
     int val = i0 - i0begin;
 
     // Find the Puppi particle matched to the input collection using the "user_index" of the object. 
-    auto puppiMatched = find_if( lCandidates.begin(), lCandidates.end(), [&val]( fastjet::PseudoJet const & i ){ return i.user_index() == val; } );
+    auto puppiMatched = find_if( lCandidates.begin(), lCandidates.end(), [&val]( PuppiCandidate const & i ){ return i.user_index() == val; } );
     if ( puppiMatched != lCandidates.end() ) {
       pVec.SetPxPyPzE(puppiMatched->px(),puppiMatched->py(),puppiMatched->pz(),puppiMatched->E());
     } else {

diff --git a/CommonTools/PileupAlgos/src/PuppiAlgo.cc b/CommonTools/PileupAlgos/src/PuppiAlgo.cc
@@ -1,7 +1,6 @@
 #include "CommonTools/PileupAlgos/interface/PuppiAlgo.h"
 #include "CommonTools/PileupAlgos/interface/PuppiContainer.h"
 #include "FWCore/Utilities/interface/Exception.h"
-#include "fastjet/internal/base.hh"
 #include "Math/QuantFuncMathCore.h"
 #include "Math/SpecFuncMathCore.h"
 #include "Math/ProbFunc.h"
@@ -87,19 +86,13 @@ void PuppiAlgo::fixAlgoEtaBin(int i_eta) {
   cur_Med = fMedian_perEta[0][i_eta]; // 0 is number of algos within this eta bin
 }
 
-void PuppiAlgo::add(const fastjet::PseudoJet &iParticle,const double &iVal,const unsigned int iAlgo) {
+void PuppiAlgo::add(const PuppiCandidate &iParticle,const double &iVal,const unsigned int iAlgo) {
     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 = std::numeric_limits<int>::lowest();
-    if ( iParticle.has_user_info() ) {
-        PuppiContainer::PuppiUserInfo const * pInfo = dynamic_cast<PuppiContainer::PuppiUserInfo const *>( iParticle.user_info_ptr() );
-        if ( pInfo != nullptr ) {
-            puppi_register = pInfo->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";
     }

diff --git a/CommonTools/PileupAlgos/src/PuppiCandidate.cc b/CommonTools/PileupAlgos/src/PuppiCandidate.cc
@@ -0,0 +1 @@
+#include "CommonTools/PileupAlgos/interface/PuppiCandidate.h"
diff --git a/CommonTools/PileupAlgos/src/PuppiContainer.cc b/CommonTools/PileupAlgos/src/PuppiContainer.cc
@@ -1,7 +1,5 @@
 #include "CommonTools/PileupAlgos/interface/PuppiContainer.h"
 #include "DataFormats/Math/interface/deltaR.h"
-#include "fastjet/internal/base.hh"
-#include "fastjet/FunctionOfPseudoJet.hh"
 #include "Math/ProbFunc.h"
 #include "TMath.h"
 #include <iostream>
@@ -10,7 +8,6 @@
 #include "FWCore/Utilities/interface/isFinite.h"
 
 using namespace std;
-using namespace fastjet;
 
 PuppiContainer::PuppiContainer(const edm::ParameterSet &iConfig) {
     fPuppiDiagnostics = iConfig.getParameter<bool>("puppiDiagnostics");
@@ -44,7 +41,7 @@ void PuppiContainer::initialize(const std::vector<RecoObj> &iRecoObjects) {
     fNPV    = 1.;
     fRecoParticles = iRecoObjects;
     for (unsigned int i = 0; i < fRecoParticles.size(); i++){
-        fastjet::PseudoJet curPseudoJet;
+        PuppiCandidate curPseudoJet;
         auto fRecoParticle = fRecoParticles[i];
         // float nom = sqrt((fRecoParticle.m)*(fRecoParticle.m) + (fRecoParticle.pt)*(fRecoParticle.pt)*(cosh(fRecoParticle.eta))*(cosh(fRecoParticle.eta))) + (fRecoParticle.pt)*sinh(fRecoParticle.eta);//hacked
         // float denom = sqrt((fRecoParticle.m)*(fRecoParticle.m) + (fRecoParticle.pt)*(fRecoParticle.pt));//hacked
@@ -59,7 +56,7 @@ void PuppiContainer::initialize(const std::vector<RecoObj> &iRecoObjects) {
         if(fRecoParticle.id == 0 or fRecoParticle.charge == 0)  puppi_register = 0; // zero is neutral hadron
         if(fRecoParticle.id == 1 and fRecoParticle.charge != 0) puppi_register = fRecoParticle.charge; // from PV use the
         if(fRecoParticle.id == 2 and fRecoParticle.charge != 0) puppi_register = fRecoParticle.charge+5; // from NPV use the charge as key +5 as key
-        curPseudoJet.set_user_info( new PuppiUserInfo( puppi_register ) );
+        curPseudoJet.set_info( puppi_register );
         // fill vector of pseudojets for internal references
         fPFParticles.push_back(curPseudoJet);
         //Take Charged particles associated to PV
@@ -75,12 +72,11 @@ void PuppiContainer::initialize(const std::vector<RecoObj> &iRecoObjects) {
 }
 PuppiContainer::~PuppiContainer(){}
 
-double PuppiContainer::goodVar(PseudoJet const &iPart,std::vector<PseudoJet> const &iParts, int iOpt,const double iRCone) {
-    double lPup = 0;
-    lPup = var_within_R(iOpt,iParts,iPart,iRCone);
-    return lPup;
+double PuppiContainer::goodVar(PuppiCandidate const &iPart,std::vector<PuppiCandidate> const &iParts, int iOpt,const double iRCone) {
+    return var_within_R(iOpt,iParts,iPart,iRCone);
 }
-double PuppiContainer::var_within_R(int iId, const vector<PseudoJet> & particles, const PseudoJet& centre, const double R){
+
+double PuppiContainer::var_within_R(int iId, const vector<PuppiCandidate> & particles, const PuppiCandidate& centre, const double R){
     if(iId == -1) return 1;
 
     //this is a circle in rapidity-phi
@@ -92,10 +88,11 @@ double PuppiContainer::var_within_R(int iId, const vector<PseudoJet> & particles
 
     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){
-      if ( part.squared_distance(centre) < R*R ){
-	near_dR2s.push_back(reco::deltaR2(part, centre));
-	near_pts.push_back(part.pt());
+      if ( part.squared_distance(centre) < r2 ){
+        near_dR2s.push_back(reco::deltaR2(part, centre));
+        near_pts.push_back(part.pt());
       }
     }
     double var = 0;
@@ -120,7 +117,7 @@ double PuppiContainer::var_within_R(int iId, const vector<PseudoJet> & particles
     return var;
 }
 //In fact takes the median not the average
-void PuppiContainer::getRMSAvg(int iOpt,std::vector<fastjet::PseudoJet> const &iConstits,std::vector<fastjet::PseudoJet> const &iParticles,std::vector<fastjet::PseudoJet> const &iChargedParticles) {
+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
@@ -160,7 +157,7 @@ void PuppiContainer::getRMSAvg(int iOpt,std::vector<fastjet::PseudoJet> const &i
     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<fastjet::PseudoJet> const &iConstits,std::vector<fastjet::PseudoJet> const &iParticles,std::vector<fastjet::PseudoJet> const &iChargedParticles) {
+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;
@@ -275,7 +272,7 @@ std::vector<double> const & PuppiContainer::puppiWeights() {
         // if(std::abs(pWeight) <= 0. ) continue; 
 
         //Produce
-        PseudoJet curjet( pWeight*fPFParticles[i0].px(), pWeight*fPFParticles[i0].py(), pWeight*fPFParticles[i0].pz(), pWeight*fPFParticles[i0].e() );
+        PuppiCandidate curjet( pWeight*fPFParticles[i0].px(), pWeight*fPFParticles[i0].py(), pWeight*fPFParticles[i0].pz(), pWeight*fPFParticles[i0].e() );
         curjet.set_user_index(i0);
         fPupParticles.push_back(curjet);
     }