Changeset 35cdc46 in git for external/fastjet/contribs/Nsubjettiness

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

Location:

external/fastjet/contribs/Nsubjettiness

Files:

• AUTHORS (modified) (1 diff)
• AxesFinder.cc (modified) (10 diffs)
• AxesFinder.hh (modified) (15 diffs)
• COPYING (added)
• ChangeLog (added)
• Makefile (added)
• MeasureFunction.cc (modified) (3 diffs)
• MeasureFunction.hh (modified) (5 diffs)
• NEWS (added)
• Njettiness.cc (modified) (7 diffs)
• Njettiness.hh (modified) (5 diffs)
• NjettinessDefinition.cc (added)
• NjettinessDefinition.hh (added)
• NjettinessPlugin.cc (modified) (3 diffs)
• NjettinessPlugin.hh (modified) (9 diffs)
• Nsubjettiness.cc (modified) (1 diff)
• Nsubjettiness.hh (modified) (7 diffs)
• README (added)
• VERSION (modified) (1 diff)
• WinnerTakeAllRecombiner.cc (modified) (2 diffs)
• WinnerTakeAllRecombiner.hh (modified) (2 diffs)

external/fastjet/contribs/Nsubjettiness/AUTHORS

@@ -18 +18 @@
    JHEP 1202:093 (2012), arXiv:1108.2701.
 
+New in v2.0 is the winner-take-all axis, described in:
+
+   Jet Shapes with the Broadening Axis. 
+   Andrew J. Larkoski, Duff Neill, and Jesse Thaler. 
+   JHEP 1404:017 (2014), arXiv:1401.2158.
+
+as well as in unpublished work by Gavin Salam.
+
 ----------------------------------------------------------------------

external/fastjet/contribs/Nsubjettiness/AxesFinder.cc

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: AxesFinder.cc 670 2014-06-06 01:24:42Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -37 +38 @@
 template <int N>
 std::vector<LightLikeAxis> AxesFinderFromOnePassMinimization::UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes, 
