Changeset 35cdc46 in git for external/fastjet/contribs/Nsubjettiness/AxesFinder.hh

Timestamp:

Sep 3, 2014, 3:18:54 PM (10 years ago)

Author:

Pavel Demin <demin@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

be2222c

Parents:

5b5a56b

Message:

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

File:

• external/fastjet/contribs/Nsubjettiness/AxesFinder.hh (modified) (15 diffs)

external/fastjet/contribs/Nsubjettiness/AxesFinder.hh

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: AxesFinder.hh 678 2014-06-12 20:43:03Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -53 +54 @@
 class AxesFinder {
    
-protected:
-   AxesFinder* _startingFinder; // storing a possible starting finder if needed
-   std::vector<fastjet::PseudoJet> _seedAxes;
-   
-   AxesFinder(AxesFinder* startingFinder = NULL) : _startingFinder(startingFinder) {}
-   
-public:
-   virtual ~AxesFinder(){
-      if (_startingFinder) delete _startingFinder;  //TODO: Convert to smart pointers to avoid this.
-   }
-   
-   // Allow setting of seedAxes from a starting finder
-   std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector<fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
-      if (_startingFinder) {
-         _seedAxes = _startingFinder->getAxes(n_jets,inputs,currentAxes);
-         return getBetterAxes(n_jets,inputs,_seedAxes);
-      } else {
-         _seedAxes = getBetterAxes(n_jets,inputs,currentAxes);
-         return _seedAxes;
-      }
-   }
-   
-   // say what the current seed axes are
-   std::vector<fastjet::PseudoJet> seedAxes() const {
-      return _seedAxes;
-   }
-   
-   // This function should be overloaded, and updates the seedAxes
-   virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector<fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& seedAxes) = 0;
-   
-};
-
+public:
+   
+   // This function should be overloaded, and updates the seedAxes to return new axes
+   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets,
+                                                   const std::vector<fastjet::PseudoJet>& inputs,
+                                                   const std::vector<fastjet::PseudoJet>& seedAxes) const = 0;
+   // convenient shorthand for squaring
+   static inline double sq(double x) {return x*x;}
+
+   //virtual destructor
+   virtual ~AxesFinder(){}
+   
+};
+
+   
 //------------------------------------------------------------------------
 /// \class AxesFinderFromExclusiveJetDefinition
@@ -90 +74 @@
 // with different jet algorithms.
 class AxesFinderFromExclusiveJetDefinition : public AxesFinder {
-
-   private:
-      fastjet::JetDefinition _def;
-   
-   public:
-      AxesFinderFromExclusiveJetDefinition(fastjet::JetDefinition def) : _def(def) {}
-      
-      virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
-         fastjet::ClusterSequence jet_clust_seq(inputs, _def);
-         return jet_clust_seq.exclusive_jets(n_jets);
-      }
+   
+public:
+   AxesFinderFromExclusiveJetDefinition(fastjet::JetDefinition def)
+   : _def(def) {}
+   
+   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets,
+                                                   const std::vector <fastjet::PseudoJet> & inputs,
+                                                   const std::vector<fastjet::PseudoJet>& /*seedAxes*/) const {
+      fastjet::ClusterSequence jet_clust_seq(inputs, _def);
+      return jet_clust_seq.exclusive_jets(n_jets);
+   }
+   
+private:
+   fastjet::JetDefinition _def;
+
 };
 
@@ -108 +96 @@
 // winner take all recombination scheme.
 class AxesFinderFromWTA_KT : public AxesFinderFromExclusiveJetDefinition { 
-   private: 
-      const WinnerTakeAllRecombiner *recomb;
-   public:
-      AxesFinderFromWTA_KT() : AxesFinderFromExclusiveJetDefinition(
-         fastjet::JetDefinition(fastjet::kt_algorithm, 
-         fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
-         recomb = new WinnerTakeAllRecombiner(), 
-         fastjet::Best)) {}
-      ~AxesFinderFromWTA_KT() {delete recomb;}
-   };
+
+public:
+   AxesFinderFromWTA_KT()
+   : AxesFinderFromExclusiveJetDefinition(
+      fastjet::JetDefinition(fastjet::kt_algorithm,
+      fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
+      &_recomb,
+      fastjet::Best)) {}
+   
+private:
+   const WinnerTakeAllRecombiner _recomb;
+
+};
    
 //------------------------------------------------------------------------
