Index: Makefile
===================================================================
--- Makefile	(revision a0431dcbf4f6030aaef06db0f34e0959b4b9c6ae)
+++ Makefile	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -494,5 +494,9 @@
 	external/fastjet/plugins/SISCone/fastjet/SISConePlugin.hh \
 	external/fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh \
-	external/fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh
+	external/fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh \
+	external/fastjet/contribs/Nsubjettiness/Nsubjettiness.hh \
+	external/fastjet/contribs/Nsubjettiness/Njettiness.hh \
+	external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh \
+	external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh
 tmp/modules/BTagging.$(ObjSuf): \
 	modules/BTagging.$(SrcSuf) \
@@ -1019,4 +1023,16 @@
 tmp/external/Hector/H_OpticalElement.$(ObjSuf): \
 	external/Hector/H_OpticalElement.$(SrcSuf)
+tmp/external/fastjet/contribs/Nsubjettiness/MeasureFunction.$(ObjSuf): \
+	external/fastjet/contribs/Nsubjettiness/MeasureFunction.$(SrcSuf)
+tmp/external/fastjet/contribs/Nsubjettiness/AxesFinder.$(ObjSuf): \
+	external/fastjet/contribs/Nsubjettiness/AxesFinder.$(SrcSuf)
+tmp/external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(ObjSuf): \
+	external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(SrcSuf)
+tmp/external/fastjet/contribs/Nsubjettiness/Nsubjettiness.$(ObjSuf): \
+	external/fastjet/contribs/Nsubjettiness/Nsubjettiness.$(SrcSuf)
+tmp/external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.$(ObjSuf): \
+	external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.$(SrcSuf)
+tmp/external/fastjet/contribs/Nsubjettiness/Njettiness.$(ObjSuf): \
+	external/fastjet/contribs/Nsubjettiness/Njettiness.$(SrcSuf)
 DELPHES_OBJ +=  \
 	tmp/classes/DelphesFormula.$(ObjSuf) \
@@ -1167,5 +1183,11 @@
 	tmp/external/Hector/H_Parameters.$(ObjSuf) \
 	tmp/external/Hector/H_Marker.$(ObjSuf) \
-	tmp/external/Hector/H_OpticalElement.$(ObjSuf)
+	tmp/external/Hector/H_OpticalElement.$(ObjSuf) \
+	tmp/external/fastjet/contribs/Nsubjettiness/MeasureFunction.$(ObjSuf) \
+	tmp/external/fastjet/contribs/Nsubjettiness/AxesFinder.$(ObjSuf) \
+	tmp/external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.$(ObjSuf) \
+	tmp/external/fastjet/contribs/Nsubjettiness/Nsubjettiness.$(ObjSuf) \
+	tmp/external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.$(ObjSuf) \
+	tmp/external/fastjet/contribs/Nsubjettiness/Njettiness.$(ObjSuf)
 
 ifeq ($(HAS_PYTHIA8),true)
@@ -1363,4 +1385,8 @@
 	@touch $@
 
+external/fastjet/contribs/Nsubjettiness/Njettiness.hh: \
+	external/fastjet/PseudoJet.hh
+	@touch $@
+
 modules/Efficiency.h: \
 	classes/DelphesModule.h
@@ -1403,4 +1429,9 @@
 	@touch $@
 
+external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh: \
+	external/fastjet/ClusterSequence.hh \
+	external/fastjet/JetDefinition.hh
+	@touch $@
+
 modules/Cloner.h: \
 	classes/DelphesModule.h
@@ -1506,4 +1537,9 @@
 	@touch $@
 
+external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh: \
+	external/fastjet/PseudoJet.hh \
+	external/fastjet/JetDefinition.hh
+	@touch $@
+
 modules/ImpactParameterSmearing.h: \
 	classes/DelphesModule.h
@@ -1558,4 +1594,8 @@
 	external/fastjet/SharedPtr.hh \
 	external/fastjet/PseudoJetStructureBase.hh
+	@touch $@
+
+external/fastjet/contribs/Nsubjettiness/Nsubjettiness.hh: \
+	external/fastjet/FunctionOfPseudoJet.hh
 	@touch $@
 
Index: classes/DelphesClasses.cc
===================================================================
--- classes/DelphesClasses.cc	(revision a0431dcbf4f6030aaef06db0f34e0959b4b9c6ae)
+++ classes/DelphesClasses.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -112,4 +112,5 @@
   Position(0.0, 0.0, 0.0, 0.0),
   Dxy(0),SDxy(0),Xd(0),Yd(0),Zd(0),
+  Tau1(0.0), Tau2(0.0), Tau3(0.0), Tau4(0.0), Tau5(0.0), 
   Area(0.0, 0.0, 0.0, 0.0),
   NCharged(0),
@@ -235,5 +236,10 @@
   object.FracPt[3] = FracPt[3];
   object.FracPt[4] = FracPt[4];
-
+  object.Tau1 = Tau1;
+  object.Tau2 = Tau2;
+  object.Tau3 = Tau3;
+  object.Tau4 = Tau4;
+  object.Tau5 = Tau5;
+ 
   object.fFactory = fFactory;
   object.fArray = 0;
@@ -292,5 +298,10 @@
   FracPt[3] = 0.0;
   FracPt[4] = 0.0;
-
+  Tau1 = 0.0;
+  Tau2 = 0.0;
+  Tau3 = 0.0;
+  Tau4 = 0.0;
+  Tau5 = 0.0;
+ 
   fArray = 0;
 }
Index: classes/DelphesClasses.h
===================================================================
--- classes/DelphesClasses.h	(revision a0431dcbf4f6030aaef06db0f34e0959b4b9c6ae)
+++ classes/DelphesClasses.h	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -318,5 +318,11 @@
   Float_t  FracPt[5];
 
-
+ // -- Nsubjettiness variables ---
+  
+  Float_t Tau1; 
+  Float_t Tau2; 
+  Float_t Tau3; 
+  Float_t Tau4; 
+  Float_t Tau5; 
 
   ClassDef(Jet, 2)
@@ -466,5 +472,13 @@
   Float_t  PTD;
   Float_t  FracPt[5];
-
+ 
+  // -- Nsubjettiness variables ---
+ 
+  Float_t Tau1;
+  Float_t Tau2;
+  Float_t Tau3;
+  Float_t Tau4;
+  Float_t Tau5;
+ 
   static CompBase *fgCompare; //!
   const CompBase *GetCompare() const { return fgCompare; }
Index: doc/genMakefile.tcl
===================================================================
--- doc/genMakefile.tcl	(revision a0431dcbf4f6030aaef06db0f34e0959b4b9c6ae)
+++ doc/genMakefile.tcl	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -273,5 +273,5 @@
 dictDeps {DISPLAY_DICT} {display/DisplayLinkDef.h}
 
-sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/fastjet/*.cc} {external/fastjet/tools/*.cc} {external/fastjet/plugins/*/*.cc} {external/Hector/*.cc}
+sourceDeps {DELPHES} {classes/*.cc} {modules/*.cc} {external/ExRootAnalysis/*.cc} {external/fastjet/*.cc} {external/fastjet/tools/*.cc} {external/fastjet/plugins/*/*.cc} {external/Hector/*.cc} {external/fastjet/contribs/Nsubjettiness/*.cc}
 
 sourceDeps {DISPLAY} {display/*.cc}
Index: external/fastjet/contribs/Nsubjettiness/AxesFinder.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/AxesFinder.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/AxesFinder.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,311 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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 "AxesFinder.hh"
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+///////
+//
+// Functions for minimization.
+//
+///////
+
+// Given starting axes, update to find better axes by using Kmeans clustering around the old axes
+template <int N>
+std::vector<LightLikeAxis> AxesFinderFromOnePassMinimization::UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes, 
+                                  const std::vector <fastjet::PseudoJet> & inputJets) {
+   assert(old_axes.size() == N);
+   
+   // some storage, declared static to save allocation/re-allocation costs
+   static LightLikeAxis new_axes[N];
+   static fastjet::PseudoJet new_jets[N];
+   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
+   }
+
+   double precision = _precision;
+   
+   /////////////// Assignment Step //////////////////////////////////////////////////////////
+   std::vector<int> assignment_index(inputJets.size()); 
+   int k_assign = -1;
+   
+   for (unsigned i = 0; i < inputJets.size(); i++){
+      double smallestDist = std::numeric_limits<double>::max();  //large number
+      for (int k = 0; k < N; k++) {
+         double thisDist = old_axes[k].DistanceSq(inputJets[i]);
+         if (thisDist < smallestDist) {
+            smallestDist = thisDist;
+            k_assign = k;
+         }
+      }
+      if (smallestDist > sq(_Rcutoff)) {k_assign = -1;}
+      assignment_index[i] = k_assign;
+   }
+   
+   //////////////// Update Step /////////////////////////////////////////////////////////////
+   double distPhi, old_dist;
+   for (unsigned i = 0; i < inputJets.size(); i++) {
+      int old_jet_i = assignment_index[i];
+      if (old_jet_i == -1) {continue;}
+
+      const fastjet::PseudoJet& inputJet_i = inputJets[i];
+      LightLikeAxis& new_axis_i = new_axes[old_jet_i];
+      double inputPhi_i = inputJet_i.phi();
+      double inputRap_i = inputJet_i.rap();
+            
+      // optimize pow() call
+      // add noise (the precision term) to make sure we don't divide by zero
+      if (_beta == 1.0) {
+         double DR = std::sqrt(sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i));
+         old_dist = 1.0/DR;
+      } else if (_beta == 2.0) {
+         old_dist = 1.0;
+      } else if (_beta == 0.0) {
+         double DRSq = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i);
+         old_dist = 1.0/DRSq;
+      } else {
+         old_dist = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i);
+         old_dist = std::pow(old_dist, (0.5*_beta-1.0));
+      }
+      
+      // TODO:  Put some of these addition functions into light-like axes
+      // rapidity sum
+      new_axis_i.set_rap(new_axis_i.rap() + inputJet_i.perp() * inputRap_i * old_dist);
+      // phi sum
+      distPhi = inputPhi_i - old_axes[old_jet_i].phi();
+      if (fabs(distPhi) <= M_PI){
+         new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * inputPhi_i * old_dist );
+      } else if (distPhi > M_PI) {
+         new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (-2*M_PI + inputPhi_i) * old_dist );
+      } else if (distPhi < -M_PI) {
+         new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (+2*M_PI + inputPhi_i) * old_dist );
+      }
+      // weights sum
+      new_axis_i.set_weight( new_axis_i.weight() + inputJet_i.perp() * old_dist );
+      // momentum magnitude sum
+      new_jets[old_jet_i] += inputJet_i;
+   }
+   // normalize sums
+   for (int k = 0; k < N; k++) {
+      if (new_axes[k].weight() == 0) {
+         // no particles were closest to this axis!  Return to old axis instead of (0,0,0,0)
+         new_axes[k] = old_axes[k];
+      } else {
+         new_axes[k].set_rap( new_axes[k].rap() / new_axes[k].weight() );
+         new_axes[k].set_phi( new_axes[k].phi() / new_axes[k].weight() );
+         new_axes[k].set_phi( std::fmod(new_axes[k].phi() + 2*M_PI, 2*M_PI) );
+         new_axes[k].set_mom( std::sqrt(new_jets[k].modp2()) );
+      }
+   }
+   std::vector<LightLikeAxis> new_axes_vec(N);
+   for (unsigned k = 0; k < N; ++k) new_axes_vec[k] = new_axes[k];
+   return new_axes_vec;
+}
+
+// Given N starting axes, this function updates all axes to find N better axes. 
+// (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) {
+   int N = old_axes.size();
+   switch (N) {
+      case 1: return UpdateAxesFast<1>(old_axes, inputJets);
+      case 2: return UpdateAxesFast<2>(old_axes, inputJets);
+      case 3: return UpdateAxesFast<3>(old_axes, inputJets);
+      case 4: return UpdateAxesFast<4>(old_axes, inputJets);
+      case 5: return UpdateAxesFast<5>(old_axes, inputJets);
+      case 6: return UpdateAxesFast<6>(old_axes, inputJets);
+      case 7: return UpdateAxesFast<7>(old_axes, inputJets);
+      case 8: return UpdateAxesFast<8>(old_axes, inputJets);
+      case 9: return UpdateAxesFast<9>(old_axes, inputJets);
+      case 10: return UpdateAxesFast<10>(old_axes, inputJets);
+      case 11: return UpdateAxesFast<11>(old_axes, inputJets);
+      case 12: return UpdateAxesFast<12>(old_axes, inputJets);
+      case 13: return UpdateAxesFast<13>(old_axes, inputJets);
+      case 14: return UpdateAxesFast<14>(old_axes, inputJets);
+      case 15: return UpdateAxesFast<15>(old_axes, inputJets);
+      case 16: return UpdateAxesFast<16>(old_axes, inputJets);
+      case 17: return UpdateAxesFast<17>(old_axes, inputJets);
+      case 18: return UpdateAxesFast<18>(old_axes, inputJets);
+      case 19: return UpdateAxesFast<19>(old_axes, inputJets);
+      case 20: return UpdateAxesFast<20>(old_axes, inputJets);
+      default: std::cout << "N-jettiness is hard-coded to only allow up to 20 jets!" << std::endl;
+         return std::vector<LightLikeAxis>();
+   }
+
+}
+
+// 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) {
+	  
+   // convert from PseudoJets to LightLikeAxes
+   std::vector< LightLikeAxis > old_axes(n_jets, LightLikeAxis(0,0,0,0));
+   for (int k = 0; k < n_jets; k++) {
+      old_axes[k].set_rap( seedAxes[k].rap() );
+      old_axes[k].set_phi( seedAxes[k].phi() );
+   }
+   
+   // Find new axes by iterating (only one pass here)
+   std::vector< LightLikeAxis > new_axes(n_jets, LightLikeAxis(0,0,0,0));
+   double cmp = std::numeric_limits<double>::max();  //large number
+   int h = 0;
+   while (cmp > _precision && h < _halt) { // Keep updating axes until near-convergence or too many update steps
+      cmp = 0.0;
+      h++;
+      new_axes = UpdateAxes(old_axes, inputJets); // Update axes
+      for (int k = 0; k < n_jets; k++) {
+         cmp += old_axes[k].Distance(new_axes[k]);
+      }
+      cmp = cmp / ((double) n_jets);
+      old_axes = new_axes;
+   }
+      
+   // Convert from internal LightLikeAxes to PseudoJet
+   std::vector<fastjet::PseudoJet> outputAxes;
+   for (int k = 0; k < n_jets; k++) {
+      fastjet::PseudoJet temp = old_axes[k].ConvertToPseudoJet();
+      outputAxes.push_back(temp);
+   }
+   
+   // this is used to debug the minimization routine to make sure that it works.
+   bool do_debug = false;
+   if (do_debug) {
+      // get this information to make sure that minimization is working properly
+      TauComponents seed_tau_components = _measureFunction.result(inputJets, seedAxes);
+      double seed_tau = seed_tau_components.tau();
+      TauComponents tau_components = _measureFunction.result(inputJets, outputAxes);
+      double outputTau = tau_components.tau();
+      assert(outputTau <= seed_tau);
+   }
+   
+   return outputAxes;
+}
+
+PseudoJet AxesFinderFromKmeansMinimization::jiggle(const PseudoJet& axis) {
+   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;
+   
+   double new_phi = axis.phi() + phi_noise;
+   if (new_phi >= 2.0*M_PI) new_phi -= 2.0*M_PI;
+   if (new_phi <= -2.0*M_PI) new_phi += 2.0*M_PI;
+
+   PseudoJet newAxis(0,0,0,0);
+   newAxis.reset_PtYPhiM(axis.perp(),axis.rap() + rap_noise,new_phi);
+   return newAxis;
+}
+   
+   
+// 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) {
+   
+   // first iteration
+	std::vector<fastjet::PseudoJet> bestAxes = _onePassFinder.getAxes(n_jets, inputJets, seedAxes);
+   
+   double bestTau = (_measureFunction.result(inputJets,bestAxes)).tau();
+   
+   for (int l = 1; l < _n_iterations; l++) { // Do minimization procedure multiple times (l = 1 to start since first iteration is done already)
+   
+      // Add noise to current best axes
+      std::vector< PseudoJet > noiseAxes(n_jets, PseudoJet(0,0,0,0));
+      for (int k = 0; k < n_jets; k++) {
+         noiseAxes[k] = jiggle(bestAxes[k]);
+      }
+
+      std::vector<fastjet::PseudoJet> testAxes = _onePassFinder.getAxes(n_jets, inputJets, noiseAxes);
+      double testTau = (_measureFunction.result(inputJets,testAxes)).tau();
+      
+      if (testTau < bestTau) {
+         bestTau = testTau;
+         bestAxes = testAxes;
+      }
+   }
+   
+   return bestAxes;
+}
+
+// Uses minimization of the geometric distance in order to find the minimum axes.
+// 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> seedAxes = currentAxes;
+   double seedTau = _function->tau(particles, seedAxes);
+   
+   for (int i = 0; i < _nAttempts; i++) {
+      
+      std::vector<fastjet::PseudoJet> newAxes(seedAxes.size(),fastjet::PseudoJet(0,0,0,0));
+      
+      // find closest axis and assign to that
+      for (unsigned int i = 0; i < particles.size(); i++) {
+         
+         // start from unclustered beam measure
+         int minJ = -1;
+         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]);
+            if (tempDist < minDist) {
+               minDist = tempDist;
+               minJ = j;
+            }
+         }
+         
+         // if not unclustered, then cluster
+         if (minJ != -1) newAxes[minJ] += particles[i];
+      }
+      
+      // calculate tau on new axes
+      seedAxes = newAxes;
+      double tempTau = _function->tau(particles, newAxes);
+      
+      // close enough to stop?
+      if (fabs(tempTau - seedTau) < _accuracy) break;
+      seedTau = tempTau;
+   }
+   
+   return seedAxes;
+}
+
+// Go from internal LightLikeAxis to PseudoJet
+fastjet::PseudoJet LightLikeAxis::ConvertToPseudoJet() {
+    double px, py, pz, E;
+    E = _mom;
+    pz = (std::exp(2.0*_rap) - 1.0) / (std::exp(2.0*_rap) + 1.0) * E;
+    px = std::cos(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) );
+    py = std::sin(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) );
+    return fastjet::PseudoJet(px,py,pz,E);
+}
+
+} //namespace contrib
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/Nsubjettiness/AxesFinder.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/AxesFinder.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/AxesFinder.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,407 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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_AXESFINDER_HH__
+#define __FASTJET_CONTRIB_AXESFINDER_HH__
+
+#include "WinnerTakeAllRecombiner.hh"
+#include "MeasureFunction.hh"
+
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/JetDefinition.hh"
+
+#include <cmath>
+#include <vector>
+#include <list>
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+///////
+//
+// Axes Finder Options
+//
+///////
+
+//------------------------------------------------------------------------
+/// \class AxesFinder
+// This is the base class for all axes finders. These axes are used along with the MeasureFunctions to calculate 
+// tau_N. There are different implementations of axes finding that are defined in derived classes below.
+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;
+   
+};
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromExclusiveJetDefinition
+// This class finds axes by clustering the particles and then finding the exclusive jets. This can be implemented
+// 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);
+      }
+};
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromWTA_KT
+// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a 
+// 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;}
+   };
+   
+//------------------------------------------------------------------------
+/// \class AxesFinderFromWTA_CA
+// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a 
+// 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;}
+//};
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromKT
+// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a 
+// 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)) {}
+};
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromCA
+// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a 
+// 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)) {}
+};
+
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromHardestJetDefinition
+// This class finds axes by clustering the particles and then finding the n hardest inclusive jets. 
+// 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;
+      }      
+};
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromAntiKT
+// This class finds axes by finding the n hardest jets after clustering the particles according 
+// 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)) {}
+};
+
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromUserInput
+// This class allows the user to manually define the axes.
+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;
+      }
+};
+
+//This is a helper class for the Minimum Axes Finders. It is defined later.
+class LightLikeAxis;                                          
+
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromOnePassMinimization
+// This class defines an AxesFinder that uses Kmeans minimization, but only on a single pass.
+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);
+
+};
+
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromKmeansMinimization
+// This class finds finds axes by using Kmeans clustering to minimizaiton N-jettiness. Given a first set of 
+// starting axes, it updates n times to get as close to the global minimum as possible. This class calls OnePass many times,
+// added noise to the axes.
+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);
+
+};
+
+//------------------------------------------------------------------------
+/// \class AxesFinderFromGeometricMinimization
+// This class finds axes by minimizing the Lorentz dot product distance between axes and particles. Given a first set of starting axes,
+// it essentially does stable cone finxing.
+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);
+      }
+
+      ~AxesFinderFromGeometricMinimization() {
+         delete _function;
+      }
+   
+      virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes);
+};
+
+//------------------------------------------------------------------------
+/// \class LightLikeAxis
+// This is a helper class for the minimum Axes Finders classes above. It creates a convenient way of defining axes
+// 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) {}
+   LightLikeAxis(double my_rap, double my_phi, double my_weight, double my_mom) :
+   _rap(my_rap), _phi(my_phi), _weight(my_weight), _mom(my_mom) {}
+   double rap() const {return _rap;}
+   double phi() const {return _phi;}
+   double weight() const {return _weight;}
+   double mom() const {return _mom;}
+   void set_rap(double my_set_rap) {_rap = my_set_rap;}
+   void set_phi(double my_set_phi) {_phi = my_set_phi;}
+   void set_weight(double my_set_weight) {_weight = my_set_weight;}
+   void set_mom(double my_set_mom) {_mom = my_set_mom;}
+   void reset(double my_rap, double my_phi, double my_weight, double my_mom) {_rap=my_rap; _phi=my_phi; _weight=my_weight; _mom=my_mom;}
+
+   // return PseudoJet with information
+   fastjet::PseudoJet ConvertToPseudoJet();
+   
+   double DistanceSq(const fastjet::PseudoJet& input) const {
+      return DistanceSq(input.rap(),input.phi());
+   }
+
+   double Distance(const fastjet::PseudoJet& input) const {
+      return std::sqrt(DistanceSq(input));
+   }
+
+   double DistanceSq(const LightLikeAxis& input) const {
+      return DistanceSq(input.rap(),input.phi());
+   }
+
+   double Distance(const LightLikeAxis& input) const {
+      return std::sqrt(DistanceSq(input));
+   }
+
+};
+
+} //namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif  // __FASTJET_CONTRIB_AXESFINDER_HH__
Index: external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/MeasureFunction.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,84 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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 "MeasureFunction.hh"
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+///////
+//
+// Measure Function
+//
+///////
+
+// 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) {
+
+   std::vector<double> jetPieces(axes.size(), 0.0);
+   double beamPiece = 0.0;
+   
+   // 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)
+   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) {
+         if (_has_beam) beamPiece += beam_numerator(particles[i]);
+         else assert(_has_beam);  // this should never happen.
+      } else {
+         jetPieces[j_min] += jet_numerator(particles[i],axes[j_min]);
+      }
+   }
+   
+   // 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
+   }
+   
+   return TauComponents(jetPieces, beamPiece, tauDen, _has_denominator, _has_beam);
+}
+   
+} //namespace contrib
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/Nsubjettiness/MeasureFunction.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/MeasureFunction.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/MeasureFunction.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,276 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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_MEASUREFUNCTION_HH__
+#define __FASTJET_CONTRIB_MEASUREFUNCTION_HH__
+
+#include "fastjet/PseudoJet.hh"
+#include <cmath>
+#include <vector>
+#include <list>
+#include <limits>
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+inline double sq(double x) {return x*x;}
+
+///////
+//
+// Measure Function
+//
+///////
+
+/// \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:
+   // empty constructor necessary to initialize tau_components in Njettiness
+   // later set correctly in Njettiness::getTau function
+   TauComponents() {
+      _jet_pieces_numerator.resize(1, 0.0);
+      _beam_piece_numerator = 0.0;
+      _denominator = 0;
+      _numerator = 0;
+      _jet_pieces.resize(1, 0.0);
+      _beam_piece = 0.0;
+      _tau = 0;
+      _has_denominator = false;
+      _has_beam = false;
+   }
+   
+   // This constructor takes input vector and double and calculates all necessary tau components
+   TauComponents(std::vector<double> jet_pieces_numerator, double beam_piece_numerator, double denominator, bool has_denominator, bool has_beam)
+   : _jet_pieces_numerator(jet_pieces_numerator),
+   _beam_piece_numerator(beam_piece_numerator),
+   _denominator(denominator),
+   _has_denominator(has_denominator),
+   _has_beam(has_beam) {
+      
+      if (!_has_denominator) assert(_denominator == 1.0); //make sure no effect from _denominator if _has_denominator is false
+      if (!_has_beam) assert (_beam_piece_numerator == 0.0); //make sure no effect from _beam_piece_numerator if _has_beam is false
+      
+      _numerator = _beam_piece_numerator;
+      _jet_pieces.resize(_jet_pieces_numerator.size(),0.0);
+      for (unsigned j = 0; j < _jet_pieces_numerator.size(); j++) {
+         _jet_pieces[j] = _jet_pieces_numerator[j]/_denominator;
+         _numerator += _jet_pieces_numerator[j];
+      }
+      
+      _beam_piece = _beam_piece_numerator/_denominator;
+      _tau = _numerator/_denominator;
+   }
+   
+   
+   // return values
+   std::vector<double> jet_pieces_numerator() const { return _jet_pieces_numerator; }
+   double beam_piece_numerator() const { return _beam_piece_numerator; }
+   double denominator() const { return _denominator; }
+   double numerator() const { return _numerator; }
+
+   bool has_denominator() const { return _has_denominator; }
+   bool has_beam() const { return _has_beam; }
+   
+   std::vector<double> jet_pieces() const { return _jet_pieces; }
+   double beam_piece() const { return _beam_piece; }
+   double tau() const { return _tau; }
+   
+};
+
+//------------------------------------------------------------------------
+/// \class MeasureFunction
+// This class calculates the tau_N of a jet given a specific measure.
+// It serves as a base class for calculating tau_N according to different measures that are implemented 
+// in further derived classes, but does not define a particular measure itself.
+class MeasureFunction {
+   
+protected:
+   //bool set by derived classes to choose whether or not to use the denominator
+   bool _has_denominator;
+   bool _has_beam;
+   
+   // This constructor allows _has_denominator to be set by derived classes
+   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
+// 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 DefaultUnnormalizedMeasure
+// 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 GeometricMeasure
+// 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;
+      }
+};
+
+
+} //namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif  // __FASTJET_CONTRIB_MEASUREFUNCTION_HH__
Index: external/fastjet/contribs/Nsubjettiness/Njettiness.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Njettiness.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/Njettiness.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,261 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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 "Njettiness.hh"
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib {
+
+
+///////
+//
+// Main Njettiness Class
+//
+///////
+
+// 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) {
+
+   // 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;
+   // The normalized measures have an R0 value.
+   double R0 = NAN;
+
+   // Find the MeasureFunction and set the parameters.
+   switch (measure_mode) {
+      case normalized_measure:
+         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); }
+         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); }
+         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); }
+         break;
+      case normalized_cutoff_measure:
+         beta = para1;
+         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); }
+         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); }
+         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); }
+         break;
+      default:
+         assert(false);
+         break;
+   }   
+
+   // Choose which AxesFinder from user input.
+   // Uses setOnePassAxesFinder helpful function to use beta and Rcutoff values about (if needed)
+   switch (axes_mode) {
+      case wta_kt_axes:
+         _axesFinder = new AxesFinderFromWTA_KT(); 
+         break;
+      case wta_ca_axes:
+         _axesFinder = new AxesFinderFromWTA_CA(); 
+         break;
+      case kt_axes:
+         _axesFinder = new AxesFinderFromKT();
+         break;
+      case ca_axes:
+         _axesFinder = new AxesFinderFromCA();
+         break;
+      case antikt_0p2_axes:
+         _axesFinder = new AxesFinderFromAntiKT(0.2);     
+         break;
+      case onepass_wta_kt_axes:
+         setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA_KT(), beta, Rcutoff);
+         break;
+      case onepass_wta_ca_axes:
+         setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA_CA(), beta, Rcutoff);
+         break;
+      case onepass_kt_axes:
+         setOnePassAxesFinder(measure_mode, new AxesFinderFromKT(), beta, Rcutoff);
+         break;
+      case onepass_ca_axes:
+         setOnePassAxesFinder(measure_mode, new AxesFinderFromCA(), beta, Rcutoff);
+         break;
+      case onepass_antikt_0p2_axes:
+         setOnePassAxesFinder(measure_mode, new AxesFinderFromAntiKT(0.2), beta, Rcutoff);
+         break;
+      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;
+      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;
+      default:
+         assert(false);
+         break;
+      }   
+
+}
+
+// setAxes for Manual mode
+void Njettiness::setAxes(std::vector<fastjet::PseudoJet> myAxes) {
+   if (_current_axes_mode == manual_axes || _current_axes_mode == onepass_manual_axes) {
+      _currentAxes = myAxes;
+   }
+   else {
+      std::cerr << "You can only use setAxes if using manual_axes or onepass_manual_axes measure mode" << std::endl;
+      exit(1);
+   }
+}
+   
+// 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) {
+   if (inputJets.size() <= n_jets) {  //if not enough particles, return zero
+      _currentAxes = inputJets;
+      _currentAxes.resize(n_jets,fastjet::PseudoJet(0.0,0.0,0.0,0.0));
+      _current_tau_components = TauComponents();
+      _seedAxes = _currentAxes;
+   } 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
+   }
+   return _current_tau_components;
+}
+   
+   
+// Partition a list of particles according to which N-jettiness axis they are closest to.
+// Return a vector of length _currentAxes.size() (which should be N).
+// 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;
+}
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/Nsubjettiness/Njettiness.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Njettiness.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/Njettiness.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,183 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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_HH__
+#define __FASTJET_CONTRIB_NJETTINESS_HH__
+
+#include "MeasureFunction.hh"
+#include "AxesFinder.hh"
+
+#include "fastjet/PseudoJet.hh"
+#include <cmath>
+#include <vector>
+#include <list>
+
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib {
+
+///////
+//
+// Main Njettiness Class
+//
+///////
+
+//------------------------------------------------------------------------
+/// \class Njettiness
+// Njettiness uses AxesFinder and MeasureFunction together in order to find tau_N for the event. The user specifies
+// which AxesFinder and which MeasureFunction to use in the calculation, and then Njettiness returns tau_N for the event.
+// It also can return information about the axes and jets it used in the calculation, as well as information about 
+// how the event was partitioned.
+class Njettiness {
+public:
+   
+   // The various axes choices available to the user
+   enum AxesMode {
+      kt_axes,             // exclusive kt axes
+      ca_axes,             // exclusive ca axes
+      antikt_0p2_axes,     // inclusive hardest axes with antikt-0.2
+      wta_kt_axes,         // Winner Take All axes with kt
+      wta_ca_axes,         // Winner Take All axes with CA
+      onepass_kt_axes,     // one-pass minimization from kt starting point
+      onepass_ca_axes,     // one-pass minimization from ca starting point
+      onepass_antikt_0p2_axes,  // one-pass minimization from antikt-0.2 starting point
+      onepass_wta_kt_axes, //one-pass minimization of WTA axes with kt 
+      onepass_wta_ca_axes, //one-pass minimization of WTA axes with ca 
+      min_axes,            // axes that minimize N-subjettiness (100 passes by default)
+      manual_axes,         // set your own axes with setAxes()
+      onepass_manual_axes  // one-pass minimization from manual starting point
+      //  These options are commented out because they have not been fully tested
+      //      wta2_kt_axes,        // Winner Take All (alpha = 2) with kt
+      //      wta2_ca_axes,         // Winner Take All (alpha = 2) with CA
+      //      onepass_wta2_kt_axes, //one-pass minimization of WTA (alpha = 2) axes with kt
+      //      onepass_wta2_ca_axes, //one-pass minimization of WTA (alpha = 2) axes with ca
+   };
+
+   // The measures available to the user.
+   // "normalized_cutoff_measure" was the default in v1.0 of Nsubjettiness
+   // "unnormalized_measure" is now the recommended default usage
+   enum MeasureMode {
+      normalized_measure,           //default normalized measure
+      unnormalized_measure,         //default unnormalized measure
+      geometric_measure,            //geometric measure
+      normalized_cutoff_measure,    //default normalized measure with explicit Rcutoff
+      unnormalized_cutoff_measure,  //default unnormalized measure with explicit Rcutoff
+      geometric_cutoff_measure      //geometric measure with explicit Rcutoff
+   };
+
+private:
+   // The chosen axes/measure modes
+   AxesFinder* _axesFinder;  // The chosen axes
+   MeasureFunction* _measureFunction; // The chosen measure
+
+   // 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; }
+
+   // 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;
+   }
+
+   // 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;
+   }
+   
+   // setAxes for Manual mode
+   void setAxes(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);
+
+   // 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) {
+      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;}
+   
+   // Return axes found by getTauComponents.
+   std::vector<fastjet::PseudoJet> currentAxes() { return _currentAxes;}
+   // Return seedAxes used if onepass minimization (otherwise, same as currentAxes)
+   std::vector<fastjet::PseudoJet> seedAxes() { return _seedAxes;}
+   
+   // partition inputs by Voronoi (each vector stores indices corresponding to inputJets)
+   std::vector<std::list<int> > getPartition(const std::vector<fastjet::PseudoJet> & inputJets);
+
+   // partition inputs by Voronoi
+   std::vector<fastjet::PseudoJet> getJets(const std::vector<fastjet::PseudoJet> & inputJets);
+
+};
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif  // __FASTJET_CONTRIB_NJETTINESS_HH__
+
Index: external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,80 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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 "NjettinessPlugin.hh"
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+
+// 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.
+void NjettinessPlugin::run_clustering(ClusterSequence& cs) const
+{
+   std::vector<fastjet::PseudoJet> particles = cs.jets();
+   _njettinessFinder.getTau(_N, particles);
+   std::vector<std::list<int> > partition = _njettinessFinder.getPartition(particles);
+
+   std::vector<fastjet::PseudoJet> jet_indices_for_extras;
+
+   // output clusterings for each jet
+   for (size_t i = 0; i < partition.size(); ++i) {
+      std::list<int>& indices = partition[i];
+      if (indices.size() == 0) continue;
+      while (indices.size() > 1) {
+         int merge_i = indices.back(); indices.pop_back();
+         int merge_j = indices.back(); indices.pop_back();
+         int newIndex;
+         double fakeDij = -1.0;
+      
+         cs.plugin_record_ij_recombination(merge_i, merge_j, fakeDij, newIndex);
+
+         indices.push_back(newIndex);
+      }
+      double fakeDib = -1.0;
+      
+      int finalJet = indices.back();
+      cs.plugin_record_iB_recombination(finalJet, fakeDib);
+      jet_indices_for_extras.push_back(cs.jets()[finalJet]);  // Get the four vector for the final jets to compare later.
+   }
+
+   NjettinessExtras * extras = new NjettinessExtras(_njettinessFinder.currentTauComponents(),jet_indices_for_extras,_njettinessFinder.currentAxes());
+   cs.plugin_associate_extras(std::auto_ptr<ClusterSequence::Extras>(extras));
+   
+}
+
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,175 @@
+// $Id$
+//
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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_NJETTINESSPLUGIN_HH__
+#define __FASTJET_CONTRIB_NJETTINESSPLUGIN_HH__
+
+#include "Njettiness.hh"
+#include "MeasureFunction.hh"
+#include "AxesFinder.hh"
+
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/JetDefinition.hh"
+
+#include <string>
+#include <climits>
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+
+namespace contrib {
+
+//------------------------------------------------------------------------
+/// \class NjettinessExtras
+// This class contains the same information as Njettiness, but redoes it in terms of the ClusterSequence::Extras class.
+// This is done in order to help improve the interface for the main NjettinessPlugin class.
+// TODO:  This class should probably be merged with TauComponents, since both have access
+// 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) {}
+      
+      double totalTau() const {return _tau_components.tau();}
+      std::vector<double> subTaus() const {return _tau_components.jet_pieces();}
+      std::vector<fastjet::PseudoJet> jets() const {return _jets;}
+      std::vector<fastjet::PseudoJet> axes() const {return _axes;}
+      
+      double totalTau(const fastjet::PseudoJet& jet) const {
+         return _tau_components.tau();
+      }
+      double subTau(const fastjet::PseudoJet& jet) const {
+         if (labelOf(jet) == -1) return NAN;
+         return _tau_components.jet_pieces()[labelOf(jet)];
+      }
+      
+      double beamTau() const {
+         return _tau_components.beam_piece();
+      }
+      
+      fastjet::PseudoJet axis(const fastjet::PseudoJet& jet) const {
+         return _axes[labelOf(jet)];
+      }
+
+      bool has_njettiness_extras(const fastjet::PseudoJet& jet) const {
+         return (labelOf(jet) >= 0);
+      }
+
+};
+
+inline const NjettinessExtras * njettiness_extras(const fastjet::PseudoJet& jet) {
+   const ClusterSequence * myCS = jet.associated_cluster_sequence();   
+   if (myCS == NULL) return NULL;
+   const NjettinessExtras* extras = dynamic_cast<const NjettinessExtras*>(myCS->extras());   
+   return extras;   
+}
+
+inline const NjettinessExtras * njettiness_extras(const fastjet::ClusterSequence& myCS) {
+   const NjettinessExtras* extras = dynamic_cast<const NjettinessExtras*>(myCS.extras());   
+   return extras;   
+}
+
+/// The Njettiness jet algorithm
+/**
+ * An exclusive jet finder that identifies N jets; first N axes are found, then
+ * particles are assigned to the nearest (DeltaR) axis and for each axis the
+ * corresponding jet is simply the four-momentum sum of these particles.
+ *
+ * Axes can be found in several ways, specified by the AxesMode argument.  The
+ * recommended choices are
+ *
+ * kt_axes              : exclusive kT
+ * wta_kt_axes          : exclusive kT with winner-take-all-recombination
+ * onepass_kt_axes      : one-pass minimization seeded by kt (pretty good)
+ * onepass_wta_kt_axes  : one-pass minimization seeded by wta_kt
+ *
+ * For the unnormalized_measure, N-jettiness is defined as:
+ *
+ * tau_N = Sum_{all particles i} p_T^i min((DR_i1)^beta, (DR_i2)^beta, ...)
+ *
+ *   DR_ij is the distance sqrt(Delta_phi^2 + Delta_rap^2) between particle i
+ *   and jet j.
+ * 
+ * The normalized_meausure 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.
+ *
+ */
+
+class NjettinessPlugin : public JetDefinition::Plugin {
+public:
+
+   NjettinessPlugin(int N,
+                    Njettiness::AxesMode axes_mode,
+                    Njettiness::MeasureMode measure_mode,
+                    double para1 = NAN,
+                    double para2 = NAN,
+                    double para3 = NAN,
+                    double para4 = NAN);
+
+   // Old constructor for backwards compatibility with v1.0,
+   // where normalized_cutoff_measure was the only option
+   NjettinessPlugin(int N,
+                    Njettiness::AxesMode mode,
+                    double beta,
+                    double R0,
+                    double Rcutoff=std::numeric_limits<double>::max());
+
+
+   // The things that are required by base class.
+   virtual std::string description () const;
+   virtual double R() const {return -1.0;} // TODO: make this not stupid
+   virtual void run_clustering(ClusterSequence&) const;
+
+   virtual ~NjettinessPlugin() {}
+
+private:
+
+   int _N;
+   mutable Njettiness _njettinessFinder; // TODO:  should muck with this so run_clustering can be const without this mutable
+
+};
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif  // __FASTJET_CONTRIB_NJETTINESSPLUGIN_HH__
Index: external/fastjet/contribs/Nsubjettiness/Nsubjettiness.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Nsubjettiness.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/Nsubjettiness.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,50 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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 "Nsubjettiness.hh"
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib {
+
+//result returns tau_N with normalization dependent on what is specified in constructor
+double Nsubjettiness::result(const PseudoJet& jet) const {
+   std::vector<fastjet::PseudoJet> particles = jet.constituents();
+   return _njettinessFinder.getTau(_N, particles);
+}
+
+TauComponents Nsubjettiness::component_result(const PseudoJet& jet) const {
+   std::vector<fastjet::PseudoJet> particles = jet.constituents();
+   return _njettinessFinder.getTauComponents(_N, particles);
+}
+
+//ratio result uses Nsubjettiness result to find the ratio tau_N/tau_M, where N and M are specified by user
+double NsubjettinessRatio::result(const PseudoJet& jet) const {
+   double numerator = _nsub_numerator.result(jet);
+   double denominator = _nsub_denominator.result(jet);
+   return numerator/denominator;
+}
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/Nsubjettiness/Nsubjettiness.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Nsubjettiness.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/Nsubjettiness.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,137 @@
+// Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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_NSUBJETTINESS_HH__
+#define __FASTJET_CONTRIB_NSUBJETTINESS_HH__
+
+#include <fastjet/internal/base.hh>
+
+#include "Njettiness.hh"
+
+#include "fastjet/FunctionOfPseudoJet.hh"
+#include <string>
+#include <climits>
+
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib {
+
+//------------------------------------------------------------------------
+/// \class Nsubjettiness
+/// Nsubjettiness extends the concept of Njettiness to a jet shape, but other
+/// than the set of particles considered, they are identical.  This class
+/// wraps the core Njettiness code to provide the fastjet::FunctionOfPseudoJet
+/// interface for convenience in larger analyses.  See NjettinessPlugin.hh for
+/// definitions of tau_N and the constructor options.
+class Nsubjettiness : public FunctionOfPseudoJet<double> {
+
+public:
+
+   // Main constructor, which takes N, axes/measure modes,
+   // and up to four parameters for parameters (i.e. beta, Rcutoff, etc depending on measure)
+   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) {}
+
+   // Old constructor for backwards compatibility with v1.0,
+   // where normalized_cutoff_measure was the only option
+   Nsubjettiness(int N,
+                 Njettiness::AxesMode axes_mode,
+                 double beta,
+                 double R0,
+                 double Rcutoff=std::numeric_limits<double>::max())
+   : _njettinessFinder(axes_mode, Njettiness::normalized_cutoff_measure, beta, R0, Rcutoff), _N(N) {}
+
+
+   /// returns tau_N, measured on the constituents of this jet 
+   double result(const PseudoJet& jet) const;
+
+   /// returns components of tau_N, so that user can find individual tau values.
+   TauComponents component_result(const PseudoJet& jet) const;
+   
+   /// returns current axes found by result() calculation
+   std::vector<fastjet::PseudoJet> currentAxes() const {
+      return _njettinessFinder.currentAxes();
+   }
+
+   /// returns seed axes used for onepass minimization (otherwise same as currentAxes)
+   std::vector<fastjet::PseudoJet> seedAxes() const {
+      return _njettinessFinder.seedAxes();
+   }
+   
+   // To set axes for manual use
+   void setAxes(std::vector<fastjet::PseudoJet> myAxes) {
+      // Cross check that manual axes are being used is in Njettiness
+    	_njettinessFinder.setAxes(myAxes);
+   }
+   
+   
+private:
+   
+   mutable Njettiness _njettinessFinder; // TODO:  should muck with this so result can be const without this mutable
+   int _N;
+
+};
+
+
+//------------------------------------------------------------------------
+/// \class NsubjettinessRatio
+// NsubjettinessRatio uses the results from Nsubjettiness to calculate the ratio
+// tau_N/tau_M, where N and M are specified by the user. The ratio of different tau values
+// is often used in analyses, so this class is helpful to streamline code.
+class NsubjettinessRatio : public FunctionOfPseudoJet<double> {
+public:
+
+   // Main constructor.  Apart from specifying both N and M, the same options as Nsubjettiness
+   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 result(const PseudoJet& jet) const;
+
+private: 
+
+   Nsubjettiness _nsub_numerator;
+   Nsubjettiness _nsub_denominator;
+
+};
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif  // __FASTJET_CONTRIB_NSUBJETTINESS_HH__
Index: external/fastjet/contribs/Nsubjettiness/VERSION
===================================================================
--- external/fastjet/contribs/Nsubjettiness/VERSION	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/VERSION	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,1 @@
+1.1.0-beta4
Index: external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,69 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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 "WinnerTakeAllRecombiner.hh"
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+std::string WinnerTakeAllRecombiner::description() const {
+   return "Winner Take All scheme recombination";
+}
+
+// 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 
+// favors central jets
+void WinnerTakeAllRecombiner::recombine(const fastjet::PseudoJet & pa, const fastjet::PseudoJet & pb, fastjet::PseudoJet & pab) const {
+   double a_pt = pa.perp(), b_pt = pb.perp(), a_rap = pa.rap(), b_rap = pb.rap();
+   
+   // special case of alpha = 1, everything is just pt (made separate so that pow function isn't called)
+   if (_alpha == 1.0) {
+      if (a_pt >= b_pt) {
+         pab.reset_PtYPhiM(a_pt + b_pt, a_rap, pa.phi());
+      }
+      else if (b_pt > a_pt) {
+         pab.reset_PtYPhiM(a_pt + b_pt, b_rap, pb.phi());
+      }
+   }
+
+   // every other case uses additional cosh(rap) term
+   else {
+      double a_metric = a_pt*pow(cosh(a_rap), 1.0-_alpha);
+      double b_metric = b_pt*pow(cosh(b_rap), 1.0-_alpha);
+      if (a_metric >= b_metric) {
+   	  double new_pt = a_pt + b_pt*pow(cosh(b_rap)/cosh(a_rap), 1.0-_alpha);
+   	  pab.reset_PtYPhiM(new_pt, a_rap, pa.phi());
+      }
+      if (b_metric > a_metric) {
+   	  double new_pt = b_pt + a_pt*pow(cosh(a_rap)/cosh(b_rap), 1.0-_alpha);
+   	  pab.reset_PtYPhiM(new_pt, b_rap, pb.phi());
+      }
+   }
+}
+
+} //namespace contrib
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ external/fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -0,0 +1,62 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//----------------------------------------------------------------------
+// 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_WINNERTAKEALLRECOMBINER_HH__
+#define __FASTJET_CONTRIB_WINNERTAKEALLRECOMBINER_HH__
+
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/JetDefinition.hh"
+
+#include <cmath>
+#include <vector>
+#include <list>
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib {
+
+//------------------------------------------------------------------------
+/// \class WinnerTakeAllRecombiner
+// WinnerTakeAllRecombiner defines a new recombination scheme by inheriting from JetDefinition::Recombiner.
+// This scheme compares the energy of two input particles, and then combines them into a particle with
+// an energy equal to the sum of the two particle energies and a direction identical to that of the harder 
+// particle. This creates a jet with an axis guaranteed to align with a particle in the event.
+class WinnerTakeAllRecombiner : public fastjet::JetDefinition::Recombiner {
+public:
+	// Constructor to choose value of alpha (defaulted to 1 for normal pT sum)
+   WinnerTakeAllRecombiner(double alpha = 1.0) : _alpha(alpha) {}
+
+   virtual std::string description() const;
+   
+   /// recombine pa and pb and put result into pab
+   virtual void recombine(const fastjet::PseudoJet & pa, const fastjet::PseudoJet & pb, fastjet::PseudoJet & pab) const;
+
+private:
+	double _alpha; //power of (pt/E) term when recombining particles
+};
+
+} //namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif  // __FASTJET_CONTRIB_WINNERTAKEALLRECOMBINER_HH__
Index: modules/FastJetFinder.cc
===================================================================
--- modules/FastJetFinder.cc	(revision a0431dcbf4f6030aaef06db0f34e0959b4b9c6ae)
+++ modules/FastJetFinder.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -47,11 +47,18 @@
 #include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
 
+#include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh" 
+#include "fastjet/contribs/Nsubjettiness/Njettiness.hh"
+#include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
+#include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh"
+
 using namespace std;
 using namespace fastjet;
+using namespace fastjet::contrib;
+
 
 //------------------------------------------------------------------------------
 
 FastJetFinder::FastJetFinder() :
-  fPlugin(0), fDefinition(0), fAreaDefinition(0), fItInputArray(0)
+  fPlugin(0), fDefinition(0), fAreaDefinition(0), fItInputArray(0), fRecomb(0), fNjettinessPlugin(0)
 {
 
@@ -69,6 +76,9 @@
 void FastJetFinder::Init()
 {
+  
   JetDefinition::Plugin *plugin = NULL;
-
+  JetDefinition::Recombiner *recomb = NULL;
+  NjettinessPlugin *njet_plugin = NULL;
+  
   // read eta ranges
 
@@ -99,7 +109,16 @@
   fJetPTMin = GetDouble("JetPTMin", 10.0);
 
+  //-- N(sub)jettiness parameters --
+  
+  fComputeNsubjettiness = GetBool("ComputeNsubjettiness", false);
+  fBeta = GetDouble("Beta", 1.0);
+  fAxisMode = GetInt("AxisMode", 1);
+  fRcutOff = GetDouble("RcutOff", 0.8); //used only if Njettiness is used as jet clustering algo (case 8)
+  fN = GetInt("N", 2);                  //used only if Njettiness is used as jet clustering algo (case 8)
+  
   // ---  Jet Area Parameters ---
   fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
   fComputeRho = GetBool("ComputeRho", false);
+  
   // - ghost based areas -
   fGhostEtaMax = GetDouble("GhostEtaMax", 5.0);
@@ -109,4 +128,5 @@
   fPtScatter = GetDouble("PtScatter", 0.1);
   fMeanGhostPt = GetDouble("MeanGhostPt", 1.0E-100);
+  
   // - voronoi based areas -
   fEffectiveRfact = GetDouble("EffectiveRfact", 1.0);
@@ -159,8 +179,19 @@
       fDefinition = new fastjet::JetDefinition(fastjet::antikt_algorithm, fParameterR);
       break;
-  }
+    case 7:
+      recomb = new fastjet::contrib::WinnerTakeAllRecombiner();
+      fDefinition = new fastjet::JetDefinition(fastjet::antikt_algorithm, fParameterR, recomb, Best);
+      break;
+    case 8:
+      njet_plugin = new fastjet::contrib::NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
+      fDefinition = new fastjet::JetDefinition(njet_plugin);
+      break;
+  }
+ 
 
   fPlugin = plugin;
-
+  fRecomb = recomb;
+  fNjettinessPlugin = njet_plugin;
+  
   ClusterSequence::print_banner();
 
@@ -184,4 +215,6 @@
   if(fAreaDefinition) delete fAreaDefinition;
   if(fPlugin) delete static_cast<JetDefinition::Plugin*>(fPlugin);
+  if(fRecomb) delete static_cast<JetDefinition::Recombiner*>(fRecomb);
+  if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin*>(fNjettinessPlugin);
 }
 
@@ -248,4 +281,5 @@
   outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
 
+
   // loop over all jets and export them
   detaMax = 0.0;
@@ -265,4 +299,5 @@
     inputList.clear();
     inputList = sequence->constituents(*itOutputList);
+      
     for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
     {
@@ -289,4 +324,31 @@
     candidate->DeltaPhi = dphiMax;
 
+    // --- compute N-subjettiness with N = 1,2,3,4,5 ----
+    
+    if(fComputeNsubjettiness)
+    {
+      Njettiness::AxesMode axisMode;
+     
+      if (fAxisMode == 1) axisMode = Njettiness::wta_kt_axes;
+      if (fAxisMode == 2) axisMode = Njettiness::onepass_wta_kt_axes;
+      if (fAxisMode == 3) axisMode = Njettiness::kt_axes;
+      if (fAxisMode == 4) axisMode = Njettiness::onepass_kt_axes;
+       
+      Njettiness::MeasureMode measureMode = Njettiness::unnormalized_measure;
+      
+      Nsubjettiness nSub1(1, axisMode, measureMode, fBeta);
+      Nsubjettiness nSub2(2, axisMode, measureMode, fBeta);
+      Nsubjettiness nSub3(3, axisMode, measureMode, fBeta); 
+      Nsubjettiness nSub4(4, axisMode, measureMode, fBeta); 
+      Nsubjettiness nSub5(5, axisMode, measureMode, fBeta); 
+    
+      candidate -> Tau1 = nSub1(*itOutputList);
+      candidate -> Tau2 = nSub2(*itOutputList);
+      candidate -> Tau3 = nSub3(*itOutputList);
+      candidate -> Tau4 = nSub4(*itOutputList);
+      candidate -> Tau5 = nSub5(*itOutputList);
+    }
+
+
     fOutputArray->Add(candidate);
   }
Index: modules/FastJetFinder.h
===================================================================
--- modules/FastJetFinder.h	(revision a0431dcbf4f6030aaef06db0f34e0959b4b9c6ae)
+++ modules/FastJetFinder.h	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -41,4 +41,7 @@
 
   void *fPlugin; //!
+  void *fRecomb; //!
+  void *fNjettinessPlugin; //!
+  
   fastjet::JetDefinition *fDefinition; //!
 
@@ -55,4 +58,12 @@
   Double_t fOverlapThreshold;
 
+  //-- N (sub)jettiness parameters --
+  
+  Bool_t fComputeNsubjettiness;
+  Double_t fBeta;
+  Int_t fAxisMode;
+  Double_t fRcutOff; 
+  Int_t fN ;
+  
   // --- FastJet Area method --------
 
Index: modules/TreeWriter.cc
===================================================================
--- modules/TreeWriter.cc	(revision a0431dcbf4f6030aaef06db0f34e0959b4b9c6ae)
+++ modules/TreeWriter.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
@@ -554,5 +554,13 @@
     entry->FracPt[3] = candidate->FracPt[3];
     entry->FracPt[4] = candidate->FracPt[4];
-
+   
+    //--- N-subjettiness variables ----
+    
+    entry->Tau1 = candidate->Tau1;
+    entry->Tau2 = candidate->Tau2;
+    entry->Tau3 = candidate->Tau3;
+    entry->Tau4 = candidate->Tau4;
+    entry->Tau5 = candidate->Tau5;
+   
     FillParticles(candidate, &entry->Particles);
   }