-                                  const std::vector <fastjet::PseudoJet> & inputJets) {
+                                  const std::vector <fastjet::PseudoJet> & inputJets) const {
    assert(old_axes.size() == N);
    
@@ -45 +46 @@
    for (int n = 0; n < N; ++n) {
       new_axes[n].reset(0.0,0.0,0.0,0.0);
-#ifdef FASTJET2
-      new_jets[n].reset(0.0,0.0,0.0,0.0);
-#else
-      // use cheaper reset if available
       new_jets[n].reset_momentum(0.0,0.0,0.0,0.0);
-#endif
    }
 
@@ -135 +131 @@
 // (This is just a wrapper for the templated version above.)
 std::vector<LightLikeAxis> AxesFinderFromOnePassMinimization::UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
-                                      const std::vector <fastjet::PseudoJet> & inputJets) {
+                                      const std::vector <fastjet::PseudoJet> & inputJets) const {
    int N = old_axes.size();
    switch (N) {
@@ -166 +162 @@
 // uses minimization of N-jettiness to continually update axes until convergence.
 // The function returns the axes found at the (local) minimum
-std::vector<fastjet::PseudoJet> AxesFinderFromOnePassMinimization::getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& seedAxes) {
+std::vector<fastjet::PseudoJet> AxesFinderFromOnePassMinimization::getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& seedAxes) const {
           
    // convert from PseudoJets to LightLikeAxes
@@ -211 +207 @@
 }
 
-PseudoJet AxesFinderFromKmeansMinimization::jiggle(const PseudoJet& axis) {
+PseudoJet AxesFinderFromKmeansMinimization::jiggle(const PseudoJet& axis) const {
    double phi_noise = ((double)rand()/(double)RAND_MAX) * _noise_range * 2.0 - _noise_range;
    double rap_noise = ((double)rand()/(double)RAND_MAX) * _noise_range * 2.0 - _noise_range;
@@ -226 +222 @@
    
 // Repeatedly calls the one pass finder to try to find global minimum
-std::vector<fastjet::PseudoJet> AxesFinderFromKmeansMinimization::getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& seedAxes) {
+std::vector<fastjet::PseudoJet> AxesFinderFromKmeansMinimization::getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& seedAxes) const {
    
    // first iteration
@@ -256 +252 @@
 // It continually updates until it reaches convergence or it reaches the maximum number of attempts.
 // This is essentially the same as a stable cone finder.
-std::vector<fastjet::PseudoJet> AxesFinderFromGeometricMinimization::getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes) {
+std::vector<fastjet::PseudoJet> AxesFinderFromGeometricMinimization::getAxes(int /*n_jets*/, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes) const {
 
    std::vector<fastjet::PseudoJet> seedAxes = currentAxes;
-   double seedTau = _function->tau(particles, seedAxes);
+   double seedTau = _function.tau(particles, seedAxes);
    
    for (int i = 0; i < _nAttempts; i++) {
@@ -270 +266 @@
          // start from unclustered beam measure
          int minJ = -1;
-         double minDist = _function->beam_distance_squared(particles[i]);
+         double minDist = _function.beam_distance_squared(particles[i]);
          
          // which axis am I closest to?
          for (unsigned int j = 0; j < seedAxes.size(); j++) {
-            double tempDist = _function->jet_distance_squared(particles[i],seedAxes[j]);
+            double tempDist = _function.jet_distance_squared(particles[i],seedAxes[j]);
             if (tempDist < minDist) {
                minDist = tempDist;
@@ -287 +283 @@
       // calculate tau on new axes
       seedAxes = newAxes;
-      double tempTau = _function->tau(particles, newAxes);
+      double tempTau = _function.tau(particles, newAxes);
       
       // close enough to stop?

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));
+   }
+   
+   
 };
 

external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: MeasureFunction.cc 670 2014-06-06 01:24:42Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -36 +37 @@
 
 // Return all of the necessary TauComponents for specific input particles and axes
-TauComponents MeasureFunction::result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) {
+TauComponents MeasureFunction::result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
+   
+   // first find partition
+   // this sets jetPartitionStorage and beamPartitionStorage
+   PseudoJet beamPartitionStorage;
+   std::vector<fastjet::PseudoJet> jetPartitionStorage = get_partition(particles,axes,&beamPartitionStorage);
+   
+   // then return result calculated from partition
+   return result_from_partition(jetPartitionStorage,axes,&beamPartitionStorage);
+}
 
-   std::vector<double> jetPieces(axes.size(), 0.0);
-   double beamPiece = 0.0;
+std::vector<fastjet::PseudoJet> MeasureFunction::get_partition(const std::vector<fastjet::PseudoJet>& particles,
+                                                               const std::vector<fastjet::PseudoJet>& axes,
+                                                               PseudoJet * beamPartitionStorage) const {
    
-   // Calculates the unnormalized sub-tau values, i.e. a std::vector of the contributions to tau_N of each Voronoi region (or region within R_0)
+   std::vector<std::vector<PseudoJet> > jetPartition(axes.size());
+   std::vector<PseudoJet> beamPartition;
+   
+   // Figures out the partiting of the input particles into the various jet pieces
+   // Based on which axis the parition is closest to
    for (unsigned i = 0; i < particles.size(); i++) {
       
@@ -60 +75 @@
       
       if (j_min == -1) {
-         if (_has_beam) beamPiece += beam_numerator(particles[i]);
+         if (_has_beam) beamPartition.push_back(particles[i]);
          else assert(_has_beam);  // this should never happen.
       } else {
-         jetPieces[j_min] += jet_numerator(particles[i],axes[j_min]);
+         jetPartition[j_min].push_back(particles[i]);
       }
    }
    
-   // Calculates normalization for tau and subTau if _has_denominator is true, otherwise returns 1.0 (i.e. no normalization)
-   double tauDen = 0.0;
-   if (_has_denominator) {
-      for (unsigned i = 0; i < particles.size(); i++) {
-         tauDen += denominator(particles[i]);
-      }
-   } else {
-      tauDen = 1.0; // if no denominator, then 1.0 for no normalization factor
+   // Store beam partition
+   if (beamPartitionStorage) {
+      *beamPartitionStorage = join(beamPartition);
+   }
+
+   // Store jet partitions
+   std::vector<PseudoJet> jetPartitionStorage(axes.size(),PseudoJet(0,0,0,0));
+   for (unsigned j = 0; j < axes.size(); j++) {
+      jetPartitionStorage[j] = join(jetPartition[j]);
    }
    
+   return jetPartitionStorage;
+}
+
+// does partition, but only stores index of PseudoJets
+std::vector<std::list<int> > MeasureFunction::get_partition_list(const std::vector<fastjet::PseudoJet>& particles,
+                                                                 const std::vector<fastjet::PseudoJet>& axes) const {
+
+   std::vector<std::list<int> > jetPartition(axes.size());
+   
+   // Figures out the partiting of the input particles into the various jet pieces
+   // Based on which axis the parition is closest to
+   for (unsigned i = 0; i < particles.size(); i++) {
+      
+      // find minimum distance; start with beam (-1) for reference
+      int j_min = -1;
+      double minRsq;
+      if (_has_beam) minRsq = beam_distance_squared(particles[i]);
+      else minRsq = std::numeric_limits<double>::max(); // make it large value
+      
+      // check to see which axis the particle is closest to
+      for (unsigned j = 0; j < axes.size(); j++) {
+         double tempRsq = jet_distance_squared(particles[i],axes[j]); // delta R distance
+         if (tempRsq < minRsq) {
+            minRsq = tempRsq;
+            j_min = j;
+         }
+      }
+      
+      if (j_min == -1) {
+         assert(_has_beam); // consistency check
+      } else {
+         jetPartition[j_min].push_back(i);
+      }
+   }
+   
+   return jetPartition;
+}
+   
+
+// Uses existing partition and calculates result
+// TODO:  Can we cache this for speed up when doing area subtraction?
+TauComponents MeasureFunction::result_from_partition(const std::vector<fastjet::PseudoJet>& jet_partition,
+                                                     const std::vector<fastjet::PseudoJet>& axes,
+                                                     PseudoJet * beamPartitionStorage) const {
+   
+   std::vector<double> jetPieces(axes.size(), 0.0);
+   double beamPiece = 0.0;
+   
+   double tauDen = 0.0;
+   if (!_has_denominator) tauDen = 1.0;  // if no denominator, then 1.0 for no normalization factor
+   
+   // first find jet pieces
+   for (unsigned j = 0; j < axes.size(); j++) {
+      std::vector<PseudoJet> thisPartition = jet_partition[j].constituents();
+      for (unsigned i = 0; i < thisPartition.size(); i++) {
+         jetPieces[j] += jet_numerator(thisPartition[i],axes[j]); //numerator jet piece
+         if (_has_denominator) tauDen += denominator(thisPartition[i]); // denominator
+      }
+   }
+   
+   // then find beam piece
+   if (_has_beam) {
+      assert(beamPartitionStorage); // make sure I have beam information
+      std::vector<PseudoJet> beamPartition = beamPartitionStorage->constituents();
+
+      for (unsigned i = 0; i < beamPartition.size(); i++) {
+         beamPiece += beam_numerator(beamPartition[i]); //numerator beam piece
+         if (_has_denominator) tauDen += denominator(beamPartition[i]); // denominator
+      }
+   }
    return TauComponents(jetPieces, beamPiece, tauDen, _has_denominator, _has_beam);
 }
+
+   
+   
+   
    
 } //namespace contrib

external/fastjet/contribs/Nsubjettiness/MeasureFunction.hh

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: MeasureFunction.hh 678 2014-06-12 20:43:03Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -31 +32 @@
 #include <limits>
 
+
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
 namespace contrib{
 
-inline double sq(double x) {return x*x;}
-
-///////
-//
-// Measure Function
-//
-///////
-
+///////
+//
+// TauComponents
+// (eventually we might want to put this in a separate header file)
+//
+///////
+   
 /// \class TauComponents
 // This class creates a wrapper for the various tau/subtau values calculated in Njettiness. This class allows Njettiness access to these variables
 // without ever having to do the calculation itself. It takes in subtau numerators and tau denominator from MeasureFunction
 // and outputs tau numerator, and normalized tau and subtau.
-// TODO:  Consider merging with NjettinessExtras.  Add axes information?
 class TauComponents {
-private:
-   
-   // these values are input in the constructor
-   std::vector<double> _jet_pieces_numerator;
-   double _beam_piece_numerator;
-   double _denominator;
-   bool _has_denominator; //added so that TauComponents knows if denominator is used or not
-   bool _has_beam; //added so that TauComponents knows if beam regions is used or not
-   
-   // these values are derived from above values
-   std::vector<double> _jet_pieces;
-   double _beam_piece;
-   double _numerator;
-   double _tau;
-   
    
 public:
@@ -115 +100 @@
    double beam_piece() const { return _beam_piece; }
    double tau() const { return _tau; }
-   
-};
+
+private:
+   
+   // these values are input in the constructor
+   std::vector<double> _jet_pieces_numerator;
+   double _beam_piece_numerator;
+   double _denominator;
+   bool _has_denominator; //added so that TauComponents knows if denominator is used or not
+   bool _has_beam; //added so that TauComponents knows if beam regions is used or not
+   
+   // these values are derived from above values
+   std::vector<double> _jet_pieces;
+   double _beam_piece;
+   double _numerator;
+   double _tau;
+   
+};
+
+///////
+//
+// Measure Function
+//
+///////
+
 
 //------------------------------------------------------------------------
@@ -125 +132 @@
 class MeasureFunction {
    
+public:
+   //These functions define the measure by which tau_N is calculated,
+   //and are overloaded by the various measures below
+   
+   // Distanes to axes.  These are called many times, so need to be as fast as possible
+   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const = 0;
+   virtual double beam_distance_squared(const fastjet::PseudoJet& particle) const = 0;
+   
+   // The actual measures used in N-(sub)jettiness
+   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const  = 0;
+   virtual double beam_numerator(const fastjet::PseudoJet& particle) const = 0;
+   
+   // a possible normalization factor
+   virtual double denominator(const fastjet::PseudoJet& particle) const = 0;
+   
+   //------
+   // The functions below call the above functions and are not virtual
+   //------
+   
+   // Return all of the necessary TauComponents for specific input particles and axes
+   // Also have optional pointers to get out information about partitioning
+   TauComponents result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;
+
+   // Just getting tau value if that is all that is needed
+   double tau(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
+      return result(particles,axes).tau();
+   }
+   
+   // Create the partitioning and stores internally
+   std::vector<fastjet::PseudoJet> get_partition(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes, PseudoJet * beamPartitionStorage = NULL) const;
+
+   // Essentially same as get_partition, but in the form needed for the jet algorithm
+   std::vector<std::list<int> > get_partition_list(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;
+
+   // calculates the tau result using an existing partition
+   TauComponents result_from_partition(const std::vector<fastjet::PseudoJet>& jet_partitioning, const std::vector<fastjet::PseudoJet>& axes, PseudoJet * beamPartitionStorage = NULL) const;
+
+   // shorthand for squaring
+   static inline double sq(double x) {return x*x;}
+
+   //virtual destructor
+   virtual ~MeasureFunction(){}
+   
 protected:
    //bool set by derived classes to choose whether or not to use the denominator
@@ -133 +183 @@
    MeasureFunction(bool has_denominator = true, bool has_beam = true) : _has_denominator(has_denominator), _has_beam(has_beam) {}
    
-public:
-   virtual ~MeasureFunction(){}
-   
-   //These functions define the measure by which tau_N is calculated,
-   //and are overloaded by the various measures below
-   
-   // Distanes to axes.  These are called many times, so need to be as fast as possible
-   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) = 0;
-   virtual double beam_distance_squared(const fastjet::PseudoJet& particle) = 0;
-   
-   // The actual measures used in N-(sub)jettiness
-   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) = 0;
-   virtual double beam_numerator(const fastjet::PseudoJet& particle) = 0;
-   
-   // a possible normalization factor
-   virtual double denominator(const fastjet::PseudoJet& particle) = 0;
-   
-   
-   // These functions call the above functions and are not virtual
-   
-   // Do I cluster a particle into a jet?
-   bool do_cluster(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
-      return (jet_distance_squared(particle,axis) <= beam_distance_squared(particle));
-   }
-   
-   // Return all of the necessary TauComponents for specific input particles and axes
-   TauComponents result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes);
-
-   double tau(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) {
-      return result(particles,axes).tau();
-   }
-
-   
-};
-
-
-/// \class DefaultNormalizedMeasure
+};
+
+
+/// \class DefaultNormalizedMeasureFunction
 // This class is the default measure, inheriting from the class above. This class will calculate tau_N 
 // of a jet according to this measure. This measure is defined as the pT of the particle multiplied by deltaR 
 // to the power of beta. This class includes the normalization factor determined by R0
