Index: external/fastjet/contribs/Nsubjettiness/AUTHORS
===================================================================
--- external/fastjet/contribs/Nsubjettiness/AUTHORS	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/AUTHORS	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -18,3 +18,11 @@
    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.
+
 ----------------------------------------------------------------------
Index: external/fastjet/contribs/Nsubjettiness/AxesFinder.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/AxesFinder.cc	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/AxesFinder.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,5 +38,5 @@
 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,10 +46,5 @@
    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,5 +131,5 @@
 // (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,5 +162,5 @@
 // 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,5 +207,5 @@
 }
 
-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,5 +222,5 @@
    
 // 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,8 +252,8 @@
 // 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,9 +266,9 @@
          // 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,5 +283,5 @@
       // calculate tau on new axes
       seedAxes = newAxes;
-      double tempTau = _function->tau(particles, newAxes);
+      double tempTau = _function.tau(particles, newAxes);
       
       // close enough to stop?
Index: external/fastjet/contribs/Nsubjettiness/AxesFinder.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/AxesFinder.hh	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/AxesFinder.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,36 +54,19 @@
 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,15 +74,19 @@
 // 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,14 +96,17 @@
 // 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,48 +115,16 @@
 // 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,10 +133,11 @@
 // 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,10 +147,11 @@
 // 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,17 +162,19 @@
 // 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,6 +184,9 @@
 // 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,11 +197,12 @@
 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,32 +216,33 @@
 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,25 +256,23 @@
 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,28 +283,23 @@
 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,22 +309,5 @@
 // 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,4 +343,22 @@
    }
 
+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));
+   }
+   
+   
 };
 