@@ -124 +115 @@
 // winner take all recombination scheme.
 class AxesFinderFromWTA_CA : public AxesFinderFromExclusiveJetDefinition {
-   private: 
-      const WinnerTakeAllRecombiner *recomb;
-   public:
-      AxesFinderFromWTA_CA() : AxesFinderFromExclusiveJetDefinition(
-         fastjet::JetDefinition(fastjet::cambridge_algorithm, 
-         fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
-         recomb = new WinnerTakeAllRecombiner(), 
-         fastjet::Best)) {}
-      ~AxesFinderFromWTA_CA() {delete recomb;}
-};
-
-//  The following classes are for testing, and are commented out for initial release
-//
-////------------------------------------------------------------------------
-///// \class AxesFinderFromWTA2_KT
-//// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a 
-//// winner take all recombination scheme with alpha = 2.
-//class AxesFinderFromWTA2_KT : public AxesFinderFromExclusiveJetDefinition { 
-//   private: 
-//      const WinnerTakeAllRecombiner *recomb;
-//   public:
-//      AxesFinderFromWTA2_KT() : AxesFinderFromExclusiveJetDefinition(
-//         fastjet::JetDefinition(fastjet::kt_algorithm, 
-//         fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
-//         recomb = new WinnerTakeAllRecombiner(2), // uses alpha = 2 here
-//         fastjet::Best)) {}
-//      ~AxesFinderFromWTA2_KT() {delete recomb;}
-//   };
-//   
-////------------------------------------------------------------------------
-///// \class AxesFinderFromWTA2_CA
-//// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a 
-//// winner take all recombination scheme with alpha = 2.
-//class AxesFinderFromWTA2_CA : public AxesFinderFromExclusiveJetDefinition {
-//   private: 
-//      const WinnerTakeAllRecombiner *recomb;
-//   public:
-//      AxesFinderFromWTA2_CA() : AxesFinderFromExclusiveJetDefinition(
-//         fastjet::JetDefinition(fastjet::cambridge_algorithm, 
-//         fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
-//         recomb = new WinnerTakeAllRecombiner(2), //uses alpha = 2 here
-//         fastjet::Best)) {}
-//      ~AxesFinderFromWTA2_CA() {delete recomb;}
-//};
+public:
+   AxesFinderFromWTA_CA()
+   : AxesFinderFromExclusiveJetDefinition(
+      fastjet::JetDefinition(fastjet::cambridge_algorithm, 
+      fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
+      &_recomb,
+      fastjet::Best)) {}
+   
+private:
+   const WinnerTakeAllRecombiner _recomb;
+};
+
 
 //------------------------------------------------------------------------
@@ -174 +133 @@
 // E_scheme recombination.
 class AxesFinderFromKT : public AxesFinderFromExclusiveJetDefinition {
-   public:
-      AxesFinderFromKT() : AxesFinderFromExclusiveJetDefinition(
-         fastjet::JetDefinition(fastjet::kt_algorithm,
-         fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
-         fastjet::E_scheme,
-         fastjet::Best)) {}
+public:
+   AxesFinderFromKT()
+   : AxesFinderFromExclusiveJetDefinition(
+      fastjet::JetDefinition(fastjet::kt_algorithm,
+      fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
+      fastjet::E_scheme,
+      fastjet::Best)) {}
 };
 
@@ -187 +147 @@
 // E_scheme recombination.
 class AxesFinderFromCA : public AxesFinderFromExclusiveJetDefinition { 
-   public:
-      AxesFinderFromCA() : AxesFinderFromExclusiveJetDefinition(
-         fastjet::JetDefinition(fastjet::cambridge_algorithm,
-                                fastjet::JetDefinition::max_allowable_R,  //maximum jet radius constant
-         fastjet::E_scheme,
-         fastjet::Best)) {}
+public:
+   AxesFinderFromCA()
+   : AxesFinderFromExclusiveJetDefinition(
+      fastjet::JetDefinition(fastjet::cambridge_algorithm,
+                             fastjet::JetDefinition::max_allowable_R,  //maximum jet radius constant
+      fastjet::E_scheme,
+      fastjet::Best)) {}
 };
 