-class DefaultNormalizedMeasure : public MeasureFunction {
-
-   private:
-      double _beta;
-      double _R0;
-      double _Rcutoff;
-
-   public:
-
-      DefaultNormalizedMeasure(double beta, double R0, double Rcutoff, bool normalized = true)
-      : MeasureFunction(normalized), _beta(beta), _R0(R0), _Rcutoff(Rcutoff) {}
-
-      virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
-         return particle.squared_distance(axis);
-      }
-   
-      virtual double beam_distance_squared(const fastjet::PseudoJet& particle) {
-         return sq(_Rcutoff);
-      }
-
-      virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
-         return particle.perp() * std::pow(jet_distance_squared(particle,axis),_beta/2.0);
-      }
-   
-      virtual double beam_numerator(const fastjet::PseudoJet& particle) {
-         return particle.perp() * std::pow(_Rcutoff,_beta);
-      }
-
-      virtual double denominator(const fastjet::PseudoJet& particle) {
-         return particle.perp() * std::pow(_R0,_beta);
-      }
-
+class DefaultNormalizedMeasureFunction : public MeasureFunction {
+
+public:
+
+   DefaultNormalizedMeasureFunction(double beta, double R0, double Rcutoff, bool normalized = true)
+   : MeasureFunction(normalized), _beta(beta), _R0(R0), _Rcutoff(Rcutoff) {}
+
+   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      return particle.squared_distance(axis);
+   }
+
+   virtual double beam_distance_squared(const fastjet::PseudoJet& /*particle*/) const {
+      return sq(_Rcutoff);
+   }
+
+   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const{
+      return particle.perp() * std::pow(jet_distance_squared(particle,axis),_beta/2.0);
+   }
+
+   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
+      return particle.perp() * std::pow(_Rcutoff,_beta);
+   }
+
+   virtual double denominator(const fastjet::PseudoJet& particle) const {
+      return particle.perp() * std::pow(_R0,_beta);
+   }
+   
+private:
+   double _beta;
+   double _R0;
+   double _Rcutoff;
+
+   
 };
 
 //------------------------------------------------------------------------
-/// \class DefaultUnnormalizedMeasure
+/// \class DefaultUnnormalizedMeasureFunction
 // This class is the unnormalized default measure, inheriting from the class above. The only difference from above
 // is that the denominator is defined to be 1.0 by setting _has_denominator to false.
-class DefaultUnnormalizedMeasure : public DefaultNormalizedMeasure {
-
-   public:
-      // Since all methods are identical, UnnormalizedMeasure inherits directly from NormalizedMeasure. R0 is defaulted to NAN since the value of R0 is unecessary for this class.
-      // the "false" flag sets _has_denominator in MeasureFunction to false so no denominator is used.
-      DefaultUnnormalizedMeasure(double beta, double Rcutoff) : DefaultNormalizedMeasure(beta, NAN, Rcutoff, false) {}
-
-      
+class DefaultUnnormalizedMeasureFunction : public DefaultNormalizedMeasureFunction {
+
+public:
+   // Since all methods are identical, UnnormalizedMeasure inherits directly
+   // from NormalizedMeasure. R0 is a dummy value since the value of R0 is unecessary for this class,
+   // and the "false" flag sets _has_denominator in MeasureFunction to false so no denominator is used.
+   DefaultUnnormalizedMeasureFunction(double beta, double Rcutoff)
+   : DefaultNormalizedMeasureFunction(beta, std::numeric_limits<double>::quiet_NaN(), Rcutoff, false) {}
 };
 
 //------------------------------------------------------------------------
-/// \class GeometricMeasure
+/// \class GeometricMeasureFunction
 // This class is the geometic measure, inheriting from the class above. This class will calculate tau_N 
 // of a jet according to this measure. This measure is defined by the Lorentz dot product between
 // the particle and the axis. This class includes normalization of tau_N.
-class GeometricMeasure : public MeasureFunction {
-
-   private:
-      double _jet_beta;
-      double _beam_beta;
-      double _Rcutoff;
-
-      // create light-like axis
-      fastjet::PseudoJet lightFrom(const fastjet::PseudoJet& input) const {
-         double length = sqrt(pow(input.px(),2) + pow(input.py(),2) + pow(input.pz(),2));
-         return fastjet::PseudoJet(input.px()/length,input.py()/length,input.pz()/length,1.0);
-      }
-
-   public:
-      // Right now, we are hard coded for beam_beta = 1.0, but that will need to change
-      GeometricMeasure(double jet_beta, double Rcutoff) : _jet_beta(jet_beta), _beam_beta(1.0), _Rcutoff(Rcutoff) {}
-   
-      virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
-         fastjet::PseudoJet lightAxis = lightFrom(axis);
-         double pseudoRsquared = 2.0*dot_product(lightFrom(axis),particle)/(lightAxis.pt()*particle.pt());
-         return pseudoRsquared;
-      }
-   
-      virtual double beam_distance_squared(const fastjet::PseudoJet& particle) {
-         return sq(_Rcutoff);
-      }
-
-      virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) {
-         fastjet::PseudoJet lightAxis = lightFrom(axis);
-         double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(lightAxis.pt(),_beam_beta - 1.0);
-         return particle.pt() * weight * std::pow(jet_distance_squared(particle,axis),_jet_beta/2.0);
-      }
-   
-      virtual double beam_numerator(const fastjet::PseudoJet& particle) {
-         double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(particle.pt()/particle.e(),_beam_beta - 1.0);
-         return particle.pt() * weight * std::pow(_Rcutoff,_jet_beta);
-      }
-
-      virtual double denominator(const fastjet::PseudoJet& particle) {
-         return 1.0;
-      }
-};
-
-
+class GeometricMeasureFunction : public MeasureFunction {
+
+public:
+   // Right now, we are hard coded for beam_beta = 1.0, but that will need to change
+   GeometricMeasureFunction(double jet_beta, double Rcutoff) : 
+     MeasureFunction(false), // doesn't have denominator
+     _jet_beta(jet_beta), _beam_beta(1.0), _Rcutoff(Rcutoff) {}
+
+   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      fastjet::PseudoJet lightAxis = lightFrom(axis);
+      double pseudoRsquared = 2.0*dot_product(lightFrom(axis),particle)/(lightAxis.pt()*particle.pt());
+      return pseudoRsquared;
+   }
+
+   virtual double beam_distance_squared(const fastjet::PseudoJet&  /*particle*/) const {
+      return sq(_Rcutoff);
+   }
+
+   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      fastjet::PseudoJet lightAxis = lightFrom(axis);
+      double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(lightAxis.pt(),_beam_beta - 1.0);
+      return particle.pt() * weight * std::pow(jet_distance_squared(particle,axis),_jet_beta/2.0);
+   }
+
+   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
+      double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(particle.pt()/particle.e(),_beam_beta - 1.0);
+      return particle.pt() * weight * std::pow(_Rcutoff,_jet_beta);
+   }
+
+   virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
+      return std::numeric_limits<double>::quiet_NaN();
+   }
+   
+   
+private:
+   double _jet_beta;
+   double _beam_beta;
+   double _Rcutoff;
+   
+   // create light-like axis
+   fastjet::PseudoJet lightFrom(const fastjet::PseudoJet& input) const {
+      double length = sqrt(pow(input.px(),2) + pow(input.py(),2) + pow(input.pz(),2));
+      return fastjet::PseudoJet(input.px()/length,input.py()/length,input.pz()/length,1.0);
+   }
+
+};
+   
+   
 } //namespace contrib
 

external/fastjet/contribs/Nsubjettiness/Njettiness.cc

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: Njettiness.cc 677 2014-06-12 18:56:46Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -35 +36 @@
 ///////
 
-// Helper function to correlate one pass minimization with appropriate measure
-void Njettiness::setOnePassAxesFinder(MeasureMode measure_mode, AxesFinder* startingFinder, double beta, double Rcutoff) {
-   if (measure_mode == normalized_measure || measure_mode == unnormalized_measure || measure_mode == normalized_cutoff_measure || measure_mode == unnormalized_cutoff_measure) {
-      _axesFinder = new AxesFinderFromOnePassMinimization(startingFinder, beta, Rcutoff);
-   }
-   else if (measure_mode == geometric_measure || measure_mode == geometric_cutoff_measure) {
-      _axesFinder = new AxesFinderFromGeometricMinimization(startingFinder, beta, Rcutoff);
-   }
-   else {
-      std::cerr << "Minimization only set up for normalized_measure, unnormalized_measure, normalized_cutoff_measure, unnormalized_cutoff_measure, geometric_measure, geometric_cutoff_measure" << std::endl;
-      exit(1); }
-}
-
-// Parsing needed for constructor to set AxesFinder and MeasureFunction
-// All of the parameter handling is here, and checking that number of parameters is correct.
-void Njettiness::setMeasureFunctionandAxesFinder(AxesMode axes_mode, MeasureMode measure_mode, double para1, double para2, double para3, double para4) {
+Njettiness::Njettiness(const AxesDefinition & axes_def, const MeasureDefinition & measure_def)
+: _axes_def(axes_def.create()), _measure_def(measure_def.create()) {
+   setMeasureFunctionAndAxesFinder();  // call helper function to do the hard work
+}
+
+Njettiness::Njettiness(AxesMode axes_mode, const MeasureDefinition & measure_def)
+: _axes_def(createAxesDef(axes_mode)), _measure_def(measure_def.create()) {
+   setMeasureFunctionAndAxesFinder();  // call helper function to do the hard work
+}
+   
+// Convert from MeasureMode enum to MeasureDefinition
+// This returns a pointer that will be claimed by a SharedPtr
+MeasureDefinition* Njettiness::createMeasureDef(MeasureMode measure_mode, int num_para, double para1, double para2, double para3) const {
 
    // definition of maximum Rcutoff for non-cutoff measures, changed later by other measures
    double Rcutoff = std::numeric_limits<double>::max();  //large number
-   // Most (but all measures have some kind of beta value)
-   double beta = NAN;
+   // Most (but not all) measures have some kind of beta value
+   double beta = std::numeric_limits<double>::quiet_NaN();
    // The normalized measures have an R0 value.
-   double R0 = NAN;
-
+   double R0 = std::numeric_limits<double>::quiet_NaN();
+   
    // Find the MeasureFunction and set the parameters.
    switch (measure_mode) {
@@ -64 +62 @@
          beta = para1;
          R0 = para2;
-         if(correctParameterCount(2, para1, para2, para3, para4)) 
-            _measureFunction = new DefaultNormalizedMeasure(beta, R0, Rcutoff); //normalized_measure requires 2 parameters, beta and R0
-         else { 
-            std::cerr << "normalized_measure needs 2 parameters (beta and R0)" << std::endl;
-            exit(1); }
+         if(num_para == 2) {
+            return new NormalizedMeasure(beta,R0);
+         } else {
+            throw Error("normalized_measure needs 2 parameters (beta and R0)");
+         }
          break;
       case unnormalized_measure:
          beta = para1;
-         if(correctParameterCount(1, para1, para2, para3, para4)) 
-            _measureFunction = new DefaultUnnormalizedMeasure(beta, Rcutoff); //unnormalized_measure requires 1 parameter, beta
-         else {
-            std::cerr << "unnormalized_measure needs 1 parameter (beta)" << std::endl;
-            exit(1); }
+         if(num_para == 1) {
+            return new UnnormalizedMeasure(beta);
+         } else {
+            throw Error("unnormalized_measure needs 1 parameter (beta)");
+         }
          break;
       case geometric_measure:
          beta = para1;
-         if(correctParameterCount(1, para1, para2, para3, para4))
-            _measureFunction = new GeometricMeasure(beta,Rcutoff); //geometric_measure requires 1 parameter, beta
-         else {
-            std::cerr << "geometric_measure needs 1 parameter (beta)" << std::endl;
-            exit(1); }
+         if (num_para == 1) {
+            return new GeometricMeasure(beta);
+         } else {
+            throw Error("geometric_measure needs 1 parameter (beta)");
+         }
          break;
       case normalized_cutoff_measure:
@@ -90 +88 @@
          R0 = para2;
          Rcutoff = para3; //Rcutoff parameter is 3rd parameter in normalized_cutoff_measure
-         if(correctParameterCount(3, para1, para2, para3, para4))
-            _measureFunction = new DefaultNormalizedMeasure(beta, R0, Rcutoff); //normalized_cutoff_measure requires 3 parameters, beta, R0, and Rcutoff
-         else { 
-            std::cerr << "normalized_cutoff_measure has 3 parameters (beta, R0, Rcutoff)" << std::endl;
-            exit(1); }
+         if (num_para == 3) {
+            return new NormalizedCutoffMeasure(beta,R0,Rcutoff);
+         } else {
+            throw Error("normalized_cutoff_measure has 3 parameters (beta, R0, Rcutoff)");
+         }
          break;
       case unnormalized_cutoff_measure:
          beta = para1;
          Rcutoff = para2; //Rcutoff parameter is 2nd parameter in normalized_cutoff_measure
-         if (correctParameterCount(2, para1, para2, para3, para4))
-            _measureFunction = new DefaultUnnormalizedMeasure(beta, Rcutoff); //unnormalized_cutoff_measure requires 2 parameters, beta and Rcutoff
-         else {
-            std::cerr << "unnormalized_cutoff_measure has 2 parameters (beta, Rcutoff)" << std::endl;
-            exit(1); }
+         if (num_para == 2) {
+            return new UnnormalizedCutoffMeasure(beta,Rcutoff);
+         } else {
+            throw Error("unnormalized_cutoff_measure has 2 parameters (beta, Rcutoff)");
+         }
          break;
       case geometric_cutoff_measure:
          beta = para1;
          Rcutoff = para2; //Rcutoff parameter is 2nd parameter in geometric_cutoff_measure
-         if(correctParameterCount(2, para1, para2, para3, para4))
-            _measureFunction = new GeometricMeasure(beta,Rcutoff); //geometric_cutoff_measure requires 2 parameters, beta and Rcutoff
-         else {
-            std::cerr << "geometric_cutoff_measure has 2 parameters (beta,Rcutoff)" << std::endl;
-            exit(1); }
+         if(num_para == 2) {
+           return new GeometricCutoffMeasure(beta,Rcutoff);
+         } else {
+            throw Error("geometric_cutoff_measure has 2 parameters (beta, Rcutoff)");
+         }
          break;
       default:
          assert(false);
          break;
-   }   
-
-   // Choose which AxesFinder from user input.
-   // Uses setOnePassAxesFinder helpful function to use beta and Rcutoff values about (if needed)
+   }
+   return NULL;
+}
+
+// Convert from AxesMode enum to AxesDefinition
+// This returns a pointer that will be claimed by a SharedPtr
+AxesDefinition* Njettiness::createAxesDef(Njettiness::AxesMode axes_mode) const {
+   
    switch (axes_mode) {
       case wta_kt_axes:
-         _axesFinder = new AxesFinderFromWTA_KT(); 
-         break;
+         return new WTA_KT_Axes();
       case wta_ca_axes:
-         _axesFinder = new AxesFinderFromWTA_CA(); 
-         break;
+         return new WTA_CA_Axes();
       case kt_axes:
-         _axesFinder = new AxesFinderFromKT();
-         break;
+         return new KT_Axes();
       case ca_axes:
-         _axesFinder = new AxesFinderFromCA();
-         break;
+         return new CA_Axes();
       case antikt_0p2_axes:
-         _axesFinder = new AxesFinderFromAntiKT(0.2);     
-         break;
+         return new AntiKT_Axes(0.2);
       case onepass_wta_kt_axes:
-         setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA_KT(), beta, Rcutoff);
-         break;
+         return new OnePass_WTA_KT_Axes();
       case onepass_wta_ca_axes:
-         setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA_CA(), beta, Rcutoff);
-         break;
+         return new OnePass_WTA_CA_Axes();
       case onepass_kt_axes:
-         setOnePassAxesFinder(measure_mode, new AxesFinderFromKT(), beta, Rcutoff);
-         break;
+         return new OnePass_KT_Axes();
       case onepass_ca_axes:
-         setOnePassAxesFinder(measure_mode, new AxesFinderFromCA(), beta, Rcutoff);
-         break;
+         return new OnePass_CA_Axes();
       case onepass_antikt_0p2_axes:
-         setOnePassAxesFinder(measure_mode, new AxesFinderFromAntiKT(0.2), beta, Rcutoff);
-         break;
+         return new OnePass_AntiKT_Axes(0.2);
       case onepass_manual_axes:
-         setOnePassAxesFinder(measure_mode, new AxesFinderFromUserInput(), beta, Rcutoff);
-         break;
-      case min_axes: //full minimization is not defined for geometric_measure.
-         if (measure_mode == normalized_measure || measure_mode == unnormalized_measure || measure_mode == normalized_cutoff_measure || measure_mode == unnormalized_cutoff_measure)
-            //Defaults to 100 iteration to find minimum
-            _axesFinder = new AxesFinderFromKmeansMinimization(new AxesFinderFromKT(), beta, Rcutoff, 100);
-         else {
-            std::cerr << "Multi-pass minimization only set up for normalized_measure, unnormalized_measure, normalized_cutoff_measure, unnormalized_cutoff_measure." << std::endl;
-            exit(1);
-         }
-         break;
+         return new OnePass_Manual_Axes();
+      case min_axes:
+         return new MultiPass_Axes(100);
       case manual_axes:
-         _axesFinder = new AxesFinderFromUserInput();
-         break;
-// These options have been commented out because they have not been fully tested
-//      case wta2_kt_axes: // option for alpha = 2 added
-//         _axesFinder = new AxesFinderFromWTA2_KT();
-//         break;
-//      case wta2_ca_axes: // option for alpha = 2 added
-//         _axesFinder = new AxesFinderFromWTA2_CA();
-//         break;
-//      case onepass_wta2_kt_axes: // option for alpha = 2 added
-//         setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA2_KT(), beta, Rcutoff);
-//         break;
-//      case onepass_wta2_ca_axes: // option for alpha = 2 added
-//         setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA2_CA(), beta, Rcutoff);
-//         break;
+         return new Manual_Axes();
       default:
          assert(false);
-         break;
-      }   
-
+         return NULL;
+   }
+}
+
+   
+// Parsing needed for constructor to set AxesFinder and MeasureFunction
+// All of the parameter handling is here, and checking that number of parameters is correct.
+void Njettiness::setMeasureFunctionAndAxesFinder() {
+   // Get the correct MeasureFunction and AxesFinders
+   _measureFunction.reset(_measure_def->createMeasureFunction());
+   _startingAxesFinder.reset(_axes_def->createStartingAxesFinder(*_measure_def));
+   _finishingAxesFinder.reset(_axes_def->createFinishingAxesFinder(*_measure_def));
 }
 
 // setAxes for Manual mode