Index: external/fastjet/contribs/Nsubjettiness/COPYING
===================================================================
--- external/fastjet/contribs/Nsubjettiness/COPYING	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/contribs/Nsubjettiness/COPYING	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,361 @@
+The Nsubjettiness contrib to FastJet is released 
+under the terms of the GNU General Public License v2 (GPLv2). 
+
+A copy of the GPLv2 is to be found at the end of this file.
+
+While the GPL license grants you considerable freedom, please bear in
+mind that this code's use falls under guidelines similar to those that
+are standard for Monte Carlo event generators
+(http://www.montecarlonet.org/GUIDELINES). In particular, if you use
+this code as part of work towards a scientific publication, whether
+directly or contained within another program you should include a citation 
+to
+
+  arXiv:1011.2268
+  arXiv:1108.2701
+
+
+
+======================================================================
+======================================================================
+======================================================================
+		    GNU GENERAL PUBLIC LICENSE
+		       Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+			    Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users.  This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it.  (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.)  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+  To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have.  You must make sure that they, too, receive or can get the
+source code.  And you must show them these terms so they know their
+rights.
+
+  We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+  Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software.  If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+  Finally, any free program is threatened constantly by software
+patents.  We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary.  To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+		    GNU GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License.  The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language.  (Hereinafter, translation is included without limitation in
+the term "modification".)  Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+  1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+  2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) You must cause the modified files to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    b) You must cause any work that you distribute or publish, that in
+    whole or in part contains or is derived from the Program or any
+    part thereof, to be licensed as a whole at no charge to all third
+    parties under the terms of this License.
+
+    c) If the modified program normally reads commands interactively
+    when run, you must cause it, when started running for such
+    interactive use in the most ordinary way, to print or display an
+    announcement including an appropriate copyright notice and a
+    notice that there is no warranty (or else, saying that you provide
+    a warranty) and that users may redistribute the program under
+    these conditions, and telling the user how to view a copy of this
+    License.  (Exception: if the Program itself is interactive but
+    does not normally print such an announcement, your work based on
+    the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+    a) Accompany it with the complete corresponding machine-readable
+    source code, which must be distributed under the terms of Sections
+    1 and 2 above on a medium customarily used for software interchange; or,
+
+    b) Accompany it with a written offer, valid for at least three
+    years, to give any third party, for a charge no more than your
+    cost of physically performing source distribution, a complete
+    machine-readable copy of the corresponding source code, to be
+    distributed under the terms of Sections 1 and 2 above on a medium
+    customarily used for software interchange; or,
+
+    c) Accompany it with the information you received as to the offer
+    to distribute corresponding source code.  (This alternative is
+    allowed only for noncommercial distribution and only if you
+    received the program in object code or executable form with such
+    an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it.  For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable.  However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License.  Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+  5. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Program or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+  6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+  7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded.  In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+  9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation.  If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+  10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission.  For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this.  Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+			    NO WARRANTY
+
+  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+		     END OF TERMS AND CONDITIONS
+
+	    How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+    Gnomovision version 69, Copyright (C) year  name of author
+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+  <signature of Ty Coon>, 1 April 1989
+  Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs.  If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library.  If this is what you want to do, use the GNU Library General
+Public License instead of this License.
Index: external/fastjet/contribs/Nsubjettiness/ChangeLog
===================================================================
--- external/fastjet/contribs/Nsubjettiness/ChangeLog	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/contribs/Nsubjettiness/ChangeLog	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,181 @@
+2014-07-09 <JDT>
+   Changed version for 2.1.0 release.
+   Updated NEWS to reflect 2.1.0 release
+2014-07-07 <JDT>
+   Added forward declaration of options in NjettinessDefinition for readability.
+   Updated README with some clarifications
+   Added usage information in the example file
+   Reran svn propset svn:keywords Id *.cc *.hh
+2014-06-25 <JDT>
+   Declaring release candidate of 2.1
+2014-06-11 <JDT>
+   Fixed virtual destructor issue in AxesFinder
+   Changing copy() to create() in NjettinessDefinition for "new" clarity
+   Converted some SharedPtr to regular pointers in NjettinessDefinition to be consistent on meaning of "create" commands.
+2014-06-10 <JDT>
+   Slight modification of example_advanced_usage
+   Fixed bug in GeometricCutoffMeasure (incorrect denominator setting)
+2014-06-05 <JDT>
+   Moved public before private in the .hh files for readability
+   Starting process of switching to SharedPtr internally
+   Clean up of AxesFinderFromGeometricMinimization
+   Simplified AxesFinder interface such that it doesn't know about starting axes finders (this is now handled in Njettiness).
+   Added const qualifiers in Njettiness
+2014-06-04 <JDT>
+   Implemented AxesDefinition class
+   Added descriptions to AxesDefinition and MeasureDefinition
+   Simplified example_advanced_usage with new Definitions
+   Made copy constructor private for Njettiness, to avoid copying 
+2014-06-03 <JDT>
+   Implemented remaining suggestions from FJ authors (Thanks!)
+   Fixed bug in example_advanced_usage where wrong beta value was used for NjettinessPlugin tests.
+   Removed NANs as signals for number of parameters in Nsubjettiness and NjettinessPlugin
+   Reduced the number of allowed parameters from 4 to 3.
+   Wrapped NEWS to 80 characters
+   Added MeasureDefinition as way to safely store information about the measures used
+   Converted a few NANs to std::numeric_limits<double>::quiet_NaN() when a parameter shouldn't be used.
+   Added AxesStruct and MeasureStruct to simplify the form of example_advanced_usage
+   Added example_v1p0p3 to check for backwards compatibility with v1.0.3
+   Changed the names of the MeasureFunctions in order to avoid conflicts with the new MeasureDefinitions
+   Fixed bug in correlation between subjets and tau values in NjettinessPlugin
+   Added currentTauComponents to Nsubjettiness
+   Added subTau information to example_basic_usage
+   Added file NjettinessDefinition to hold MeasureDefinition
+   Changed Njettiness constructors to treat MeasureSpecification as primary object
+   Fixed segmentation fault with ClusterSequenceAreas
+2014-06-02 <JDT>
+   Implemented many suggestions from FJ authors (Thanks!)
+   Removed FastJet 2 specific code
+   Made sq() function into internal namespace (as "inline static" to avoid conflicts with other packages)
+   Made setAxes() take const reference argument
+   Rewrapped README to 80 characters and updated/improved some of the descriptions
+   Clarified NEWS about what parts of the Nsubjettiness code is backwards compatible with v1.0.3
+   Clarified the para choices in Nsubjettiness constructor
+2014-04-30 <JDT>
+   Added (void)(n_jets) in AxesFinder.hh to fix unused-parameter warning
+2014-04-29 <JDT>
+   Added manual definition of NAN for compilers that don't have it.
+   Removed a few more unused parameters for compilation
+2014-04-22 <JDT>
+   Turned on -Wunused-parameter compiler flag to fix ATLAS compile issues.
+2014-04-18 <JDT>
+   Tweaks to NEWS and README.  Preparing for 2.0.0-rc1 release.
+2014-04-16 <JDT>
+   Decided that enough has changed that this should be v2.0
+   Added Id tags
+2014-04-14 <JDT>
+   Added get_partition_list to MeasureFunction
+   Removed do_cluster from MeasureFunction (no longer needed)
+   Fixed bug with NjettinessPlugin where jets were listed in backwards order from axes.
+   Removed various commented out pieces of code.
+2014-03-16 <JDT>
+   Added partitioning information to Nsubjettiness
+   Partitioning is now calculated in MeasureFunction and stored by Njettiness.
+   Rewrote MeasureFunction result() to call result_from_partition()
+   Added subjet (and constituent counting) information to example_basic_usage
+   Commented out redundant "getJets" function
+2014-02-25 <JDT>
+   Added access to seedAxes used for one-pass minimization routines.
+   Added axes print out to example_basic_usage, and fixed too many PrintJets declarations
+2014-02-24 <JDT>
+   Fixed embarrassing bug with min_axes (error introduced after v1.0 to make it the same as onepass_kt)
+   Simplified GeometricMeasure and added possibility of beta dependence
+   Commented out WTA2 options, since those have not been fully tested (nor do they seem particularly useful at the moment).  They can be reinstated if the physics case can be made to use them.
+   Split example into example_basic_usage and example_advanced_usage
+2014-01-28 <TJ>
+   Added new options in WinnerTakeAllRecombiner to use either pT or pT^2/E to recombine particles
+2014-01-24 <JDT>
+   Added access to currentAxes from Nsubjettiness.
+2014-01-18 <JDT>
+   Added beam regions to MeasureFunction, correspondingly renamed functions to have jet and beam regions
+   Renamed functions in TauComponents for consistency with MeasureFunction
+   Adding debugging code to AxesFinderFromOnePassMinimization::getAxes
+   Worked extensively on example.cc to make sure that it tests all available options.
+   Rewrote PrintJets command in example.cc for later improvements
+   Converted some magic numbers to std::numeric_limits<double>::max()
+2014-01-17 <JDT>
+   Rewrote KMeansMinimization to call OnePassMinimization, adding noise explicitly.
+   Removed any nothing of noise from OnePassMinimization
+   Removed Double32_t for root usage is Nsubjettiness
+   Clean up of many comments throughout the code, updating of README file
+   Removed unnecessary establishAxes in Njettiness
+   Removed bare constructor for Njettiness to avoid incompatibility with enum choices, may reinstate later.  Also removed setMeasureFunction, setAxesFinder for same reason
+   NjettinessExtras now calls TauComponents
+2014-01-16 <TJ>
+   Moved minimization functions to OnePassMinimization, changed KMeansMinimization class to simply call OnePassMinimization a specified number of times
+   Added extra tau function in TauComponents for users to get tau directly
+   Changed radius parameter in AxesFinderFromExclusiveJet subclasses to use max_allowable_R
+   Updated example.ref to account for changes due to change in radius parameter
+2014-01-15 <TJ>
+   Changed NjettinessComponents to TauComponents
+   Updated MeasureFunction with "result" function that returns TauComponents object
+   TauComponents changed to calculate all tau components given subtaus_numerator and tau_denominator
+   Njettiness updated to return TauComponents object rather than individual components
+   Nsubjettiness and NjettinessPlugin updated to have option for 4th parameter
+2014-01-14 <TJ>
+   Added NjettinessComponents class so Njettiness does not recalculate tau values
+   Removed old Njettiness constructors, updated Nsubjettiness and NjettinessPlugin constructors to use new constructor
+   Added geometric minimization to OnePassAxesFinders
+   Created new Njettiness function to set OnePassAxesFinders to reduce code
+   Updated LightLikeAxis with ConvertToPseudoJet function
+   Updated README with new functionality of code
+2014-01-12 <TJ>
+   Removed NsubGeometricParameters in all functions/constructors, replaced with Rcutoff double
+   Added three new measure mode options where Rcutoff is declared explicitly in parameters
+   Added checks so minimization axes finders are not used for geometric measures
+   AxesFinderFromOnePassMinimization class created as child of AxesFinderFromKmeansMinimization
+   Added new NsubjettinessRatio constructor to include MeasureMode option
+   Moved AxesFinder and MeasureFunction declarations from AxesMode and MeasureMode into separate Njettiness function
+   Removed R0 from AxesFinderFromKmeansMinimization
+   Changed example.cc to get rid of use of NsubGeometricParameters
+2014-01-9 <TJ>
+   Removed NsubParameters in all functions/constructors, replaced with three separate parameters
+   Added checks for correct number of parameters in Njettiness constructor
+2014-01-8 <TJ>
+   Removed normalization information from Nsubjettiness
+   Added flag to MeasureFunction to give option of using the denominator
+   Split DefaultMeasure into separate normalized and unnormalized classes
+2014-01-7 <TJ>
+   Added capability of choosing a specific Measure in Njettiness
+   Added new Nsubjettiness constructor to allow choice of both AxesMode and MeasureMode
+2014-01-6 <TJ>
+   Updated copyright information
+   Fixed bug in WinnerTakeAllRecombiner
+   Moved KMeansParameters to AxesFinder
+   Updated README with descriptions of new header files
+2013-12-30 <TJ>
+   Changed name of MeasureFunctor to MeasureFunction
+   Created separate .hh/.cc files for MeasureFunction, AxesFinder, and WinnerTakeAllRecombiner
+   Updated Makefile to account for new files
+   Removed getMinimumAxes in AxesFinderFromKMeansMinimization, consolidated with getAxes
+   Updated comments on classes and major functions
+2013-12-22 <TJ>
+   Created .cc files and moved all function definitions into .cc files
+   Updated Makefile to account for new .cc files
+2013-11-12 <TJ>
+   Added to fjcontrib svn
+2013-11-12 <jthaler>
+   Debugging svn
+2013-11-11 <TJ>
+   Changed MeasureFunctor to separately treat tau numerator and denominator
+   Changed some of the function names in MeasureFunctor.  Should not affect users
+   Added more informative function names to Njettiness.
+   Njettiness now allows finding unnormalized tau values
+   Added WTARecombiner to define winner-take-all axes
+   Added various WTA options to AxesMode
+   Added setAxes to Nsubjettiness
+   Added NsubjettinessRatio function
+2013-08-26 <jthaler>
+   Added inlines to fix compile issue
+   Put some of the minimization code inside of the AxesFinderFromKmeansMinimization class
+2013-02-23 <jthaler>
+   Fixed dependency issue (now using make depend)
+2013-02-22 <jthaler>
+   Fixed memory management and failed make check issues.
+2013-02-21 <jthaler>
+   First version submitted to fjcontrib
+2013-02-20 <jthaler>
+   Initial creation based on previous plugin hosted at http://www.jthaler.net/jets/
+
+
+
Index: external/fastjet/contribs/Nsubjettiness/Makefile
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Makefile	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/contribs/Nsubjettiness/Makefile	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,101 @@
+# If you are using this Makefile standalone and fastjet-config is not
+# in your path, edit this line to specify the full path
+FASTJETCONFIG=fastjet-config
+PREFIX=`$(FASTJETCONFIG) --prefix`
+CXX=g++
+CXXFLAGS= -O3 -Wall -Woverloaded-virtual -g -Wunused-parameter
+install_script = $(SHELL) ../utils/install-sh
+check_script = ../utils/check.sh
+
+# global contrib-wide Makefile include may override some of the above
+# variables (leading "-" means don't give an error if you can't find
+# the file)
+-include ../.Makefile.inc
+
+#------------------------------------------------------------------------
+# things that are specific to this contrib
+NAME=Nsubjettiness
+SRCS=Nsubjettiness.cc Njettiness.cc NjettinessPlugin.cc MeasureFunction.cc AxesFinder.cc WinnerTakeAllRecombiner.cc NjettinessDefinition.cc
+EXAMPLES=example_basic_usage example_advanced_usage example_v1p0p3
+INSTALLED_HEADERS=Nsubjettiness.hh Njettiness.hh NjettinessPlugin.hh MeasureFunction.hh AxesFinder.hh WinnerTakeAllRecombiner.hh NjettinessDefinition.hh
+#------------------------------------------------------------------------
+
+CXXFLAGS+= $(shell $(FASTJETCONFIG) --cxxflags)
+LDFLAGS += -lm $(shell $(FASTJETCONFIG) --libs)
+
+OBJS  = $(SRCS:.cc=.o)
+EXAMPLES_SRCS  = $(EXAMPLES:=.cc)
+
+install_HEADER  = $(install_script) -c -m 644
+install_LIB     = $(install_script) -c -m 644
+install_DIR     = $(install_script) -d
+install_DATA    = $(install_script) -c -m 644
+install_PROGRAM = $(install_script) -c -s
+install_SCRIPT  = $(install_script) -c
+
+.PHONY: clean distclean examples check install
+
+# compilation of the code (default target)
+all: lib$(NAME).a
+
+lib$(NAME).a: $(OBJS)
+	ar cru lib$(NAME).a $(OBJS)
+	ranlib lib$(NAME).a
+
+# building the examples
+examples: $(EXAMPLES)
+
+# the following construct alloews to build each of the examples listed
+# in $EXAMPLES automatically
+$(EXAMPLES): % : %.o all
+	$(CXX) -o $@ $< -L. -l$(NAME) $(LDFLAGS)
+
+# check that everything went fine
+check: examples
+	@for prog in $(EXAMPLES); do\
+	  $(check_script) $${prog} ../data/single-event.dat || exit 1; \
+	done
+	@echo "All tests successful"
+
+# cleaning the directory
+clean:
+	rm -f *~ *.o *.a
+
+distclean: clean
+	rm -f lib$(NAME).a $(EXAMPLES)
+
+# install things in PREFIX/...
+install: all
+	$(install_DIR) $(PREFIX)/include/fastjet/contrib
+	for header in $(INSTALLED_HEADERS); do\
+	  $(install_HEADER) $$header $(PREFIX)/include/fastjet/contrib/;\
+	done
+	$(install_DIR) $(PREFIX)/lib
+	$(install_LIB) lib$(NAME).a $(PREFIX)/lib
+
+depend:
+	makedepend -Y --   -- $(SRCS) $(EXAMPLES_SRCS)
+# DO NOT DELETE
+
+Nsubjettiness.o: Nsubjettiness.hh Njettiness.hh MeasureFunction.hh
+Nsubjettiness.o: AxesFinder.hh WinnerTakeAllRecombiner.hh
+Nsubjettiness.o: NjettinessDefinition.hh
+Njettiness.o: Njettiness.hh MeasureFunction.hh AxesFinder.hh
+Njettiness.o: WinnerTakeAllRecombiner.hh NjettinessDefinition.hh
+NjettinessPlugin.o: NjettinessPlugin.hh Njettiness.hh MeasureFunction.hh
+NjettinessPlugin.o: AxesFinder.hh WinnerTakeAllRecombiner.hh
+NjettinessPlugin.o: NjettinessDefinition.hh
+MeasureFunction.o: MeasureFunction.hh
+AxesFinder.o: AxesFinder.hh WinnerTakeAllRecombiner.hh MeasureFunction.hh
+WinnerTakeAllRecombiner.o: WinnerTakeAllRecombiner.hh
+NjettinessDefinition.o: NjettinessDefinition.hh MeasureFunction.hh
+NjettinessDefinition.o: AxesFinder.hh WinnerTakeAllRecombiner.hh
+example_basic_usage.o: Nsubjettiness.hh Njettiness.hh MeasureFunction.hh
+example_basic_usage.o: AxesFinder.hh WinnerTakeAllRecombiner.hh
+example_basic_usage.o: NjettinessDefinition.hh NjettinessPlugin.hh
+example_advanced_usage.o: Nsubjettiness.hh Njettiness.hh MeasureFunction.hh
+example_advanced_usage.o: AxesFinder.hh WinnerTakeAllRecombiner.hh
+example_advanced_usage.o: NjettinessDefinition.hh NjettinessPlugin.hh
+example_v1p0p3.o: Nsubjettiness.hh Njettiness.hh MeasureFunction.hh
+example_v1p0p3.o: AxesFinder.hh WinnerTakeAllRecombiner.hh
+example_v1p0p3.o: NjettinessDefinition.hh NjettinessPlugin.hh
Index: external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,10 +37,24 @@
 
 // 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,23 +75,98 @@
       
       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
Index: external/fastjet/contribs/Nsubjettiness/MeasureFunction.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/MeasureFunction.hh	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/MeasureFunction.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,37 +32,21 @@
 #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,6 +100,28 @@
    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,4 +132,47 @@
 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,141 +183,113 @@
    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
 
Index: external/fastjet/contribs/Nsubjettiness/NEWS
===================================================================
--- external/fastjet/contribs/Nsubjettiness/NEWS	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/contribs/Nsubjettiness/NEWS	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,83 @@
+-------------------------
+Version 2
+-------------------------
+
+This is a streamlining of the N-subjettiness code, developed mainly by TJ
+Wilkason.  The core functionality is unchanged, but classes have been
+dramatically reorganized to allow for later expansion.  Because the API for
+Njettiness has changed, we have called this v2 (http://semver.org).
+
+Note that we have maintain backwards compatibility for the typical ways that
+Nsubjettiness was used.  In particular, all of the Nsubjettiness class code in
+the example file v1.0.3 from still compiles, as does the NjettinessPlugin class
+code that uses the default measure.
+
+The key new features are:
+
+  * NsubjettinessRatio:  Direct access to tau_N / tau_M (the most requested
+    feature)
+  * MeasureDefinition to allow access to both normalized and unnormalized measures
+  * AxesDefinition to allow for access to more general axes modes
+  * Winner-Take-All recombination axes:  a faster way to find axes than beta = 1
+    minimization, but with comparable performance.
+  * TauComponents to get access to the pieces of the N-(sub)jettiness
+    calculation.
+  * For clarity, split the example file into an example_basic_usage and
+    example_advanced_usage
+  * In Nsubjettiness, access to seedAxes() and currentAxes() to figure out the
+    axes used before and after minimization.
+  * In Nsubjettiness, access to currentSubjets() to get the subjet fourvectors.
+
+
+-- 2.1.0:  (July 9, 2014) Inclusion of MeasureDefinition/AxesDefinition interface.
+           This was the first publicly available version of Nsubjettiness v2.
+-- 2.0.0:  Initial release of v2.0.  This was never officially made public.
+
+-------------------------
+Version 1
+-------------------------
+
+This was a new release using FastJet contrib framework, primary developed by
+Jesse Thaler.
+
+-- 1.0.3:  Added inlines to fix compile issue (thanks Matthew Low)
+-- 1.0.2:  Fixed potential dependency issue (thanks FJ authors)
+-- 1.0.1:  Fixed memory leak in Njettiness.hh (thanks FJ authors)
+-- 1.0.0:  New release using FastJet contrib framework.  This includes a new
+makefile and a simplified example program.
+
+-------------------------
+Previous Versions
+-------------------------
+
+The previous versions of this code were developed initially by Ken Van Tilburg,
+tweaked by Jesse Thaler, and made into a robust FastJet add on by Chris
+Vermilion.
+
+Previous versions available from:
+   http://jthaler.net/jets/Njettiness-0.5.1.tar.gz (Experimental Version)
+   http://jthaler.net/jets/Njettiness-0.4.1.tar.gz (Stable Version)
+
+Previous version history:
+-- 0.5.1:  For Njettiness Plugin, added access to currentTau values and axes via
+           ClusterSequence::Extras class.  (thanks to Dinko Ferencek and John
+           Paul Chou)
+-- 0.5.0:  Corrected fatal error in ConstituentTauValue (TauValue unaffected).
+           Started process of allowing for more general measures and alternative
+           minimization schemes. Extremely preliminary inclusion of alternative
+           "geometric" measure.
+-- 0.4.1:  Corrected bug where a too-small value of Rcut would cause the
+           minimization procedure to fail (thanks Marat Freytsis, Brian Shuve)
+-- 0.4.0:  Adding Nsubjettiness FunctionOfPseudoJet<float>.  Re-organizing file
+           structure and doing some re-naming to clarify Njettiness vs.
+           Nsubjettiness.  Some speedup in UpdateAxes code. (CKV)
+-- 0.3.2:  Returns zero instead of a segmentation fault when the number of
+           particles in a jet is smaller than the N value in tau_N (thanks
+           Grigory Ovanesyan)
+-- 0.3.2:  Fixed -Wshadow errors (thanks Grigory Ovanesyan)
+-- 0.3.1:  Fixed stray comma/semicolon compiler error (thanks Grigory Ovanesyan)
+-- 0.3.1:  Corrected tarbomb issue (thanks Jonathan Walsh)
+-- 0.3.1:  Added anti-kT seeds as option
+-- 0.3.1:  Fixed bug in minimization code with R_cutoff (thanks Chris Vermilion)
+-- 0.3.1:  Added getPartition() and getJets() functions as helper functions for
+           Chris Vermilion. (JT)
Index: external/fastjet/contribs/Nsubjettiness/Njettiness.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Njettiness.cc	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/Njettiness.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,28 +36,25 @@
 ///////
 
-// 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,25 +62,25 @@
          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,109 +88,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,5 +177,5 @@
 // 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,8 +183,21 @@
       _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,42 +209,10 @@
 // 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
 
Index: external/fastjet/contribs/Nsubjettiness/Njettiness.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Njettiness.hh	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/Njettiness.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,17 +26,21 @@
 #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,4 +59,5 @@
    
    // 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,4 +84,5 @@
    // "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,93 +94,86 @@
    };
 
-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
 
Index: external/fastjet/contribs/Nsubjettiness/NjettinessDefinition.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/NjettinessDefinition.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/contribs/Nsubjettiness/NjettinessDefinition.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,76 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//  $Id: NjettinessDefinition.cc 704 2014-07-07 14:30:43Z jthaler $
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+
+#include "NjettinessDefinition.hh"
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib {
+
+   std::string NormalizedMeasure::description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Normalized Measure (beta = " << _beta << ", R0 = " << _R0 << ")";
+      return stream.str();
+   };
+   
+   std::string UnnormalizedMeasure::description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Unnormalized Measure (beta = " << _beta << ", in GeV)";
+      return stream.str();
+   };
+   
+   std::string GeometricMeasure::description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Geometric Measure (beta = " << _beta << ", in GeV)";
+      return stream.str();
+   };
+   
+   std::string NormalizedCutoffMeasure::description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Normalized Cutoff Measure (beta = " << _beta << ", R0 = " << _R0 << ", Rcut = " << _Rcutoff << ")";
+      return stream.str();
+   };
+   
+   std::string UnnormalizedCutoffMeasure::description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Unnormalized Cutoff Measure (beta = " << _beta << ", Rcut = " << _Rcutoff << ", in GeV)";
+      return stream.str();
+   };
+   
+   std::string GeometricCutoffMeasure::description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Geometric Cutoff Measure (beta = " << _beta << ", Rcut = " << _Rcutoff << ", in GeV)";
+      return stream.str();
+   };
+   
+   
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/Nsubjettiness/NjettinessDefinition.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/NjettinessDefinition.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/contribs/Nsubjettiness/NjettinessDefinition.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,721 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//  $Id: NjettinessDefinition.hh 704 2014-07-07 14:30:43Z jthaler $
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+
+#ifndef __FASTJET_CONTRIB_NJETTINESS_DEFINITION_HH__
+#define __FASTJET_CONTRIB_NJETTINESS_DEFINITION_HH__
+
+
+#include "MeasureFunction.hh"
+#include "AxesFinder.hh"
+
+#include "fastjet/PseudoJet.hh"
+#include <fastjet/LimitedWarning.hh>
+
+#include <iomanip>
+#include <cmath>
+#include <vector>
+#include <list>
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib {
+   
+// Eventually this file might contains an NjettinessDefinition that combines an
+// AxesDefinition with a MeasureDefinition. It's not clear how useful that would be, though
+  
+// The following Measures are available (and the relevant arguments):
+class NormalizedMeasure;         // (beta,R0)
+class UnnormalizedMeasure;       // (beta)
+class GeometricMeasure;          // (beta)
+class NormalizedCutoffMeasure;   // (beta,R0,Rcutoff)
+class UnnormalizedCutoffMeasure; // (beta,Rcutoff)
+class GeometricCutoffMeasure;    // (beta,Rcutoff)
+
+// The following Axes are available (and the relevant arguments, if needed)
+class KT_Axes;
+class CA_Axes;
+class AntiKT_Axes;   // (R0)
+class WTA_KT_Axes;
+class WTA_CA_Axes;
+class Manual_Axes;
+class OnePass_KT_Axes;
+class OnePass_CA_Axes;
+class OnePass_AntiKT_Axes;   // (R0)
+class OnePass_WTA_KT_Axes;
+class OnePass_WTA_CA_Axes;
+class OnePass_Manual_Axes;
+class MultiPass_Axes;
+
+// Below are just technical implementations of the variable axes and measures.
+  
+///////
+//
+// MeasureDefinition
+//
+///////
+
+//The MeasureDefinition is a way to set the MeasureMode and
+//parameters used by Njettiness, with errors in the number of parameters given at
+//compile time (instead of at run time).  The MeasureDefintion knows which core objects
+//to call to make the measurement, as well as properties of the MeasureFunction
+class MeasureDefinition {
+   
+public:
+
+   // Description of measure and parameters
+   virtual std::string description() const = 0;
+   
+   // In derived classes, this should return a copy of the corresponding
+   // derived class
+   virtual MeasureDefinition* create() const = 0;
+   
+   //Return the MeasureFunction corresponding to this definition
+   virtual MeasureFunction* createMeasureFunction() const = 0;
+   
+   //Return the AxesFinder that should be used for one-pass minimization
+   virtual AxesFinder* createOnePassAxesFinder() const = 0;
+
+   //Specify whether multi-pass minimization makes sense, and if so, return the AxesFinder that should be used for multi-pass minimization
+   virtual bool supportsMultiPassMinimization() const = 0;
+   virtual AxesFinder* createMultiPassAxesFinder( unsigned int) const = 0;
+   
+   virtual ~MeasureDefinition(){};
+};
+
+// The normalized measure, with two parameters: beta and R0
+class NormalizedMeasure : public MeasureDefinition {
+   
+public:
+   NormalizedMeasure(double beta, double R0)
+   : _beta(beta), _R0(R0) {}
+   
+   virtual std::string description() const;
+
+   virtual MeasureFunction* createMeasureFunction() const {
+      return (new DefaultNormalizedMeasureFunction(_beta,_R0,std::numeric_limits<double>::max()));
+   }
+
+   virtual NormalizedMeasure* create() const {return new NormalizedMeasure(*this);}
+   
+   virtual AxesFinder* createOnePassAxesFinder() const {
+      return (new AxesFinderFromOnePassMinimization(_beta, std::numeric_limits<double>::max()));
+   }
+   virtual AxesFinder* createMultiPassAxesFinder( unsigned int Npass) const {
+      return (new AxesFinderFromKmeansMinimization(_beta, std::numeric_limits<double>::max(),Npass));
+   }
+
+   virtual bool supportsMultiPassMinimization() const { return true; }
+   
+private:
+   double _beta;
+   double _R0;
+};
+   
+// The unnormalized measure, with just one parameter: beta
+class UnnormalizedMeasure : public MeasureDefinition {
+
+public:
+   UnnormalizedMeasure(double beta)
+   : _beta(beta) {}
+
+   virtual UnnormalizedMeasure* create() const {return new UnnormalizedMeasure(*this);}
+   
+   virtual std::string description() const;
+   
+   virtual MeasureFunction* createMeasureFunction() const {
+      return (new DefaultUnnormalizedMeasureFunction(_beta,std::numeric_limits<double>::max()));
+   }
+
+   
+   virtual AxesFinder* createOnePassAxesFinder() const {
+      return (new AxesFinderFromOnePassMinimization(_beta, std::numeric_limits<double>::max()));
+   }
+   virtual AxesFinder* createMultiPassAxesFinder( unsigned int Npass) const {
+      return (new AxesFinderFromKmeansMinimization(_beta, std::numeric_limits<double>::max(),Npass));
+   }
+
+   virtual bool supportsMultiPassMinimization() const { return true; }
+
+private:
+   double _beta;
+};
+   
+   
+// The geometric  measure, with 1 parameter: beta
+// This measure is still evolving and shouldn't be used for any critial applications yet
+class GeometricMeasure : public MeasureDefinition {
+   
+public:
+   GeometricMeasure(double beta)
+   : _beta(beta) {}
+   
+   virtual GeometricMeasure* create() const {return new GeometricMeasure(*this);}
+   
+   virtual std::string description() const;
+   
+   virtual MeasureFunction* createMeasureFunction() const {
+      return (new GeometricMeasureFunction(_beta,std::numeric_limits<double>::max()));
+   }
+   
+   virtual AxesFinder* createOnePassAxesFinder() const {
+      return (new AxesFinderFromGeometricMinimization(_beta, std::numeric_limits<double>::max()));
+   }
+   virtual AxesFinder* createMultiPassAxesFinder(unsigned int) const {
+      throw Error("GeometricMeasure does not support multi-pass minimization.");
+      return NULL;
+   }
+   
+   virtual bool supportsMultiPassMinimization() const { return false; }
+
+private:
+   double _beta;
+
+};
+
+
+// The normalized cutoff measure, with 3 parameters: beta, R0, Rcutoff
+class NormalizedCutoffMeasure : public MeasureDefinition {
+
+public:
+   NormalizedCutoffMeasure(double beta, double R0, double Rcutoff)
+   : _beta(beta), _R0(R0), _Rcutoff(Rcutoff) {}
+
+   virtual std::string description() const;
+   
+   virtual NormalizedCutoffMeasure* create() const {return new NormalizedCutoffMeasure(*this);}
+   
+   virtual MeasureFunction* createMeasureFunction() const {
+      return (new DefaultNormalizedMeasureFunction(_beta,_R0,_Rcutoff));
+   }
+   
+   virtual AxesFinder* createOnePassAxesFinder() const {
+      return (new AxesFinderFromOnePassMinimization(_beta, _Rcutoff));
+   }
+   virtual AxesFinder* createMultiPassAxesFinder( unsigned int Npass) const {
+      return (new AxesFinderFromKmeansMinimization(_beta, std::numeric_limits<double>::max(),Npass));
+   }
+  
+   virtual bool supportsMultiPassMinimization() const { return true; }
+
+private:
+   double _beta;
+   double _R0;
+   double _Rcutoff;
+   
+};
+
+// The unnormalized cutoff measure, with 2 parameters: beta, Rcutoff
+class UnnormalizedCutoffMeasure : public MeasureDefinition {
+   
+public:
+   UnnormalizedCutoffMeasure(double beta, double Rcutoff)
+   : _beta(beta), _Rcutoff(Rcutoff) {}
+
+   virtual std::string description() const;
+   
+   virtual UnnormalizedCutoffMeasure* create() const {return new UnnormalizedCutoffMeasure(*this);}
+
+   virtual MeasureFunction* createMeasureFunction() const {
+      return (new DefaultUnnormalizedMeasureFunction(_beta,_Rcutoff));
+   }
+   
+   virtual AxesFinder* createOnePassAxesFinder() const {
+      return (new AxesFinderFromOnePassMinimization(_beta, _Rcutoff));
+   }
+   virtual AxesFinder* createMultiPassAxesFinder( unsigned int Npass) const {
+      return (new AxesFinderFromKmeansMinimization(_beta, std::numeric_limits<double>::max(),Npass));
+   }
+
+   virtual bool supportsMultiPassMinimization() const { return true; }
+
+private:
+   double _beta;
+   double _Rcutoff;
+   
+};
+
+// The Geometric  measure, with 2 parameters: beta, Rcutoff
+// This measure is still evolving and shouldn't be used for any critial applications yet
+class GeometricCutoffMeasure : public MeasureDefinition {
+   
+public:
+   GeometricCutoffMeasure(double beta, double Rcutoff)
+   : _beta(beta), _Rcutoff(Rcutoff) {}
+
+   virtual GeometricCutoffMeasure* create() const {return new GeometricCutoffMeasure(*this);}
+
+   virtual std::string description() const;
+   
+   virtual MeasureFunction* createMeasureFunction() const {
+      return (new GeometricMeasureFunction(_beta,_Rcutoff));
+   }
+   
+   virtual AxesFinder* createOnePassAxesFinder() const {
+      return (new AxesFinderFromGeometricMinimization(_beta, _Rcutoff));
+   }
+   virtual AxesFinder* createMultiPassAxesFinder(unsigned int) const {
+      throw Error("GeometricCutoffMeasure does not support multi-pass minimization.");
+      return NULL;
+   }
+   
+   virtual bool supportsMultiPassMinimization() const { return false; }
+
+private:
+   double _beta;
+   double _Rcutoff;
+   
+};
+   
+   
+///////
+//
+// AxesDefinition
+//
+///////
+
+//Analogous to MeasureDefinition, AxesDefinition defines which AxesFinder to use
+//At the moment, most AxesDefinition do not have an arugment (except the Anti-KT ones)
+class AxesDefinition {
+   
+public:
+   // description of axes (and any parameters)
+   virtual std::string short_description() const = 0;
+   virtual std::string description() const = 0;
+   
+   // In derived classes, this should return a copy of the corresponding
+   // derived class
+   virtual AxesDefinition* create() const = 0;
+   
+   // These describe how the axes finder works
+   virtual bool givesRandomizedResults() const = 0;
+   virtual bool supportsManualAxes() const = 0;
+   
+   //for best testing
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition &) const = 0;
+   virtual AxesFinder* createFinishingAxesFinder(const MeasureDefinition &) const {
+      return NULL;  //By default, nothing.
+   };
+
+   virtual ~AxesDefinition() {};
+};
+
+// kt axes
+class KT_Axes : public AxesDefinition {
+public:
+   KT_Axes() {}
+
+   virtual std::string short_description() const {
+      return "KT";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "KT Axes";
+      return stream.str();
+   };
+   
+   virtual KT_Axes* create() const {return new KT_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition &) const {
+      return (new AxesFinderFromKT());
+   }
+   
+};
+
+// ca axes
+class CA_Axes : public AxesDefinition {
+public:
+   CA_Axes() {}
+
+   virtual std::string short_description() const {
+      return "CA";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "CA Axes";
+      return stream.str();
+   };
+   
+   virtual CA_Axes* create() const {return new CA_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition &) const {
+      return (new AxesFinderFromCA());
+   }
+   
+};
+
+// anti-kt axes, one parameter R0 is subjet radius
+class AntiKT_Axes : public AxesDefinition {
+
+public:
+   AntiKT_Axes(double R0 = 0.2): _R0(R0) {}
+
+   virtual std::string short_description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "AKT" << _R0;
+      return stream.str();
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Anti-KT Axes (R0 = " << _R0 << ")";
+      return stream.str();
+   };
+   
+   virtual AntiKT_Axes* create() const {return new AntiKT_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition &) const {
+      return (new AxesFinderFromAntiKT(_R0));
+   }
+   
+   
+private:
+   double _R0;
+
+};
+
+// winner-take-all recombination with kt axes
+class WTA_KT_Axes : public AxesDefinition {
+public:
+   WTA_KT_Axes() {}
+
+   virtual std::string short_description() const {
+      return "WTA KT";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Winner-Take-All KT Axes";
+      return stream.str();
+   };
+   
+   virtual WTA_KT_Axes* create() const {return new WTA_KT_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition &) const {
+      return (new AxesFinderFromWTA_KT());
+   }
+
+   
+};
+
+// winner-take-all recombination with CA axes
+class WTA_CA_Axes : public AxesDefinition {
+public:
+   WTA_CA_Axes() {}
+
+   virtual std::string short_description() const {
+      return "WTA CA";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Winner-Take-All CA Axes";
+      return stream.str();
+   };
+   
+   virtual WTA_CA_Axes* create() const {return new WTA_CA_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition &) const {
+      return (new AxesFinderFromWTA_CA());
+   }
+   
+};
+   
+// Onepass minimization from kt axes
+class OnePass_KT_Axes : public AxesDefinition {
+public:
+   OnePass_KT_Axes() {}
+   
+   virtual std::string short_description() const {
+      return "OnePass KT";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "One-Pass Minimization from KT Axes";
+      return stream.str();
+   };
+   
+   virtual OnePass_KT_Axes* create() const {return new OnePass_KT_Axes(*this);}
+   
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition & ) const {
+      return (new AxesFinderFromKT());
+   }
+   virtual AxesFinder* createFinishingAxesFinder(const MeasureDefinition & measure_def) const {
+      return measure_def.createOnePassAxesFinder();
+   }
+
+};
+
+// Onepass minimization from CA axes
+class OnePass_CA_Axes : public AxesDefinition {
+public:
+   OnePass_CA_Axes() {}
+
+   virtual std::string short_description() const {
+      return "OnePass CA";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "One-Pass Minimization from CA Axes";
+      return stream.str();
+   };
+   
+   virtual OnePass_CA_Axes* create() const {return new OnePass_CA_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition & ) const {
+      return (new AxesFinderFromCA());
+   }
+   virtual AxesFinder* createFinishingAxesFinder(const MeasureDefinition & measure_def) const {
+      return measure_def.createOnePassAxesFinder();
+   }
+};
+
+// Onepass minimization from AntiKT axes, one parameter R0
+class OnePass_AntiKT_Axes : public AxesDefinition {
+
+public:
+   OnePass_AntiKT_Axes(double R0): _R0(R0) {}
+   
+   virtual std::string short_description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "OnePassAKT" << _R0;
+      return stream.str();
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "One-Pass Minimization from Anti-KT Axes (R0 = " << _R0 << ")";
+      return stream.str();
+   };
+   
+   virtual OnePass_AntiKT_Axes* create() const {return new OnePass_AntiKT_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition & ) const {
+      return (new AxesFinderFromAntiKT(_R0));
+   }
+   virtual AxesFinder* createFinishingAxesFinder(const MeasureDefinition & measure_def) const {
+      return measure_def.createOnePassAxesFinder();
+   }
+
+private:
+   double _R0;
+};
+
+// Onepass minimization from winner-take-all kt axes
+class OnePass_WTA_KT_Axes : public AxesDefinition {
+public:
+   OnePass_WTA_KT_Axes() {}
+   
+   virtual std::string short_description() const {
+      return "OnePass WTA KT";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "One-Pass Minimization from Winner-Take-All KT Axes";
+      return stream.str();
+   };
+   
+   virtual OnePass_WTA_KT_Axes* create() const {return new OnePass_WTA_KT_Axes(*this);}
+   
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition & ) const {
+      return (new AxesFinderFromWTA_KT());
+   }
+   virtual AxesFinder* createFinishingAxesFinder(const MeasureDefinition & measure_def) const {
+      return measure_def.createOnePassAxesFinder();
+   }
+};
+
+
+// Onepass minimization from winner-take-all CA axes
+class OnePass_WTA_CA_Axes : public AxesDefinition {
+public:
+   OnePass_WTA_CA_Axes() {}
+
+   virtual std::string short_description() const {
+      return "OnePass WTA CA";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "One-Pass Minimization from Winner-Take-All CA Axes";
+      return stream.str();
+   };
+   
+   virtual OnePass_WTA_CA_Axes* create() const {return new OnePass_WTA_CA_Axes(*this);}
+   
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition & ) const {
+      return (new AxesFinderFromWTA_CA());
+   }
+   virtual AxesFinder* createFinishingAxesFinder(const MeasureDefinition & measure_def) const {
+      return measure_def.createOnePassAxesFinder();
+   }
+
+   
+};
+   
+// set axes manually
+class Manual_Axes : public AxesDefinition {
+public:
+   Manual_Axes() {}
+   
+   virtual std::string short_description() const {
+      return "Manual";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Manual Axes";
+      return stream.str();
+   };
+   
+   virtual Manual_Axes* create() const {return new Manual_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return true;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition &) const {
+      return (new AxesFinderFromUserInput());
+   }
+};
+
+// one pass minimization from manual starting point
+class OnePass_Manual_Axes : public AxesDefinition {
+public:
+   OnePass_Manual_Axes() {}
+
+   virtual std::string short_description() const {
+      return "OnePass Manual";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "One-Pass Minimization from Manual Axes";
+      return stream.str();
+   };
+   
+   virtual OnePass_Manual_Axes* create() const {return new OnePass_Manual_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return false;}
+   virtual bool supportsManualAxes() const {return true;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition & ) const {
+      return (new AxesFinderFromUserInput());
+   }
+   virtual AxesFinder* createFinishingAxesFinder(const MeasureDefinition & measure_def) const {
+      return measure_def.createOnePassAxesFinder();
+   }
+
+   
+};
+   
+// multi-pass minimization from kT starting point
+class MultiPass_Axes : public AxesDefinition {
+
+public:
+   MultiPass_Axes(unsigned int Npass) : _Npass(Npass) {}
+
+   virtual std::string short_description() const {
+      return "MultiPass";
+   };
+   
+   virtual std::string description() const {
+      std::stringstream stream;
+      stream << std::fixed << std::setprecision(2)
+      << "Multi-Pass Axes (Npass = " << _Npass << ")";
+      return stream.str();
+   };
+   
+   virtual MultiPass_Axes* create() const {return new MultiPass_Axes(*this);}
+
+   virtual bool givesRandomizedResults() const {return true;}
+   virtual bool supportsManualAxes() const {return false;}
+
+   virtual AxesFinder* createStartingAxesFinder(const MeasureDefinition & ) const {
+      return (new AxesFinderFromKT());
+   }
+   virtual AxesFinder* createFnishingAxesFinder(const MeasureDefinition & measure_def) const {
+      return measure_def.createMultiPassAxesFinder(_Npass);
+   }
+   
+private:
+   unsigned int _Npass;
+   
+};
+   
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif  // __FASTJET_CONTRIB_NJETTINESS_HH__
+
Index: external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.cc	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,26 +30,31 @@
 
 
-// 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,4 +76,7 @@
    }
 