@@ -201 +162 @@
 // This can be implemented with different jet algorithms.
 class AxesFinderFromHardestJetDefinition : public AxesFinder {
-
-   private:
-      fastjet::JetDefinition _def;
-   
-   public:
-      AxesFinderFromHardestJetDefinition(fastjet::JetDefinition def) : _def(def) {}
-      
-      virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
-         fastjet::ClusterSequence jet_clust_seq(inputs, _def);
-         std::vector<fastjet::PseudoJet> myJets = sorted_by_pt(jet_clust_seq.inclusive_jets());
-         myJets.resize(n_jets);  // only keep n hardest
-         return myJets;
-      }      
+public:
+   AxesFinderFromHardestJetDefinition(fastjet::JetDefinition def)
+   : _def(def) {}
+   
+   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets,
+                                                   const std::vector <fastjet::PseudoJet> & inputs,
+                                                   const std::vector<fastjet::PseudoJet>& /*seedAxes*/) const {
+      fastjet::ClusterSequence jet_clust_seq(inputs, _def);
+      std::vector<fastjet::PseudoJet> myJets = sorted_by_pt(jet_clust_seq.inclusive_jets());
+      myJets.resize(n_jets);  // only keep n hardest
+      return myJets;
+   }
+   
+private:
+   fastjet::JetDefinition _def;
 };
 
@@ -221 +184 @@
 // to an anti kT algorithm and E_scheme.
 class AxesFinderFromAntiKT : public AxesFinderFromHardestJetDefinition {
-   public:
-      AxesFinderFromAntiKT(double R0) : AxesFinderFromHardestJetDefinition(fastjet::JetDefinition(fastjet::antikt_algorithm,R0,fastjet::E_scheme,fastjet::Best)) {}
+public:
+   AxesFinderFromAntiKT(double R0)
+   : AxesFinderFromHardestJetDefinition(
+      fastjet::JetDefinition(fastjet::antikt_algorithm,
+                             R0,fastjet::E_scheme,fastjet::Best)) {}
 };
 
@@ -231 +197 @@
 class AxesFinderFromUserInput : public AxesFinder {
 
-   public:
-      AxesFinderFromUserInput() {}
-      
-      virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
-         assert(currentAxes.size() == (unsigned int) n_jets);
-         return currentAxes;
-      }
+public:
+   AxesFinderFromUserInput() {}
+   
+   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & /*inputs*/, const std::vector<fastjet::PseudoJet>& currentAxes) const {
+      assert(currentAxes.size() == (unsigned int) n_jets);
+      (void)(n_jets);  // adding this line to fix unused-parameter warning
+      return currentAxes;
+   }
 };
 
@@ -249 +216 @@
 class AxesFinderFromOnePassMinimization : public AxesFinder {
 
-   private:
-      double _precision;  // Desired precision in axes alignment
-      int _halt;  // maximum number of steps per iteration
-      
-      double _beta;
-      double _Rcutoff;
-      
-      DefaultUnnormalizedMeasure _measureFunction;
-   
-   public:
-
-      // From a startingFinder, try to minimize the unnormalized_measure
-      AxesFinderFromOnePassMinimization(AxesFinder* startingFinder, double beta, double Rcutoff)
-         : AxesFinder(startingFinder),
-           _precision(0.0001), //hard coded for now
-           _halt(1000), //hard coded for now
-           _beta(beta),
-           _Rcutoff(Rcutoff),
-           _measureFunction(beta, Rcutoff)
-           {}
-   
-      virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes);
-
-      template <int N> std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
-                                  const std::vector <fastjet::PseudoJet> & inputJets);
-   
-      std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes, 
-                                      const std::vector <fastjet::PseudoJet> & inputJets);
+public:
+
+   // From a startingFinder, try to minimize the unnormalized_measure
+   AxesFinderFromOnePassMinimization(double beta, double Rcutoff)
+      : _precision(0.0001), //hard coded for now
+        _halt(1000), //hard coded for now
+        _beta(beta),
+        _Rcutoff(Rcutoff),
+        _measureFunction(beta, Rcutoff)
+        {}
+
+   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets,
+                                                   const std::vector <fastjet::PseudoJet> & inputJets,
+                                                   const std::vector<fastjet::PseudoJet>& currentAxes) const;
+   
+private:
+   double _precision;  // Desired precision in axes alignment
+   int _halt;  // maximum number of steps per iteration
+   
+   double _beta;
+   double _Rcutoff;
+   
+   DefaultUnnormalizedMeasureFunction _measureFunction;
+   
+   template <int N> std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
+                                                              const std::vector <fastjet::PseudoJet> & inputJets) const;
+   
+   std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
+                                         const std::vector <fastjet::PseudoJet> & inputJets) const;
 
 };
@@ -288 +256 @@
 class AxesFinderFromKmeansMinimization : public AxesFinder{
 
-   private:
-      int _n_iterations;   // Number of iterations to run  (0 for no minimization, 1 for one-pass, >>1 for global minimum)
-      double _noise_range; // noise range for random initialization
-   
-      DefaultUnnormalizedMeasure _measureFunction; //function to test whether minimum is reached
-   
-      AxesFinderFromOnePassMinimization _onePassFinder;  //one pass finder for minimization
-
-      PseudoJet jiggle(const PseudoJet& axis);
-   
-   public:
-      AxesFinderFromKmeansMinimization(AxesFinder *startingFinder, double beta, double Rcutoff, int n_iterations) :
-         AxesFinder(startingFinder),
-         _n_iterations(n_iterations),
-         _noise_range(1.0), // hard coded for the time being
-         _measureFunction(beta, Rcutoff),
-         _onePassFinder(NULL, beta, Rcutoff)
-         {}
-
-      virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes);
-
+public:
+   AxesFinderFromKmeansMinimization(double beta, double Rcutoff, int n_iterations)
+   :  _n_iterations(n_iterations),
+      _noise_range(1.0), // hard coded for the time being
+      _measureFunction(beta, Rcutoff),
+      _onePassFinder(beta, Rcutoff)
+      {}
+
+   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes) const;
+   
+private:
+   int _n_iterations;   // Number of iterations to run  (0 for no minimization, 1 for one-pass, >>1 for global minimum)
+   double _noise_range; // noise range for random initialization
+   
+   DefaultUnnormalizedMeasureFunction _measureFunction; //function to test whether minimum is reached
+   
+   AxesFinderFromOnePassMinimization _onePassFinder;  //one pass finder that is repeatedly called
+   
+   PseudoJet jiggle(const PseudoJet& axis) const;
 };
 
@@ -317 +283 @@
 class AxesFinderFromGeometricMinimization : public AxesFinder {
 
-   private:
-      MeasureFunction* _function;
-      double _Rcutoff;
-      double _nAttempts;
-      double _accuracy;
-
-   
-   public:
-      AxesFinderFromGeometricMinimization(AxesFinder* startingFinder, double beta, double Rcutoff) : AxesFinder(startingFinder), _Rcutoff(Rcutoff) {
-         if (beta != 2.0) {
-            std::cerr << "Geometric minimization is currently only defined for beta = 2.0." << std::endl;
-            exit(1);
-         }
-         
-         _nAttempts = 100;
-         _accuracy = 0.000000001;
-         _function = new GeometricMeasure(beta,_Rcutoff);
+public:
+   AxesFinderFromGeometricMinimization(double beta, double Rcutoff)
+   :  _nAttempts(100),
+      _accuracy(0.000000001),
+      _function(beta,Rcutoff)
+   {
+      if (beta != 2.0) {
+         throw Error("Geometric minimization is currently only defined for beta = 2.0.");
       }
-
-      ~AxesFinderFromGeometricMinimization() {
-         delete _function;
-      }
-   
-      virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes);
+   }
+
+   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes) const;
+
+private:
+   double _nAttempts;
+   double _accuracy;
+   GeometricMeasureFunction _function;
+
+
 };
 
@@ -348 +309 @@
 // in order to better facilitate calculations.
 class LightLikeAxis {
-private:
-   double _rap, _phi, _weight, _mom;
-   
-   double DistanceSq(double rap2, double phi2) const {
-      double rap1 = _rap;
-      double phi1 = _phi;
-      
-      double distRap = rap1-rap2;
-      double distPhi = std::fabs(phi1-phi2);
-      if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
-      return sq(distRap) + sq(distPhi);
-   }
-   
-   double Distance(double rap2, double phi2) const {
-      return std::sqrt(DistanceSq(rap2,phi2));
-   }
-   
-   
+
 public:
    LightLikeAxis() : _rap(0.0), _phi(0.0), _weight(0.0), _mom(0.0) {}
@@ -399 +343 @@
    }
 
+private:
+   double _rap, _phi, _weight, _mom;
+   
+   double DistanceSq(double rap2, double phi2) const {
+      double rap1 = _rap;
+      double phi1 = _phi;
+      
+      double distRap = rap1-rap2;
+      double distPhi = std::fabs(phi1-phi2);
+      if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
+      return distRap*distRap + distPhi*distPhi;
+   }
+   
+   double Distance(double rap2, double phi2) const {
+      return std::sqrt(DistanceSq(rap2,phi2));
+   }
+   
+   
 };