-void Njettiness::setAxes(std::vector<fastjet::PseudoJet> myAxes) {
-   if (_current_axes_mode == manual_axes || _current_axes_mode == onepass_manual_axes) {
+void Njettiness::setAxes(const std::vector<fastjet::PseudoJet> & myAxes) {
+   if (_axes_def->supportsManualAxes()) {
       _currentAxes = myAxes;
-   }
-   else {
-      std::cerr << "You can only use setAxes if using manual_axes or onepass_manual_axes measure mode" << std::endl;
-      exit(1);
+   } else {
+      throw Error("You can only use setAxes for manual AxesDefinitions");
    }
 }
@@ -200 +177 @@
 // Calculates and returns all TauComponents that user would want.
 // This information is stored in _current_tau_components for later access as well.
-TauComponents Njettiness::getTauComponents(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) {
+TauComponents Njettiness::getTauComponents(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) const {
    if (inputJets.size() <= n_jets) {  //if not enough particles, return zero
       _currentAxes = inputJets;
@@ -206 +183 @@
       _current_tau_components = TauComponents();
       _seedAxes = _currentAxes;
+      _currentJets = _currentAxes;
+      _currentBeam = PseudoJet(0.0,0.0,0.0,0.0);
    } else {
-      _currentAxes = _axesFinder->getAxes(n_jets,inputJets,_currentAxes); // sets current Axes
-      _seedAxes = _axesFinder->seedAxes(); // sets seed Axes (if one pass minimization was used)
-      _current_tau_components = _measureFunction->result(inputJets, _currentAxes);  // sets current Tau Values
+
+      _seedAxes = _startingAxesFinder->getAxes(n_jets,inputJets,_currentAxes); //sets starting point for minimization
+      if (_finishingAxesFinder) {
+         _currentAxes = _finishingAxesFinder->getAxes(n_jets,inputJets,_seedAxes);
+      } else {
+         _currentAxes = _seedAxes;
+      }
+      
+      // Find partition and store information
+      // (jet information in _currentJets, beam in _currentBeam)
+      _currentJets = _measureFunction->get_partition(inputJets,_currentAxes,&_currentBeam);
+      
+      // Find tau value and store information
+      _current_tau_components = _measureFunction->result_from_partition(_currentJets, _currentAxes,&_currentBeam);  // sets current Tau Values
    }
    return _current_tau_components;
@@ -219 +209 @@
 // Each vector element is a list of ints corresponding to the indices in
 // particles of the particles belonging to that jet.
-// TODO:  Consider moving to MeasureFunction
-std::vector<std::list<int> > Njettiness::getPartition(const std::vector<fastjet::PseudoJet> & particles) {
-   std::vector<std::list<int> > partitions(_currentAxes.size());
-
-   for (unsigned i = 0; i < particles.size(); i++) {
-      
-      int j_min = -1;
-      // find minimum distance
-      double minR = std::numeric_limits<double>::max();  //large number
-      for (unsigned j = 0; j < _currentAxes.size(); j++) {
-         double tempR = _measureFunction->jet_distance_squared(particles[i],_currentAxes[j]); // delta R distance
-         if (tempR < minR) {
-            minR = tempR;
-            j_min = j;
-         }
-      }
-      if (_measureFunction->do_cluster(particles[i],_currentAxes[j_min])) partitions[j_min].push_back(i);
-   }
-   return partitions;
-}
-
-// Having found axes, assign each particle in particles to an axis, and return a set of jets.
-// Each jet is the sum of particles closest to an axis (Njet = Naxes).
-// TODO:  Consider moving to MeasureFunction
-std::vector<fastjet::PseudoJet> Njettiness::getJets(const std::vector<fastjet::PseudoJet> & particles) {
-   
-   std::vector<fastjet::PseudoJet> jets(_currentAxes.size());
-
-   std::vector<std::list<int> > partition = getPartition(particles);
-   for (unsigned j = 0; j < partition.size(); ++j) {
-      std::list<int>::const_iterator it, itE;
-      for (it = partition[j].begin(), itE = partition[j].end(); it != itE; ++it) {
-         jets[j] += particles[*it];
-      }
-   }
-   return jets;
-}
-
+std::vector<std::list<int> > Njettiness::getPartitionList(const std::vector<fastjet::PseudoJet> & particles) const {
+   // core code is in MeasureFunction
+   return _measureFunction->get_partition_list(particles,_currentAxes);
+}
+
+   
 } // namespace contrib
 

external/fastjet/contribs/Nsubjettiness/Njettiness.hh

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: Njettiness.hh 670 2014-06-06 01:24:42Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -25 +26 @@
 #define __FASTJET_CONTRIB_NJETTINESS_HH__
 
+
 #include "MeasureFunction.hh"
 #include "AxesFinder.hh"
+#include "NjettinessDefinition.hh"
 
 #include "fastjet/PseudoJet.hh"
+#include "fastjet/SharedPtr.hh"
+#include <fastjet/LimitedWarning.hh>
+
 #include <cmath>
 #include <vector>
 #include <list>
 
-
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
 namespace contrib {
-
+   
 ///////
 //
@@ -54 +59 @@
    
    // The various axes choices available to the user
+   // It is recommended to use AxesDefinition instead of these.
    enum AxesMode {
       kt_axes,             // exclusive kt axes
@@ -78 +84 @@
    // "normalized_cutoff_measure" was the default in v1.0 of Nsubjettiness
    // "unnormalized_measure" is now the recommended default usage
+   // But it is recommended to use MeasureDefinition instead of these.
    enum MeasureMode {
       normalized_measure,           //default normalized measure
@@ -87 +94 @@
    };
 
-private:
-   // The chosen axes/measure modes
-   AxesFinder* _axesFinder;  // The chosen axes
-   MeasureFunction* _measureFunction; // The chosen measure
+   // Main constructor that uses AxesMode and MeasureDefinition to specify measure
+   // Unlike Nsubjettiness or NjettinessPlugin, the value N is not chosen
+   Njettiness(const AxesDefinition & axes_def, const MeasureDefinition & measure_def);
 
-   // Enum information so functions can specify output based on specific options, primarily for setAxes
-   AxesMode _current_axes_mode;
-   MeasureMode _current_measure_mode;
-   
-   // Information about the current information
-   TauComponents _current_tau_components; //automatically set to have components of 0; these values will be set by the getTau function call
-   std::vector<fastjet::PseudoJet> _currentAxes;
-   std::vector<fastjet::PseudoJet> _seedAxes; // axes used prior to minimization (if applicable)
-   
-   // Needed for compilation of non C++11 users
-   bool isnan(double para) { return para != para; }
+   // Intermediate constructor (needed to enable v1.0.3 backwards compatibility?)
+   Njettiness(AxesMode axes_mode, const MeasureDefinition & measure_def);
 
-   // Helpful function to check to make sure input has correct number of parameters
-   bool correctParameterCount(int n, double para1, double para2, double para3, double para4){
-      int numpara;
-      if (!isnan(para1) && !isnan(para2) && !isnan(para3) && !isnan(para4)) numpara = 4;
-      else if (!isnan(para1) && !isnan(para2) && !isnan(para3) && isnan(para4)) numpara = 3;
-      else if (!isnan(para1) && !isnan(para2) && isnan(para3) && isnan(para4)) numpara = 2;
-      else if (!isnan(para1) && isnan(para2) && isnan(para3) && isnan(para4)) numpara = 1;
-      else numpara = 0;
-      return n == numpara;
+   // Alternative constructor which takes axes/measure information as enums with measure parameters
+   // This version is not recommended
+   Njettiness(AxesMode axes_mode,
+              MeasureMode measure_mode,
+              int num_para,
+              double para1 = std::numeric_limits<double>::quiet_NaN(),
+              double para2 = std::numeric_limits<double>::quiet_NaN(),
+              double para3 = std::numeric_limits<double>::quiet_NaN())
+   : _axes_def(createAxesDef(axes_mode)), _measure_def(createMeasureDef(measure_mode, num_para, para1, para2, para3)) {
+      setMeasureFunctionAndAxesFinder();  // call helper function to do the hard work
    }
 
-   // Helper function to set onepass_axes depending on input measure_mode and startingFinder
-   void setOnePassAxesFinder(MeasureMode measure_mode, AxesFinder* startingFinder, double para1, double Rcutoff);
- 
-   // created separate function to set MeasureFunction and AxesFinder in order to keep constructor cleaner.
-   void setMeasureFunctionandAxesFinder(AxesMode axes_mode, MeasureMode measure_mode, double para1, double para2, double para3, double para4);
-
-public:
-
-
-   // Main constructor which takes axes/measure information, and possible parameters.
-   // Unlike Nsubjettiness or NjettinessPlugin, the value N is not chosen
-   Njettiness(AxesMode axes_mode,
-              MeasureMode measure_mode,
-              double para1 = NAN,
-              double para2 = NAN,
-              double para3 = NAN,
-              double para4 = NAN)
-   : _current_axes_mode(axes_mode),
-   _current_measure_mode(measure_mode) {
-      setMeasureFunctionandAxesFinder(axes_mode, measure_mode, para1, para2, para3, para4);  // call helper function to do the hard work
-   }
-
-   ~Njettiness() {
-      // clean house
-      delete _measureFunction;
-      delete _axesFinder;
-   }
+   // destructor
+   ~Njettiness() {};
    
    // setAxes for Manual mode
-   void setAxes(std::vector<fastjet::PseudoJet> myAxes);
+   void setAxes(const std::vector<fastjet::PseudoJet> & myAxes);
    
    // Calculates and returns all TauComponents that user would want.
    // This information is stored in _current_tau_components for later access as well.
-   TauComponents getTauComponents(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets);
+   TauComponents getTauComponents(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) const;
 
    // Calculates the value of N-subjettiness,
    // but only returns the tau value from _current_tau_components
-   double getTau(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) {
+   double getTau(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) const {
       return getTauComponents(n_jets, inputJets).tau();
    }
-   
-   // returns enum information
-   MeasureMode currentMeasureMode() { return _current_measure_mode;}
-   AxesMode currentAxesMode() { return _current_axes_mode;}
 
    // Return all relevant information about tau components
-   TauComponents currentTauComponents() {return _current_tau_components;}
-   
+   TauComponents currentTauComponents() const {return _current_tau_components;}
    // Return axes found by getTauComponents.
-   std::vector<fastjet::PseudoJet> currentAxes() { return _currentAxes;}
+   std::vector<fastjet::PseudoJet> currentAxes() const { return _currentAxes;}
    // Return seedAxes used if onepass minimization (otherwise, same as currentAxes)
-   std::vector<fastjet::PseudoJet> seedAxes() { return _seedAxes;}
+   std::vector<fastjet::PseudoJet> seedAxes() const { return _seedAxes;}
+   // Return jet partition found by getTauComponents.
+   std::vector<fastjet::PseudoJet> currentJets() const {return _currentJets;}
+   // Return beam partition found by getTauComponents.
+   fastjet::PseudoJet currentBeam() const {return _currentBeam;}
    
    // partition inputs by Voronoi (each vector stores indices corresponding to inputJets)
-   std::vector<std::list<int> > getPartition(const std::vector<fastjet::PseudoJet> & inputJets);
+   std::vector<std::list<int> > getPartitionList(const std::vector<fastjet::PseudoJet> & inputJets) const;
 
-   // partition inputs by Voronoi
-   std::vector<fastjet::PseudoJet> getJets(const std::vector<fastjet::PseudoJet> & inputJets);
+private:
+   
+   // Information about Axes and Measures to be Used
+   // Implemented as SharedPtrs to avoid memory management headaches
+   SharedPtr<const AxesDefinition> _axes_def;
+   SharedPtr<const MeasureDefinition> _measure_def;
+   
+   // The chosen axes/measure mode workers
+   // Implemented as SharedPtrs to avoid memory management headaches
+   // TODO: make into a SharedPtr<const AxesFinder>?
+   SharedPtr<MeasureFunction> _measureFunction;  // The chosen measure
+   SharedPtr<AxesFinder> _startingAxesFinder;    // The initial axes finder
+   SharedPtr<AxesFinder> _finishingAxesFinder;   // A possible minimization step
+   
+   // Information about the current information
+   // Defined as mutables, so user should be aware that these change when getTau is called.
+   mutable TauComponents _current_tau_components; //automatically set to have components of 0; these values will be set by the getTau function call
+   mutable std::vector<fastjet::PseudoJet> _currentAxes; //axes found after minimization
+   mutable std::vector<fastjet::PseudoJet> _seedAxes; // axes used prior to minimization (if applicable)
+   mutable std::vector<fastjet::PseudoJet> _currentJets; //partitioning information
+   mutable fastjet::PseudoJet _currentBeam; //return beam, if requested
+   
+   // created separate function to set MeasureFunction and AxesFinder in order to keep constructor cleaner.
+   void setMeasureFunctionAndAxesFinder();
+   
+   // Convert old style enums into new style MeasureDefinition
+   AxesDefinition* createAxesDef(AxesMode axes_mode) const;
+   
+   // Convert old style enums into new style MeasureDefinition
+   MeasureDefinition* createMeasureDef(MeasureMode measure_mode, int num_para, double para1, double para2, double para3) const;
 
 };
-
+   
 } // namespace contrib
 

external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.cc

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: NjettinessPlugin.cc 663 2014-06-03 21:26:41Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -29 +30 @@
 
 
-// Constructor with same arguments as Nsubjettiness.
-NjettinessPlugin::NjettinessPlugin(int N, Njettiness::AxesMode axes_mode, Njettiness::MeasureMode measure_mode, double para1, double para2, double para3, double para4)
-  : _N(N), _njettinessFinder(axes_mode, measure_mode, para1, para2, para3, para4) {}
-
-// Old constructor for compatibility
-NjettinessPlugin::NjettinessPlugin(int N, Njettiness::AxesMode mode, double beta, double R0, double Rcutoff)
-   : _N(N), _njettinessFinder(mode, Njettiness::normalized_cutoff_measure, beta, R0, Rcutoff) {}
 
 std::string NjettinessPlugin::description() const {return "N-jettiness jet finder";}
 
+
 // Clusters the particles according to the Njettiness jet algorithm
-// TODO: this code should be revisited to see if if can be made more clear.
+// Apologies for the complication with this code, but we need to make
+// a fake jet clustering tree.  The partitioning is done by getPartitionList
 void NjettinessPlugin::run_clustering(ClusterSequence& cs) const
 {
    std::vector<fastjet::PseudoJet> particles = cs.jets();
+
+   // HACK: remove area information from particles (in case this is called by
+   // a ClusterSequenceArea.  Will be fixed in a future FastJet release)
+   for (unsigned i = 0; i < particles.size(); i++) {
+      particles[i].set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>());
+   }
+   
+   
    _njettinessFinder.getTau(_N, particles);
-   std::vector<std::list<int> > partition = _njettinessFinder.getPartition(particles);
+
+   std::vector<std::list<int> > partition = _njettinessFinder.getPartitionList(particles);
 
    std::vector<fastjet::PseudoJet> jet_indices_for_extras;
 
    // output clusterings for each jet
-   for (size_t i = 0; i < partition.size(); ++i) {
+   for (size_t i0 = 0; i0 < partition.size(); ++i0) {
+      size_t i = partition.size() - 1 - i0; // reversed order of reading to match axes order
       std::list<int>& indices = partition[i];
       if (indices.size() == 0) continue;
@@ -70 +76 @@
    }
 
+   //HACK:  Re-reverse order of reading to match CS order
+   reverse(jet_indices_for_extras.begin(),jet_indices_for_extras.end());
+
    NjettinessExtras * extras = new NjettinessExtras(_njettinessFinder.currentTauComponents(),jet_indices_for_extras,_njettinessFinder.currentAxes());
    cs.plugin_associate_extras(std::auto_ptr<ClusterSequence::Extras>(extras));

external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh

@@ -1 @@
-// $Id$
-//
 //  Nsubjettiness Package
 //  Questions/Comments?  jthaler@jthaler.net
@@ -7 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: NjettinessPlugin.hh 671 2014-06-10 17:47:52Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -49 +48 @@
 // to similar information
 class NjettinessExtras : public ClusterSequence::Extras {
-   private:
-   
-      TauComponents _tau_components;
-      std::vector<fastjet::PseudoJet> _jets;
-      std::vector<fastjet::PseudoJet> _axes;
-      
-      int labelOf(const fastjet::PseudoJet& jet) const {
-         int thisJet = -1;
-         for (unsigned int i = 0; i < _jets.size(); i++) {
-            if (_jets[i].cluster_hist_index() == jet.cluster_hist_index()) {
-               thisJet = i;
-               break;
-            }
-         }
-         return thisJet;
-      }
-      
+   
    public:
       NjettinessExtras(TauComponents tau_components, std::vector<fastjet::PseudoJet> jets, std::vector<fastjet::PseudoJet> axes) : _tau_components(tau_components), _jets(jets), _axes(axes) {}
@@ -74 +57 @@
       std::vector<fastjet::PseudoJet> axes() const {return _axes;}
       
-      double totalTau(const fastjet::PseudoJet& jet) const {
+      double totalTau(const fastjet::PseudoJet& /*jet*/) const {
          return _tau_components.tau();
       }
+      
       double subTau(const fastjet::PseudoJet& jet) const {
-         if (labelOf(jet) == -1) return NAN;
+         if (labelOf(jet) == -1) return std::numeric_limits<double>::quiet_NaN(); // nonsense
          return _tau_components.jet_pieces()[labelOf(jet)];
       }
@@ -93 +77 @@
          return (labelOf(jet) >= 0);
       }
-
+   
+private:
+   
+   TauComponents _tau_components;
+   std::vector<fastjet::PseudoJet> _jets;
+   std::vector<fastjet::PseudoJet> _axes;
+   
+   int labelOf(const fastjet::PseudoJet& jet) const {
+      int thisJet = -1;
+      for (unsigned int i = 0; i < _jets.size(); i++) {
+         if (_jets[i].cluster_hist_index() == jet.cluster_hist_index()) {
+            thisJet = i;
+            break;
+         }
+      }
+      return thisJet;
+   }
 };
 
@@ -122 +122 @@
  * onepass_wta_kt_axes  : one-pass minimization seeded by wta_kt
  *
- * For the unnormalized_measure, N-jettiness is defined as:
+ * For the UnnormalizedMeasure(beta), N-jettiness is defined as:
  *
  * tau_N = Sum_{all particles i} p_T^i min((DR_i1)^beta, (DR_i2)^beta, ...)
@@ -129 +129 @@
  *   and jet j.
  * 
- * The normalized_meausure include an extra parameter R0, and the various cutoff
+ * The NormalizedMeausure include an extra parameter R0, and the various cutoff
  * measures include an Rcutoff, which effectively defines an angular cutoff
  * similar in effect to a cone-jet radius.
@@ -138 +138 @@
 public:
 
-   NjettinessPlugin(int N,
-                    Njettiness::AxesMode axes_mode,
-                    Njettiness::MeasureMode measure_mode,
-                    double para1 = NAN,
-                    double para2 = NAN,
-                    double para3 = NAN,
-                    double para4 = NAN);
+   // Constructor with same arguments as Nsubjettiness.
+   NjettinessPlugin(int N,
+                    const AxesDefinition & axes_def,
+                    const MeasureDefinition & measure_def)
+   : _njettinessFinder(axes_def, measure_def), _N(N) {}
+   
+   
+   // Alternative constructors that define the measure via enums and parameters
+   // These constructors are likely be removed
+   NjettinessPlugin(int N,
+                 Njettiness::AxesMode axes_mode,
+                 Njettiness::MeasureMode measure_mode)
+   : _njettinessFinder(axes_mode, measure_mode, 0), _N(N) {}
+   
+   
+   NjettinessPlugin(int N,
+                 Njettiness::AxesMode axes_mode,
+                 Njettiness::MeasureMode measure_mode,
+                 double para1)
+   : _njettinessFinder(axes_mode, measure_mode, 1, para1), _N(N) {}
+   
+   
+   NjettinessPlugin(int N,
+                 Njettiness::AxesMode axes_mode,
+                 Njettiness::MeasureMode measure_mode,
+                 double para1,
+                 double para2)
+   : _njettinessFinder(axes_mode, measure_mode, 2, para1, para2), _N(N) {}
+   
+   
+   NjettinessPlugin(int N,
+                 Njettiness::AxesMode axes_mode,
+                 Njettiness::MeasureMode measure_mode,
+                 double para1,
+                 double para2,
+                 double para3)
+   : _njettinessFinder(axes_mode, measure_mode, 3, para1, para2, para3), _N(N) {}
+
 
    // Old constructor for backwards compatibility with v1.0,
-   // where normalized_cutoff_measure was the only option
+   // where NormalizedCutoffMeasure was the only option
    NjettinessPlugin(int N,
                     Njettiness::AxesMode mode,
                     double beta,
                     double R0,
-                    double Rcutoff=std::numeric_limits<double>::max());
+                    double Rcutoff=std::numeric_limits<double>::max())
+   : _njettinessFinder(mode, NormalizedCutoffMeasure(beta, R0, Rcutoff)), _N(N) {}
+
 
 
@@ -164 +197 @@
 private:
 
+   Njettiness _njettinessFinder;
    int _N;
-   mutable Njettiness _njettinessFinder; // TODO:  should muck with this so run_clustering can be const without this mutable
 
 };

external/fastjet/contribs/Nsubjettiness/Nsubjettiness.cc

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: Nsubjettiness.cc 597 2014-04-16 23:07:55Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.

external/fastjet/contribs/Nsubjettiness/Nsubjettiness.hh

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: Nsubjettiness.hh 670 2014-06-06 01:24:42Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -32 +33 @@
 #include <string>
 #include <climits>
-
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
@@ -49 +49 @@
 public:
 
-   // Main constructor, which takes N, axes/measure modes,
-   // and up to four parameters for parameters (i.e. beta, Rcutoff, etc depending on measure)
+   
+   // Main constructor, which takes N, the AxesDefiniation, and the MeasureDefinition.
+   // The Definitions are given in NjettinessDefinition.hh
+   //
+   // The recommended AxesDefinitions are (more are available as listed in the README
+   // and defined in NjettinessDefinition.hh):
+   //   KT_Axes             : exclusive kt axes
+   //   WTA_KT_Axes         : exclusive kt with winner-take-all recombination
+   //   OnePass_KT_Axes     : one-pass minimization from kt starting point
+   //   OnePass_WTA_KT_Axes : one-pass min. from wta_kt starting point
+   //
+   // The recommended measure definitions are (with the corresponding parameters)
+   //   NormalizedMeasure(beta,R0)
+   //      :  This was the original N-subjettiness measure (dimensionless)
+   //   UnnormalizedMeasure(beta)
+   //      :  This is the new recommended default, same as above but without
+   //      :  the normalization factor, and hence has units of GeV
+   //   NormalizedCutoffMeasure(beta,R0,Rcutoff)
+   //      :  Same as normalized_measure, but cuts off at Rcutoff
+   //   UnnormalizedCutoffMeasure(beta,Rcutoff)
+   //      :  Same as unnormalized_measure, but cuts off at Rcutoff
+   Nsubjettiness(int N,
+                 const AxesDefinition& axes_def,
+                 const MeasureDefinition& measure_def)
+   : _njettinessFinder(axes_def,measure_def), _N(N) {}
+   
+   
+   // Alternative constructors that define the measure via enums and parameters
+   // These constructors are likely be removed
+   // Zero parameter arguments
+   // (Currently, no measure uses this)
+   Nsubjettiness(int N,
+                 Njettiness::AxesMode axes_mode,
+                 Njettiness::MeasureMode measure_mode)
+   : _njettinessFinder(axes_mode, measure_mode, 0), _N(N) {}
+
+   // One parameter argument
+   // (for unnormalized_measure, para1=beta)
    Nsubjettiness(int N,
                  Njettiness::AxesMode axes_mode,
                  Njettiness::MeasureMode measure_mode,
-                 double para1 = NAN,
-                 double para2 = NAN,
-                 double para3 = NAN,
-                 double para4 = NAN)
-   : _njettinessFinder(axes_mode, measure_mode, para1, para2, para3, para4), _N(N) {}
+                 double para1)
+   : _njettinessFinder(axes_mode, measure_mode, 1, para1), _N(N) {}
+
+   // Two parameter arguments
+   // (for normalized_measure, para1=beta, para2=R0)
+   // (for unnormalized_cutoff_measure, para1=beta, para2=Rcutoff)
+   Nsubjettiness(int N,
+                 Njettiness::AxesMode axes_mode,
+                 Njettiness::MeasureMode measure_mode,
+                 double para1,
+                 double para2)
+   : _njettinessFinder(axes_mode, measure_mode, 2, para1, para2), _N(N) {}
+
+   // Three parameter arguments
+   // (for unnormalized_cutoff_measure, para1=beta, para2=R0, para3=Rcutoff)
+   Nsubjettiness(int N,
+                 Njettiness::AxesMode axes_mode,
+                 Njettiness::MeasureMode measure_mode,
+                 double para1,
+                 double para2,
+                 double para3)
+   : _njettinessFinder(axes_mode, measure_mode, 3, para1, para2, para3), _N(N) {}
 
    // Old constructor for backwards compatibility with v1.0,
@@ -67 +120 @@
                  double R0,
                  double Rcutoff=std::numeric_limits<double>::max())
-   : _njettinessFinder(axes_mode, Njettiness::normalized_cutoff_measure, beta, R0, Rcutoff), _N(N) {}
-
-
+   : _njettinessFinder(axes_mode, NormalizedCutoffMeasure(beta,R0,Rcutoff)), _N(N) {}
+   
    /// returns tau_N, measured on the constituents of this jet 
    double result(const PseudoJet& jet) const;
@@ -86 +138 @@
    }
    
+   /// returns subjet regions found by result() calculation (these have valid constituents)
+   /// Note that the axes and the subjets are not the same
+   std::vector<fastjet::PseudoJet> currentSubjets() const {
+      return _njettinessFinder.currentJets();
+   }
+
+   /// returns components of tau_N without recalculating anything
+   TauComponents currentTauComponents() const {
+      return _njettinessFinder.currentTauComponents();
+   }
+   
    // To set axes for manual use
-   void setAxes(std::vector<fastjet::PseudoJet> myAxes) {
+   void setAxes(const std::vector<fastjet::PseudoJet> & myAxes) {
       // Cross check that manual axes are being used is in Njettiness
         _njettinessFinder.setAxes(myAxes);
@@ -95 +158 @@
 private:
    
-   mutable Njettiness _njettinessFinder; // TODO:  should muck with this so result can be const without this mutable
+   Njettiness _njettinessFinder; // TODO:  should muck with this so result can be const without this mutable
    int _N;
 
@@ -112 +175 @@
    NsubjettinessRatio(int N,
                       int M,
+                      const AxesDefinition & axes_def,
+                      const MeasureDefinition & measure_def)
+   : _nsub_numerator(N,axes_def,measure_def),
+   _nsub_denominator(M,axes_def,measure_def) {}
+   
+   // Alternative constructor with enums and parameters
+   // Again, likely to be removed
+   NsubjettinessRatio(int N,
+                      int M,
+                      Njettiness::AxesMode axes_mode,
+                      Njettiness::MeasureMode measure_mode)
+   : _nsub_numerator(N, axes_mode, measure_mode),
+   _nsub_denominator(M, axes_mode, measure_mode) {}
+
+   
+   NsubjettinessRatio(int N,
+                      int M,
                       Njettiness::AxesMode axes_mode,
                       Njettiness::MeasureMode measure_mode,
-                      double para1 = NAN,
-                      double para2 = NAN,
-                      double para3 = NAN,
-                      double para4 = NAN)
-   : _nsub_numerator(N, axes_mode, measure_mode, para1, para2, para3, para4),
-   _nsub_denominator(M, axes_mode, measure_mode, para1, para2, para3, para4) {}
-
-   //returns tau_N/tau_M based off the input jet using result function from Nsubjettiness 
+                      double para1)
+   : _nsub_numerator(N, axes_mode, measure_mode, para1),
+   _nsub_denominator(M, axes_mode, measure_mode, para1) {}
+
+   NsubjettinessRatio(int N,
+                      int M,
+                      Njettiness::AxesMode axes_mode,
+                      Njettiness::MeasureMode measure_mode,
+                      double para1,
+                      double para2)
+   : _nsub_numerator(N, axes_mode, measure_mode, para1, para2),
+   _nsub_denominator(M, axes_mode, measure_mode, para1, para2) {}
+   
+   NsubjettinessRatio(int N,
+                      int M,
+                      Njettiness::AxesMode axes_mode,
+                      Njettiness::MeasureMode measure_mode,
+                      double para1,
+                      double para2,
+                      double para3)
+   : _nsub_numerator(N, axes_mode, measure_mode, para1, para2, para3),
+   _nsub_denominator(M, axes_mode, measure_mode, para1, para2, para3) {}
+
+   //returns tau_N/tau_M based off the input jet using result function from Nsubjettiness
    double result(const PseudoJet& jet) const;
 

external/fastjet/contribs/Nsubjettiness/VERSION

@@ -1 @@
-1.1.0-beta4
+2.1.0

external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.cc

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: WinnerTakeAllRecombiner.cc 597 2014-04-16 23:07:55Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -33 +34 @@
 
 // recombine pa and pb by creating pab with energy of the sum of particle energies in the direction of the harder particle
-
 // updated recombiner to use more general form of a metric equal to E*(pT/E)^(alpha), which reduces to pT*cosh(rap)^(1-alpha)
 // alpha is specified by the user. The default is alpha = 1, which is the typical behavior. alpha = 2 provides a metric which more 

external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh

@@ -5 +5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
+//  $Id: WinnerTakeAllRecombiner.hh 670 2014-06-06 01:24:42Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -50 +51 @@
    
    /// recombine pa and pb and put result into pab
-   virtual void recombine(const fastjet::PseudoJet & pa, const fastjet::PseudoJet & pb, fastjet::PseudoJet & pab) const;
+   virtual void recombine(const fastjet::PseudoJet & pa,
+                          const fastjet::PseudoJet & pb,
+                          fastjet::PseudoJet & pab) const;
 
 private: