Index: external/fastjet/contribs/Nsubjettiness/MeasureDefinition.hh
===================================================================
--- external/fastjet/contribs/Nsubjettiness/MeasureDefinition.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/contribs/Nsubjettiness/MeasureDefinition.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
@@ -0,0 +1,853 @@
+//  Nsubjettiness Package
+//  Questions/Comments?  jthaler@jthaler.net
+//
+//  Copyright (c) 2011-14
+//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
+//
+//  $Id: MeasureDefinition.hh 828 2015-07-20 14:52:06Z jthaler $
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+
+#ifndef __FASTJET_CONTRIB_MEASUREDEFINITION_HH__
+#define __FASTJET_CONTRIB_MEASUREDEFINITION_HH__
+
+#include "fastjet/PseudoJet.hh"
+#include <cmath>
+#include <vector>
+#include <list>
+#include <limits>
+
+#include "TauComponents.hh"
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+
+
+// The following Measures are available (and their relevant arguments):
+// Recommended for usage as jet shapes
+class DefaultMeasure;               // Default measure from which next classes derive (should not be called directly)
+class NormalizedMeasure;            // (beta,R0)
+class UnnormalizedMeasure;          // (beta)
+class NormalizedCutoffMeasure;      // (beta,R0,Rcutoff)
+class UnnormalizedCutoffMeasure;    // (beta,Rcutoff)
+
+// New measures as of v2.2
+// Recommended for usage as event shapes (or for jet finding)
+class ConicalMeasure;               // (beta,R)
+class OriginalGeometricMeasure;     // (R)
+class ModifiedGeometricMeasure;     // (R)
+class ConicalGeometricMeasure;      // (beta, gamma, R)
+class XConeMeasure;                 // (beta, R)
+   
+// Formerly GeometricMeasure, now no longer recommended, kept commented out only for cross-check purposes
+//class DeprecatedGeometricMeasure;         // (beta)
+//class DeprecatedGeometricCutoffMeasure;   // (beta,Rcutoff)
+
+   
+///////
+//
+// MeasureDefinition
+//
+///////
+
+//This is a helper class for the Minimum Axes Finders. It is defined later.
+class LightLikeAxis;                                          
+   
+///------------------------------------------------------------------------
+/// \class MeasureDefinition
+/// \brief Base class for measure definitions
+///
+/// This is the base class for measure definitions.  Derived classes will calculate
+/// the tau_N of a jet given a specific measure and a set of axes.  The measure is
+/// determined by various jet and beam distances (and possible normalization factors).
+///------------------------------------------------------------------------
+class MeasureDefinition {
+   
+public:
+   
+   /// Description of measure and parameters
+   virtual std::string description() const = 0;
+   
+   /// In derived classes, this should return a copy of the corresponding derived class
+   virtual MeasureDefinition* create() const = 0;
+
+   //The following five functions define the measure by which tau_N is calculated,
+   //and are overloaded by the various measures below
+   
+   /// Distanes to jet axis.  This is called many times, so needs to be as fast as possible
+   /// Unless overloaded, it just calls jet_numerator
+   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      return jet_numerator(particle,axis);
+   }
+
+   /// Distanes to beam.  This is called many times, so needs to be as fast as possible
+   /// Unless overloaded, it just calls beam_numerator
+   virtual double beam_distance_squared(const fastjet::PseudoJet& particle) const {
+      return beam_numerator(particle);
+   }
+   
+   /// The jet measure used in N-(sub)jettiness
+   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const = 0;
+   /// The beam measure used in N-(sub)jettiness
+   virtual double beam_numerator(const fastjet::PseudoJet& particle) const = 0;
+   
+   /// A possible normalization factor
+   virtual double denominator(const fastjet::PseudoJet& particle) const = 0;
+   
+   /// Run a one-pass minimization routine.  There is a generic one-pass minimization that works for a wide variety of measures.
+   /// This should be overloaded to create a measure-specific minimization scheme
+   virtual std::vector<fastjet::PseudoJet> get_one_pass_axes(int n_jets,
+                                                             const std::vector<fastjet::PseudoJet>& inputs,
+                                                             const std::vector<fastjet::PseudoJet>& seedAxes,
+                                                             int nAttempts = 1000,         // cap number of iterations
+                                                             double accuracy = 0.0001      // cap distance of closest approach
+   ) const;
+   
+public:
+   
+   /// Returns the tau value for a give set of particles and axes
+   double result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
+      return component_result(particles,axes).tau();
+   }
+   
+   /// Short-hand for the result() function
+   inline double operator() (const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const {
+      return result(particles,axes);
+   }
+   
+   /// Return all of the TauComponents for specific input particles and axes
+   TauComponents component_result(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;
+   
+   /// Create the partitioning according the jet/beam distances and stores them a TauPartition
+   TauPartition get_partition(const std::vector<fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes) const;
+
+   /// Calculate the tau result using an existing partition
+   TauComponents component_result_from_partition(const TauPartition& partition, const std::vector<fastjet::PseudoJet>& axes) const;
+
+   
+
+   /// virtual destructor
+   virtual ~MeasureDefinition(){}
+   
+protected:
+   
+   /// Flag set by derived classes to choose whether or not to use beam/denominator
+   TauMode _tau_mode;
+   
+   /// Flag set by derived classes to say whether cheap get_one_pass_axes method can be used (true by default)
+   bool _useAxisScaling;
+
+   /// This is the only constructor, which requires _tau_mode and _useAxisScaling to be manually set by derived classes.
+   MeasureDefinition() : _tau_mode(UNDEFINED_SHAPE), _useAxisScaling(true) {}
+
+   
+   /// Used by derived classes to set whether or not to use beam/denominator information
+   void setTauMode(TauMode tau_mode) {
+      _tau_mode = tau_mode;
+   }
+
+   /// Used by derived classes to say whether one can use cheap get_one_pass_axes
+   void setAxisScaling(bool useAxisScaling) {
+      _useAxisScaling = useAxisScaling;
+   }
+
+   /// Uses denominator information?
+   bool has_denominator() const { return (_tau_mode == NORMALIZED_JET_SHAPE || _tau_mode == NORMALIZED_EVENT_SHAPE);}
+   /// Uses beam information?
+   bool has_beam() const {return (_tau_mode == UNNORMALIZED_EVENT_SHAPE || _tau_mode == NORMALIZED_EVENT_SHAPE);}
+
+   /// Create light-like axis (used in default one-pass minimization routine)
+   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);
+   }
+   
+   /// Shorthand for squaring
+   static inline double sq(double x) {return x*x;}
+
+};
+   
+
+///////
+//
+// Default Measures
+//
+///////
+   
+   
+///------------------------------------------------------------------------
+/// \enum DefaultMeasureType
+/// \brief Options for default measure
+///
+/// Can be used to switch between pp and ee measure types in the DefaultMeasure
+///------------------------------------------------------------------------
+enum DefaultMeasureType {
+   pt_R,       ///  use transverse momenta and boost-invariant angles,
+   E_theta,    ///  use energies and angles,
+   lorentz_dot, ///  use dot product inspired measure
+   perp_lorentz_dot /// use conical geometric inspired measures
+};
+
+///------------------------------------------------------------------------
+/// \class DefaultMeasure
+/// \brief Base class for default N-subjettiness measure definitions
+///
+/// This class is the default measure as defined in the original N-subjettiness papers.
+/// Based on the conical measure, but with a normalization factor
+/// 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 DefaultMeasure : public MeasureDefinition {
+   
+public:
+
+   /// Description
+   virtual std::string description() const;
+   /// To allow copying around of these objects
+   virtual DefaultMeasure* create() const {return new DefaultMeasure(*this);}
+
+   /// fast jet distance
+   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      return angleSquared(particle, axis);
+   }
+   
+   /// fast beam distance
+   virtual double beam_distance_squared(const fastjet::PseudoJet& /*particle*/) const {
+      return _RcutoffSq;
+   }
+   
+   /// true jet distance (given by general definitions of energy and angle)
+   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const{
+      double jet_dist = angleSquared(particle, axis);
+      if (jet_dist > 0.0) {
+         return energy(particle) * std::pow(jet_dist,_beta/2.0);
+      } else {
+         return 0.0;
+      }
+   }
+   
+   /// true beam distance
+   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
+      return energy(particle) * std::pow(_Rcutoff,_beta);
+   }
+   
+   /// possible denominator for normalization
+   virtual double denominator(const fastjet::PseudoJet& particle) const {
+      return energy(particle) * std::pow(_R0,_beta);
+   }
+
+   /// Special minimization routine (from v1.0 of N-subjettiness)
+   virtual std::vector<fastjet::PseudoJet> get_one_pass_axes(int n_jets,
+                                                             const std::vector<fastjet::PseudoJet>& inputs,
+                                                             const std::vector<fastjet::PseudoJet>& seedAxes,
+                                                             int nAttempts,   // cap number of iterations
+                                                             double accuracy  // cap distance of closest approach
+                                                             ) const;
+   
+protected:
+   double _beta;      ///< Angular exponent
+   double _R0;        ///< Normalization factor
+   double _Rcutoff;   ///< Cutoff radius
+   double _RcutoffSq; ///< Cutoff radius squared
+   DefaultMeasureType _measure_type;  ///< Type of measure used (i.e. pp style or ee style)
+   
+   
+   /// Constructor is protected so that no one tries to call this directly.
+   DefaultMeasure(double beta, double R0, double Rcutoff, DefaultMeasureType measure_type = pt_R)
+   : MeasureDefinition(), _beta(beta), _R0(R0), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)), _measure_type(measure_type)
+   {
+      if (beta <= 0) throw Error("DefaultMeasure:  You must choose beta > 0.");
+      if (R0 <= 0) throw Error("DefaultMeasure:  You must choose R0 > 0.");
+      if (Rcutoff <= 0) throw Error("DefaultMeasure:  You must choose Rcutoff > 0.");
+   }
+   
+   /// Added set measure method in case it becomes useful later
+   void setDefaultMeasureType(DefaultMeasureType measure_type) {
+      _measure_type = measure_type;
+   }
+
+   /// Generalized energy value (determined by _measure_type)
+   double energy(const PseudoJet& jet) const;
+   /// Generalized angle value (determined by _measure_type)
+   double angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const;
+
+   /// Name of _measure_type, so description will include the measure type
+   std::string measure_type_name() const {
+      if (_measure_type == pt_R) return "pt_R";
+      else if (_measure_type == E_theta) return "E_theta";
+      else if (_measure_type == lorentz_dot) return "lorentz_dot";
+      else if (_measure_type == perp_lorentz_dot) return "perp_lorentz_dot";
+      else return "Measure Type Undefined";
+   }      
+
+   /// templated for speed (TODO: probably should remove, since not clear that there is a speed gain)
+   template <int N> std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
+                                                              const std::vector <fastjet::PseudoJet> & inputJets,
+                                                              double accuracy) const;
+   
+   /// called by get_one_pass_axes to update axes iteratively
+   std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
+                                         const std::vector <fastjet::PseudoJet> & inputJets,
+                                         double accuracy) const;
+};
+   
+
+///------------------------------------------------------------------------
+/// \class NormalizedCutoffMeasure
+/// \brief Dimensionless default measure, with radius cutoff
+///
+/// This measure is just a wrapper for DefaultMeasure
+///------------------------------------------------------------------------
+class NormalizedCutoffMeasure : public DefaultMeasure {
+
+public:
+
+   /// Constructor
+   NormalizedCutoffMeasure(double beta, double R0, double Rcutoff, DefaultMeasureType measure_type = pt_R) 
+   : DefaultMeasure(beta, R0, Rcutoff, measure_type) {
+      setTauMode(NORMALIZED_JET_SHAPE);
+   }
+
+   /// Description
+   virtual std::string description() const;
+
+   /// For copying purposes
+   virtual NormalizedCutoffMeasure* create() const {return new NormalizedCutoffMeasure(*this);}
+
+};
+
+///------------------------------------------------------------------------
+/// \class NormalizedMeasure
+/// \brief Dimensionless default measure, with no cutoff
+///
+/// This measure is the same as NormalizedCutoffMeasure, with Rcutoff taken to infinity.
+///------------------------------------------------------------------------
+class NormalizedMeasure : public NormalizedCutoffMeasure {
+
+public:
+
+   /// Constructor
+   NormalizedMeasure(double beta, double R0, DefaultMeasureType measure_type = pt_R)
+   : NormalizedCutoffMeasure(beta, R0, std::numeric_limits<double>::max(), measure_type) {
+      _RcutoffSq = std::numeric_limits<double>::max();
+      setTauMode(NORMALIZED_JET_SHAPE);
+   }
+   
+   /// Description
+   virtual std::string description() const;
+   /// For copying purposes
+   virtual NormalizedMeasure* create() const {return new NormalizedMeasure(*this);}
+
+};
+   
+
+///------------------------------------------------------------------------
+/// \class UnnormalizedCutoffMeasure
+/// \brief Dimensionful default measure, with radius cutoff
+///
+/// This class is the unnormalized conical measure. The only difference from NormalizedCutoffMeasure
+/// is that the denominator is defined to be 1.0 by setting _has_denominator to false.
+/// class UnnormalizedCutoffMeasure : public NormalizedCutoffMeasure {
+///------------------------------------------------------------------------
+class UnnormalizedCutoffMeasure : public DefaultMeasure {
+   
+public:
+   
+   /// Since all methods are identical, UnnormalizedMeasure inherits directly
+   /// from NormalizedMeasure. R0 is a dummy value since the value of R0 is unecessary for this class,
+   /// and the "false" flag sets _has_denominator in MeasureDefinition to false so no denominator is used.
+   UnnormalizedCutoffMeasure(double beta, double Rcutoff, DefaultMeasureType measure_type = pt_R)
+   : DefaultMeasure(beta, std::numeric_limits<double>::quiet_NaN(), Rcutoff, measure_type) {
+      setTauMode(UNNORMALIZED_EVENT_SHAPE);
+   }
+
+   /// Description
+   virtual std::string description() const;
+   /// For copying purposes
+   virtual UnnormalizedCutoffMeasure* create() const {return new UnnormalizedCutoffMeasure(*this);}
+
+};
+
+   
+///------------------------------------------------------------------------
+/// \class UnnormalizedMeasure
+/// \brief Dimensionless default measure, with no cutoff
+///
+/// This measure is the same as UnnormalizedCutoffMeasure, with Rcutoff taken to infinity.
+///------------------------------------------------------------------------
+class UnnormalizedMeasure : public UnnormalizedCutoffMeasure {
+   
+public:
+   /// Since all methods are identical, UnnormalizedMeasure inherits directly
+   /// from NormalizedMeasure. R0 is a dummy value since the value of R0 is unecessary for this class,
+   /// and the "false" flag sets _has_denominator in MeasureDefinition to false so no denominator is used.
+   UnnormalizedMeasure(double beta, DefaultMeasureType measure_type = pt_R)
+   : UnnormalizedCutoffMeasure(beta, std::numeric_limits<double>::max(), measure_type) {
+      _RcutoffSq = std::numeric_limits<double>::max();
+      setTauMode(UNNORMALIZED_JET_SHAPE);
+   }
+
+   /// Description
+   virtual std::string description() const;
+   
+   /// For copying purposes
+   virtual UnnormalizedMeasure* create() const {return new UnnormalizedMeasure(*this);}
+   
+};
+
+
+///------------------------------------------------------------------------
+/// \class ConicalMeasure
+/// \brief Dimensionful event-shape measure, with radius cutoff
+///
+/// Very similar to UnnormalizedCutoffMeasure, but with different normalization convention
+/// and using the new default one-pass minimization algorithm.
+/// Axes are also made to be light-like to ensure sensible behavior
+/// Intended to be used as an event shape.
+///------------------------------------------------------------------------
+class ConicalMeasure : public MeasureDefinition {
+   
+public:
+   
+   /// Constructor
+   ConicalMeasure(double beta, double Rcutoff)
+   :   MeasureDefinition(), _beta(beta), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)) {
+      if (beta <= 0) throw Error("ConicalMeasure:  You must choose beta > 0.");
+      if (Rcutoff <= 0) throw Error("ConicalMeasure:  You must choose Rcutoff > 0.");
+      setTauMode(UNNORMALIZED_EVENT_SHAPE);
+   }
+   
+   /// Description
+   virtual std::string description() const;
+   /// For copying purposes
+   virtual ConicalMeasure* create() const {return new ConicalMeasure(*this);}
+   
+   /// fast jet distance
+   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      PseudoJet lightAxis = lightFrom(axis);
+      return particle.squared_distance(lightAxis);
+   }
+   
+   /// fast beam distance
+   virtual double beam_distance_squared(const fastjet::PseudoJet&  /*particle*/) const {
+      return _RcutoffSq;
+   }
+   
+   
+   /// true jet distance
+   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      PseudoJet lightAxis = lightFrom(axis);
+      double jet_dist = particle.squared_distance(lightAxis)/_RcutoffSq;
+      double jet_perp = particle.perp();
+      
+      if (_beta == 2.0) {
+         return jet_perp * jet_dist;
+      } else {
+         return jet_perp * pow(jet_dist,_beta/2.0);
+      }
+   }
+   
+   /// true beam distance
+   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
+      return particle.perp();
+   }
+   
+   /// no denominator used for this measure
+   virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
+      return std::numeric_limits<double>::quiet_NaN();
+   }
+   
+protected:
+   double _beta;     ///< angular exponent
+   double _Rcutoff;  ///< effective jet radius
+   double _RcutoffSq;///< effective jet radius squared
+};
+   
+
+
+///------------------------------------------------------------------------
+/// \class OriginalGeometricMeasure
+/// \brief Dimensionful event-shape measure, with dot-product distances
+///
+/// This class is the original (and hopefully now correctly coded) geometric measure.
+/// This measure is defined by the Lorentz dot product between
+/// the particle and the axis.  This class does not include normalization of tau_N.
+/// New in Nsubjettiness v2.2
+/// NOTE: This is defined differently from the DeprecatedGeometricMeasure which are now commented out.
+///------------------------------------------------------------------------
+class OriginalGeometricMeasure : public MeasureDefinition {
+
+public:
+   /// Constructor
+   OriginalGeometricMeasure(double Rcutoff)
+   : MeasureDefinition(), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)) {
+      if (Rcutoff <= 0) throw Error("OriginalGeometricMeasure:  You must choose Rcutoff > 0.");
+      setTauMode(UNNORMALIZED_EVENT_SHAPE);
+      setAxisScaling(false);  // No need to rescale axes (for speed up in one-pass minimization)
+   }
+
+   /// Description
+   virtual std::string description() const;
+   /// For copying purposes
+   virtual OriginalGeometricMeasure* create() const {return new OriginalGeometricMeasure(*this);}
+   
+   // This class uses the default jet_distance_squared and beam_distance_squared
+
+   /// true jet measure
+   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      return dot_product(lightFrom(axis), particle)/_RcutoffSq;
+   }
+
+   /// true beam measure
+   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
+      fastjet::PseudoJet beam_a(0,0,1,1);
+      fastjet::PseudoJet beam_b(0,0,-1,1);
+      double min_perp = std::min(dot_product(beam_a, particle),dot_product(beam_b, particle));
+      return min_perp;
+   }
+
+   /// no denominator needed for this measure.
+   virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
+      return std::numeric_limits<double>::quiet_NaN();
+   }
+   
+protected:
+   double _Rcutoff;   ///< Effective jet radius (rho = R^2)
+   double _RcutoffSq; ///< Effective jet radius squared
+
+};
+
+
+///------------------------------------------------------------------------
+/// \class ModifiedGeometricMeasure
+/// \brief Dimensionful event-shape measure, with dot-product distances, modified beam measure
+///
+/// This class is the Modified geometric measure.  This jet measure is defined by the Lorentz dot product between
+/// the particle and the axis, as in the Original Geometric Measure. The beam measure is defined differently from
+/// the above OriginalGeometric to allow for more conical jets. New in Nsubjettiness v2.2
+///------------------------------------------------------------------------
+class ModifiedGeometricMeasure : public MeasureDefinition {
+
+public:
+   /// Constructor
+   ModifiedGeometricMeasure(double Rcutoff)
+   :  MeasureDefinition(), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)) {
+      if (Rcutoff <= 0) throw Error("ModifiedGeometricMeasure:  You must choose Rcutoff > 0.");
+      setTauMode(UNNORMALIZED_EVENT_SHAPE);
+      setAxisScaling(false);  // No need to rescale axes (for speed up in one-pass minimization)
+   }
+
+   /// Description
+   virtual std::string description() const;
+   /// For copying purposes
+   virtual ModifiedGeometricMeasure* create() const {return new ModifiedGeometricMeasure(*this);}
+   
+   // This class uses the default jet_distance_squared and beam_distance_squared
+
+   /// True jet measure
+   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      return dot_product(lightFrom(axis), particle)/_RcutoffSq;
+   }
+
+   /// True beam measure
+   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
+      fastjet::PseudoJet lightParticle = lightFrom(particle);
+      return 0.5*particle.mperp()*lightParticle.pt();
+   }
+
+   /// This measure does not require a denominator
+   virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
+      return std::numeric_limits<double>::quiet_NaN();
+   }
+   
+protected:
+   double _Rcutoff;   ///< Effective jet radius (rho = R^2)
+   double _RcutoffSq; ///< Effective jet radius squared
+
+
+};
+
+///------------------------------------------------------------------------
+/// \class ConicalGeometricMeasure
+/// \brief Dimensionful event-shape measure, basis for XCone jet algorithm
+///
+/// This class is the Conical Geometric measure.  This measure is defined by the Lorentz dot product between
+/// the particle and the axis normalized by the axis and particle pT, as well as a factor of cosh(y) to vary
+/// the rapidity depepdence of the beam. New in Nsubjettiness v2.2, and the basis for the XCone jet algorithm
+///------------------------------------------------------------------------
+class ConicalGeometricMeasure : public MeasureDefinition {
+
+public:
+   
+   /// Constructor
+   ConicalGeometricMeasure(double jet_beta, double beam_gamma, double Rcutoff)
+   :  MeasureDefinition(),  _jet_beta(jet_beta), _beam_gamma(beam_gamma), _Rcutoff(Rcutoff), _RcutoffSq(sq(Rcutoff)){
+      if (jet_beta <= 0)   throw Error("ConicalGeometricMeasure:  You must choose beta > 0.");
+      if (beam_gamma <= 0) throw Error("ConicalGeometricMeasure:  You must choose gamma > 0.");
+      if (Rcutoff <= 0)   throw Error("ConicalGeometricMeasure:  You must choose Rcutoff > 0.");
+      setTauMode(UNNORMALIZED_EVENT_SHAPE);
+   }
+
+   /// Description
+   virtual std::string description() const;
+   /// For copying purposes
+   virtual ConicalGeometricMeasure* create() const {return new ConicalGeometricMeasure(*this);}
+   
+   /// fast jet measure
+   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      fastjet::PseudoJet lightAxis = lightFrom(axis);
+      double pseudoRsquared = 2.0*dot_product(lightFrom(axis),particle)/(lightAxis.pt()*particle.pt());
+      return pseudoRsquared;
+   }
+
+   /// fast beam measure
+   virtual double beam_distance_squared(const fastjet::PseudoJet&  /*particle*/) const {
+      return _RcutoffSq;
+   }
+
+   /// true jet measure
+   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+      double jet_dist = jet_distance_squared(particle,axis)/_RcutoffSq;
+      if (jet_dist > 0.0) {
+         fastjet::PseudoJet lightParticle = lightFrom(particle);
+         double weight = (_beam_gamma == 1.0) ? 1.0 : std::pow(0.5*lightParticle.pt(),_beam_gamma - 1.0);
+         return particle.pt() * weight * std::pow(jet_dist,_jet_beta/2.0);
+      } else {
+         return 0.0;
+      }
+   }
+
+   /// true beam measure
+   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
+      fastjet::PseudoJet lightParticle = lightFrom(particle);
+      double weight = (_beam_gamma == 1.0) ? 1.0 : std::pow(0.5*lightParticle.pt(),_beam_gamma - 1.0);
+      return particle.pt() * weight;
+   }
+
+   /// no denominator needed
+   virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
+      return std::numeric_limits<double>::quiet_NaN();
+   }
+      
+protected:
+   double _jet_beta;   ///< jet angular exponent
+   double _beam_gamma; ///< beam angular exponent (gamma = 1.0 is recommended)
+   double _Rcutoff;    ///< effective jet radius
+   double _RcutoffSq;  ///< effective jet radius squared
+   
+};
+
+///------------------------------------------------------------------------
+/// \class XConeMeasure
+/// \brief Dimensionful event-shape measure used in XCone jet algorithm
+///
+/// This class is the XCone Measure.  This is the default measure for use with the
+/// XCone algorithm. It is identical to the conical geometric measure but with gamma = 1.0.
+///------------------------------------------------------------------------
+class XConeMeasure : public ConicalGeometricMeasure {
+
+public:
+   /// Constructor
+   XConeMeasure(double jet_beta, double R)
+   :   ConicalGeometricMeasure(jet_beta,
+                               1.0, // beam_gamma, hard coded at gamma = 1.0 default
+                               R    // Rcutoff scale
+                               ) { }
+
+   /// Description
+   virtual std::string description() const;
+   /// For copying purposes
+   virtual XConeMeasure* create() const {return new XConeMeasure(*this);}
+
+};
+
+///------------------------------------------------------------------------
+/// \class LightLikeAxis
+/// \brief Helper class to define light-like axes directions
+///
+/// This is a helper class for the minimization routines.
+/// It creates a convenient way of defining axes in order to better facilitate calculations.
+///------------------------------------------------------------------------
+class LightLikeAxis {
+
+public:
+   /// Bare constructor
+   LightLikeAxis() : _rap(0.0), _phi(0.0), _weight(0.0), _mom(0.0) {}
+   /// Constructor
+   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) {}
+   
+   /// Rapidity
+   double rap() const {return _rap;}
+   /// Azimuth
+   double phi() const {return _phi;}
+   /// weight factor
+   double weight() const {return _weight;}
+   /// pt momentum
+   double mom() const {return _mom;}
+   
+   /// set rapidity
+   void set_rap(double my_set_rap) {_rap = my_set_rap;}
+   /// set azimuth
+   void set_phi(double my_set_phi) {_phi = my_set_phi;}
+   /// set weight factor
+   void set_weight(double my_set_weight) {_weight = my_set_weight;}
+   /// set pt momentum
+   void set_mom(double my_set_mom) {_mom = my_set_mom;}
+   /// set all kinematics
+   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 version
+   fastjet::PseudoJet ConvertToPseudoJet();
+   
+   /// Squared distance to PseudoJet
+   double DistanceSq(const fastjet::PseudoJet& input) const {
+      return DistanceSq(input.rap(),input.phi());
+   }
+
+   /// Distance to PseudoJet
+   double Distance(const fastjet::PseudoJet& input) const {
+      return std::sqrt(DistanceSq(input));
+   }
+
+   /// Squared distance to Lightlikeaxis
+   double DistanceSq(const LightLikeAxis& input) const {
+      return DistanceSq(input.rap(),input.phi());
+   }
+
+   /// Distance to Lightlikeaxis
+   double Distance(const LightLikeAxis& input) const {
+      return std::sqrt(DistanceSq(input));
+   }
+
+private:
+   double _rap;    ///< rapidity
+   double _phi;    ///< azimuth
+   double _weight; ///< weight factor
+   double _mom;    ///< pt momentum
+   
+   /// Internal squared distance calculation
+   double DistanceSq(double rap2, double phi2) const {
+      double rap1 = _rap;
+      double phi1 = _phi;
+      
+      double distRap = rap1-rap2;
+      double distPhi = std::fabs(phi1-phi2);
+      if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
+      return distRap*distRap + distPhi*distPhi;
+   }
+   
+   /// Internal distance calculation
+   double Distance(double rap2, double phi2) const {
+      return std::sqrt(DistanceSq(rap2,phi2));
+   }
+      
+};
+   
+   
+////------------------------------------------------------------------------
+///// \class DeprecatedGeometricCutoffMeasure
+//// This class is the old, incorrectly coded geometric measure.
+//// It is kept in case anyone wants to check old code, but should not be used for production purposes.
+//class DeprecatedGeometricCutoffMeasure : public MeasureDefinition {
+//
+//public:
+//
+//   // Please, please don't use this.
+//   DeprecatedGeometricCutoffMeasure(double jet_beta, double Rcutoff)
+//   :   MeasureDefinition(),
+//      _jet_beta(jet_beta),
+//      _beam_beta(1.0), // This is hard coded, since alternative beta_beam values were never checked.
+//      _Rcutoff(Rcutoff),
+//      _RcutoffSq(sq(Rcutoff)) {
+//         setTauMode(UNNORMALIZED_EVENT_SHAPE);
+//         setAxisScaling(false);
+//         if (jet_beta != 2.0) {
+//         throw Error("Geometric minimization is currently only defined for beta = 2.0.");
+//      }
+//   }
+//
+//   virtual std::string description() const;
+//
+//   virtual DeprecatedGeometricCutoffMeasure* create() const {return new DeprecatedGeometricCutoffMeasure(*this);}
+//
+//   virtual double jet_distance_squared(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+//      fastjet::PseudoJet lightAxis = lightFrom(axis);
+//      double pseudoRsquared = 2.0*dot_product(lightFrom(axis),particle)/(lightAxis.pt()*particle.pt());
+//      return pseudoRsquared;
+//   }
+//
+//   virtual double beam_distance_squared(const fastjet::PseudoJet&  /*particle*/) const {
+//      return _RcutoffSq;
+//   }
+//
+//   virtual double jet_numerator(const fastjet::PseudoJet& particle, const fastjet::PseudoJet& axis) const {
+//      fastjet::PseudoJet lightAxis = lightFrom(axis);
+//      double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(lightAxis.pt(),_beam_beta - 1.0);
+//      return particle.pt() * weight * std::pow(jet_distance_squared(particle,axis),_jet_beta/2.0);
+//   }
+//
+//   virtual double beam_numerator(const fastjet::PseudoJet& particle) const {
+//      double weight = (_beam_beta == 1.0) ? 1.0 : std::pow(particle.pt()/particle.e(),_beam_beta - 1.0);
+//      return particle.pt() * weight * std::pow(_Rcutoff,_jet_beta);
+//   }
+//
+//   virtual double denominator(const fastjet::PseudoJet&  /*particle*/) const {
+//      return std::numeric_limits<double>::quiet_NaN();
+//   }
+//
+//   virtual std::vector<fastjet::PseudoJet> get_one_pass_axes(int n_jets,
+//                                                             const std::vector<fastjet::PseudoJet>& inputs,
+//                                                             const std::vector<fastjet::PseudoJet>& seedAxes,
+//                                                             int nAttempts,    // cap number of iterations
+//                                                             double accuracy   // cap distance of closest approach
+//                                                            ) const;
+//
+//protected:
+//   double _jet_beta;
+//   double _beam_beta;
+//   double _Rcutoff;
+//   double _RcutoffSq;
+//
+//};
+//
+//// ------------------------------------------------------------------------
+//// / \class DeprecatedGeometricMeasure
+//// Same as DeprecatedGeometricMeasureCutoffMeasure, but with Rcutoff taken to infinity.
+//// NOTE:  This class should not be used for production purposes.
+//class DeprecatedGeometricMeasure : public DeprecatedGeometricCutoffMeasure {
+//
+//public:
+//   DeprecatedGeometricMeasure(double beta)
+//   : DeprecatedGeometricCutoffMeasure(beta,std::numeric_limits<double>::max()) {
+//      _RcutoffSq = std::numeric_limits<double>::max();
+//      setTauMode(UNNORMALIZED_JET_SHAPE);
+//   }
+//
+//   virtual std::string description() const;
+//
+//   virtual DeprecatedGeometricMeasure* create() const {return new DeprecatedGeometricMeasure(*this);}
+//};
+
+
+} //namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif  // __FASTJET_CONTRIB_MEASUREDEFINITION_HH__