+   //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));
Index: external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,4 +1,2 @@
-// $Id$
-//
 //  Nsubjettiness Package
 //  Questions/Comments?  jthaler@jthaler.net
@@ -7,4 +5,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,21 +48,5 @@
 // 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,9 +57,10 @@
       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,5 +77,21 @@
          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,5 +122,5 @@
  * 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,5 +129,5 @@
  *   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,19 +138,52 @@
 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,6 +197,6 @@
 private:
 
+   Njettiness _njettinessFinder;
    int _N;
-   mutable Njettiness _njettinessFinder; // TODO:  should muck with this so run_clustering can be const without this mutable
 
 };
Index: external/fastjet/contribs/Nsubjettiness/Nsubjettiness.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Nsubjettiness.cc	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/Nsubjettiness.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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.
Index: external/fastjet/contribs/Nsubjettiness/Nsubjettiness.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Nsubjettiness.hh	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/Nsubjettiness.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,5 +33,4 @@
 #include <string>
 #include <climits>
-
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
@@ -49,14 +49,67 @@
 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,7 +120,6 @@
                  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,6 +138,17 @@
    }
    
+   /// 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,5 +158,5 @@
 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,14 +175,47 @@
    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;
 
Index: external/fastjet/contribs/Nsubjettiness/README
===================================================================
--- external/fastjet/contribs/Nsubjettiness/README	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/contribs/Nsubjettiness/README	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,244 @@
+--------------------------------------------------------------------------------
+Nsubjettiness Package
+--------------------------------------------------------------------------------
+
+The Nsubjettiness package is based on the physics described in:
+
+   Identifying Boosted Objects with N-subjettiness.
+   Jesse Thaler and Ken Van Tilburg.
+   JHEP 1103:015 (2011), arXiv:1011.2268.
+
+   Maximizing Boosted Top Identification by Minimizing N-subjettiness. 
+   Jesse Thaler and Ken Van Tilburg. 
+   JHEP 1202:093 (2012), arXiv:1108.2701.
+
+New in v2.0 is the winner-take-all axis, which is 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.
+
+--------------------------------------------------------------------------------
+Core Classes
+--------------------------------------------------------------------------------
+
+There are various ways to access N-(sub)jettiness variables, described
+in more detail below:
+
+Nsubjettiness       [Nsubjettiness.hh]:
+   A FunctionOfPseudoJet<double> interface to measure N-subjettiness
+   (Recommended for most users)
+NsubjettinessRatio  [Nsubjettiness.hh]:
+   A FunctionOfPseudoJet<double> interface to measure ratios of
+   two different N-subjettiness (i.e. tau3/tau2)
+NjettinessPlugin    [NjettinessPlugin.hh]:
+   A FastJet plugin for finding jets by minimizing N-jettiness
+   (Recommended for advanced users)
+Njettiness          [Njettiness.hh]:
+   Access to the core Njettiness code.
+   (Not recommended for users, since the interface might change)
+
+The code assumes that you have FastJet 3.
+
+--------------------------------------------------------------------------------
+Basic Usage: Nsubjettiness and NsubjettinessRatio  [Nsubjettiness.hh]
+--------------------------------------------------------------------------------
+
+Most users will only need to use the Nsubjettiness class.  The basic
+functionality is given by:
+
+   Nsubjettiness nSub(N, AxesDefinition, MeasureDefinition)
+   // N specifies the number of (sub) jets to measure
+   // AxesDefinition is WTA_KT_Axes, OnePass_KT_Axes, etc.
+   // MeasureDefinition is UnnormalizedMeasure(beta),
+   //                      NormalizedMeasure(beta,R0), etc.
+      
+   // get tau value
+   double tauN = nSub.result(PseudoJet);
+
+Also available are ratios of N-subjettiness values
+    NsubjettinessRatio nSubRatio(N, M, AxesDefinition,
+                                 MeasureDefinition)
+    // N and M give tau_N / tau_M, all other options the same
+
+--------------------------------------------------------------------------------
+AxesDefinition  [NjettinessDefinition.hh]
+--------------------------------------------------------------------------------
+
+N-(sub)jettiness requires choosing axes as well as a measure (see below). There
+are a number of axes choices available to the user, though modes with a (*) are
+recommended. Arguments in parentheses are parameters that the user must set.
+
+Axes can be found using standard recursive clustering procedures.  New is the
+option to use the "winner-take-all" recombination scheme:
+(*) KT_Axes           // exclusive kt axes
+    CA_Axes           // exclusive ca axes
+    AntiKT_Axes(R0)   // inclusive hardest axes with antikt, R0 = radius
+(*) WTA_KT_Axes       // exclusive kt with winner-take-all recombination
+    WTA_CA_Axes       // exclusive ca with winner-take-all recombination
+
+One can also run a minimization routine to find a (local) minimum of
+N-(sub)jettiness:
+(*) OnePass_KT_Axes          // one-pass minimization from kt starting point
+    OnePass_CA_Axes          // one-pass min. from ca starting point
+    OnePass_AntiKT(R0)       // one-pass min. from antikt starting point,R0=rad
+(*) OnePass_WTA_KT_Axes      // one-pass min. from wta_kt starting point
+    OnePass_WTA_CA_Axes      // one-pass min. from wta_ca starting point
+
+In general, it is difficult to finding the global minimum, but this mode
+attempts to do so
+    MultiPass_Axes(Npass)    // axes that (attempt to) minimize N-subjettiness
+                             // (NPass = 100 is typical)
+
+Finally, one can set manual axes:
+    Manual_Axes          // set your own axes with setAxes()
+    OnePass_Manual_Axes  // one-pass minimization from manual starting point
+
+For most cases, running with OnePass_KT_Axes or OnePass_WTA_KT_Axes gives
+reasonable results (and the results are IRC safe).  Because it uses random
+number seeds, MultiPass_Axes is not IRC safe (and the code is rather slow).  Note
+that for the minimization routines, beta = 1.1 is faster than beta = 1, with
+comparable performance.
+
+--------------------------------------------------------------------------------
+MeasureDefinition  [NjettinessDefinition.hh]
+--------------------------------------------------------------------------------
+
+At the moment, there are only a few measures.  Note that each one has a
+different number of parameters.  The one indicated by (*)
+is the one recommended for use by users new to Nsubjettiness.
+
+The original N-subjettiness measures are:
+    NormalizedMeasure(beta,R0)  //default normalized measure with
+                                //parameters beta and R0 (dimensionless)
+(*) UnnormalizedMeasure(beta)   //default unnormalized measure with just
+                                //parameter beta (dimensionful)
+
+There are also measures that incorporate a radial cutoff:
+    NormalizedCutoffMeasure(beta,R0,Rcutoff) //normalized measure with
+                                             //additional Rcutoff
+    UnnormalizedCutoffMeasure(beta,Rcutoff)  //unnormalized measure with
+                                             //additional Rcutoff
+
+In beta testing are "geometric" measures where distances are measured using the
+Lorentz dot product (N.B. the formula for the geometric measure is likely to
+change since there should be separate beam and jet beta factors.)
+    GeometricMeasure(beta)               //geometric measure with exponent beta
+    GeometricCutoffMeasure(beta,Rcutoff) //geometric measure with Rcutoff
+
+--------------------------------------------------------------------------------
+A note on beta dependence
+--------------------------------------------------------------------------------
+
+The angular exponent in N-subjettiness is called beta.  The original
+N-subjettiness paper advocated beta = 1, but it is now understood that different
+beta values can be useful in different contexts.  The two main choices are:
+
+beta = 1:  aka broadening/girth/width measure
+   wta_kt_axes are approximately the same as minimizing beta = 1 measure
+
+beta = 2:  aka thrust/mass measure
+   kt_axes are approximately the same as minimizing beta = 2 measure
+
+N.B. The minimization routines are only valid for 1 < beta < 3.
+
+For quark/gluon discrimination with N = 1, beta~0.2 with wta_kt_axes appears
+to be a good choice. 
+
+--------------------------------------------------------------------------------
+TauComponents  [MeasureFunction.hh]
+--------------------------------------------------------------------------------
+
+For most users, they will only need the value of N-subjettiness (i.e. tau)
+itself.  For advanced users, they can access individual tau components (i.e. 
+the individual numerator pieces, the denominator, etc.)  
+
+   TauComponents tauComp = nSub.component_result(jet);
+   vector<double> numer = tauComp.jet_pieces_numerator(); //tau for each subjet
+   double denom = tauComp.denominator();  //normalization factor
+
+--------------------------------------------------------------------------------
+WinnerTakeAllRecombiner  [WinnerTakeAllRecombiner.hh]
+--------------------------------------------------------------------------------
+
+New for version 2.0 of Nsubjettiness are winner-take-all axes.  They are found
+with the help of the WinnerTakeAllRecombiner. This class defines a new
+recombination scheme for clustering particles. This scheme recombines two
+PseudoJets into a PseudoJet with pT of the sum of the two input PseudoJet pTs
+and direction of the harder PseudoJet.  This is a "recoil-free" recombination
+scheme that guarantees that the axes is aligned with one of the input particles.
+It is IRC safe.  Axes found with the standard E-scheme recombiner at similar to
+the beta = 2 minimization, while winner-take-all is similar to the beta = 1
+measure.
+
+Note that the WinnerTakeAllRecombiner can be used outside of Nsubjettiness
+itself for jet finding.  For example, the direction of anti-kT jets found
+with the WinnerTakeAllRecombiner is particularly robust against soft jet
+contamination.
+
+--------------------------------------------------------------------------------
+Advanced Usage:  NjettinessPlugin  [NjettinessPlugin.hh]
+--------------------------------------------------------------------------------
+
+The Njettiness FastJet plugin represents an exclusive jet finder (yielding a
+fixed N number of jets). The algorithm finds N axes, and jets are simply the sum
+of particles closest to a given axis (or unclustered if they are closest to the
+beam).  The axes finding methods and measures are the same as for Nsubjettiness.
+
+   NjettinessPlugin plugin(N, AxesDefinition, MeasureDefinition);
+   JetDefinition def(&plugin);
+   ClusterSequence cs(vector<PseudoJet>,def);
+   vector<PseudoJet> jets = cs.inclusive_jets();
+
+Note that despite being an exclusive jet algorithm, one finds the jets using the
+inclusive_jets() call.
+
+--------------------------------------------------------------------------------
+Very Advanced Usage:  Njettiness  [Njettiness.hh]
+--------------------------------------------------------------------------------
+
+Most users will want to use the Nsubjettiness or NjettinessPlugin classes to
+access N-(sub)jettiness information.  For direct access to the Njettiness class,
+one can use Njettiness.hh directly. This class is still evolving, so users who
+wish to extend its functionality should contact the authors first.
+
+--------------------------------------------------------------------------------
+Technical Details
+--------------------------------------------------------------------------------
+
+In general, the user will never need access to these header files.  Here is a
+brief description about how they are used to help the calculation of
+N-(sub)jettiness:
+
+AxesFinder.hh:
+
+The AxesFinder class (and derived classes) defines the axes used in the
+calculation of N-(sub)jettiness.  These axes can be defined from the exclusive
+jets from a kT or CA algorithm, the hardest jets from an anti-kT algorithm,
+manually, or from minimization of N-jettiness.  In the future, the user will be
+able to write their own axes finder, though currently the interface is still
+evolving.  At the moment, the user should stick to the options allowed by
+AxesDefinition.
+
+MeasureFunction.hh:
+
+The MeasureFunction class (and derived classes) defines the measure by which
+N-(sub)jettiness is calculated. This measure is calculated between each
+particle and its corresponding axis, and then summed and normalized to
+produce N-(sub)jettiness. The default measure for this calculation is
+pT*dR^beta, where dR is the rapidity-azimuth distance between the particle
+and its axis, and beta is the angular exponent.  Again, in the future the user
+will be able to write their own measures, but for the time being, only the
+predefined MeasureDefinition values should be used.
+
+--------------------------------------------------------------------------------
+Known Issues
+--------------------------------------------------------------------------------
+
+-- The MultiPass_Axes mode gives different answers on different runs, since
+   random numbers are used.
+-- In rare cases, one pass minimization can give a larger value of Njettiness
+   than without minimization.
+-- Nsubjettiness is not thread safe, since there are mutables in Njettiness.
Index: external/fastjet/contribs/Nsubjettiness/VERSION
===================================================================
--- external/fastjet/contribs/Nsubjettiness/VERSION	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/VERSION	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,1 +1,1 @@
-1.1.0-beta4
+2.1.0
Index: external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.cc	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,5 +34,4 @@
 
 // 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 
Index: external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh	(revision f1436343d526e830f774e5fde60b4f480d19d6d7)
+++ external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -5,4 +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,5 +51,7 @@
    
    /// 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:
