Index: external/fastjet/CircularRange.hh
===================================================================
--- external/fastjet/CircularRange.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/CircularRange.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: CircularRange.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: CircularRange.hh 4074 2016-03-08 09:09:25Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -35,4 +35,5 @@
 #include "fastjet/RangeDefinition.hh"
 #include "fastjet/Error.hh"
+#include "fastjet/internal/deprecated.hh"
 
 // for backwards compatibility: one should now use SelectorCircle,
@@ -50,7 +51,9 @@
 public:
   /// constructor
+  FASTJET_DEPRECATED_MSG("CircularRange is deprecated since FastJet 3.0. Use SelectorCircle instead")
   CircularRange() {_set_invalid_rapphi();}
   
   /// initialise CircularRange with a jet
+  FASTJET_DEPRECATED_MSG("CircularRange is deprecated since FastJet 3.0. Use SelectorCircle instead")
   CircularRange(const fastjet::PseudoJet & jet, double distance) {
                 _distance = distance;
@@ -60,4 +63,5 @@
 
   /// initialise CircularRange with a (rap,phi) point
+  FASTJET_DEPRECATED_MSG("CircularRange is deprecated since FastJet 3.0. Use SelectorCircle instead")
   CircularRange(double rap, double phi, double distance) {
                 _distance = distance;
@@ -67,4 +71,5 @@
 
   /// initialise CircularRange with just the radius parameter
+  FASTJET_DEPRECATED_MSG("CircularRange is deprecated since FastJet 3.0. Use SelectorCircle instead")
   CircularRange(double distance) {
                 _set_invalid_rapphi();
@@ -76,5 +81,5 @@
   
   /// return description of range
-  virtual inline std::string description() const {
+  virtual inline std::string description() const FASTJET_OVERRIDE {
     std::ostringstream ostr;
     ostr << "CircularRange: within distance "<< _distance << " of given jet or point." ;
@@ -83,8 +88,8 @@
   /// returns true since this range is localizable (i.e. set_position
   /// does something meaningful)
-  virtual inline bool is_localizable() const { return true; }
+  virtual inline bool is_localizable() const FASTJET_OVERRIDE { return true; }
   
   /// return bool according to whether (rap,phi) is in range
-  virtual inline bool is_in_range(double rap, double phi) const {
+  virtual inline bool is_in_range(double rap, double phi) const FASTJET_OVERRIDE {
      if (! _rapphi_are_valid()) {
        throw Error("Circular range used without a center having being defined (use set_position())");
@@ -98,5 +103,5 @@
 
   /// return the minimal and maximal rapidity of this range
-  virtual inline void get_rap_limits(double & rapmin, double & rapmax) const {
+  virtual inline void get_rap_limits(double & rapmin, double & rapmax) const FASTJET_OVERRIDE {
      rapmin = _rapjet - _distance;
      rapmax = _rapjet + _distance; }
Index: external/fastjet/ClosestPair2D.cc
===================================================================
--- external/fastjet/ClosestPair2D.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClosestPair2D.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClosestPair2D.cc 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClosestPair2D.cc 4059 2016-03-03 20:49:48Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -145,5 +145,5 @@
 
     // and create the search tree
-    _trees[ishift] = auto_ptr<Tree>(new Tree(shuffles, max_size));
+    _trees[ishift] = SharedPtr<Tree>(new Tree(shuffles, max_size));
 
     // now we look for the closest-pair candidates on this tree
@@ -174,5 +174,5 @@
     mindists2[i] = _points[i].neighbour_dist2;}
   
-  _heap = auto_ptr<MinHeap>(new MinHeap(mindists2, max_size));
+  _heap = SharedPtr<MinHeap>(new MinHeap(mindists2, max_size));
 }
 
Index: external/fastjet/ClusterSequence.cc
===================================================================
--- external/fastjet/ClusterSequence.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequence.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequence.cc 3809 2015-02-20 13:05:13Z soyez $
+// $Id: ClusterSequence.cc 4154 2016-07-20 16:20:48Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -152,6 +152,6 @@
   // set the pointer in the wrapper to this object to NULL to say that
   // we're going out of scope
-  if (_structure_shared_ptr()){
-    ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr()); 
+  if (_structure_shared_ptr){
+    ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr.get()); 
     // normally the csi is purely internal so it really should not be
     // NULL i.e assert should be OK
@@ -724,6 +724,9 @@
   }
   
-  bool code_should_never_reach_here = false;
-  assert(code_should_never_reach_here); 
+  //bool code_should_never_reach_here = false;
+  //assert(code_should_never_reach_here);
+
+  assert(0 && "Code should never reach here");
+
   return N2MHTLazy9;
 
@@ -775,4 +778,13 @@
 
 
+ClusterSequence & ClusterSequence::operator=(const ClusterSequence & cs) {
+  // self assignment is trivial
+  if (&cs != this) {
+    _deletes_self_when_unused = false;
+    transfer_from_sequence(cs);
+  }
+  return *this;
+}
+
 //----------------------------------------------------------------------
 // transfer the sequence contained in other_seq into our own;
@@ -815,5 +827,5 @@
 
   // clean up existing structure
-  if (_structure_shared_ptr()) {
+  if (_structure_shared_ptr) {
     // If there are jets associated with an old version of the CS and
     // a new one, keeping track of when to delete the CS becomes more
@@ -823,5 +835,5 @@
     // anything that is currently associated with the cluster sequence
     // should be told that its cluster sequence no longer exists
-    ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr()); 
+    ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr.get()); 
     assert(csi != NULL);
     csi->set_associated_cs(NULL);
@@ -1455,5 +1467,6 @@
 // initialise the history in a standard way
 void ClusterSequence::_add_step_to_history (
-	       const int step_number, const int parent1, 
+               //NO_LONGER_USED: const int step_number,
+               const int parent1, 
 	       const int parent2, const int jetp_index,
 	       const double dij) {
@@ -1469,5 +1482,7 @@
 
   int local_step = _history.size()-1;
-  assert(local_step == step_number);
+  //#ifndef __NO_ASSERTS__
+  //assert(local_step == step_number);
+  //#endif
 
   // sanity check: make sure the particles have not already been recombined
@@ -1669,6 +1684,10 @@
   int hist_j = _jets[jet_j].cluster_hist_index();
 
-  _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
+  _add_step_to_history(min(hist_i, hist_j), max(hist_i,hist_j),
 		       newjet_k, dij);
+
+  //  _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
+  //		       newjet_k, dij);
+
 
 }
@@ -1680,10 +1699,13 @@
 void ClusterSequence::_do_iB_recombination_step(
 				  const int jet_i, const double diB) {
-  // get history index
-  int newstep_k = _history.size();
-
   // recombine the jet with the beam
-  _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
+  _add_step_to_history(_jets[jet_i].cluster_hist_index(),BeamJet,
 		       Invalid, diB);
+
+  // // get history index
+  // int newstep_k = _history.size();
+  // 
+  // _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
+  //         	       Invalid, diB);
 
 }
Index: external/fastjet/ClusterSequence.hh
===================================================================
--- external/fastjet/ClusterSequence.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequence.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: ClusterSequence.hh 3911 2015-07-02 12:09:58Z salam $
+// $Id: ClusterSequence.hh 4154 2016-07-20 16:20:48Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -49,4 +49,6 @@
 #include "fastjet/ClusterSequenceStructure.hh"
 
+#include "fastjet/internal/deprecated.hh"
+
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
@@ -80,4 +82,7 @@
   }
 
+  /// explicit assignment operator for a ClusterSequence
+  ClusterSequence & operator=(const ClusterSequence & cs);
+  
   // virtual ClusterSequence destructor, in case any derived class
   // thinks of needing a destructor at some point
@@ -362,7 +367,10 @@
   /// As of FJ v3.1, this is deprecated, in line with the deprecation
   /// of auto_ptr in C++11
-  inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
+#ifdef FASTJET_HAVE_AUTO_PTR_INTERFACE
+  FASTJET_DEPRECATED_MSG("Please use ClusterSequence::plugin_associate_extras(Extras * extras_in)) instead")
+  inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in){
     _extras.reset(extras_in.release());
   }
+#endif
 
   /// returns true when the plugin is allowed to run the show.
@@ -730,7 +738,8 @@
   void _fast_NsqrtN_cluster();
 
-  void _add_step_to_history(const int step_number, const int parent1, 
-			       const int parent2, const int jetp_index,
-			       const double dij);
+  void _add_step_to_history( //const int step_number,
+                            const int parent1, 
+			    const int parent2, const int jetp_index,
+			    const double dij);
 
   /// internal routine associated with the construction of the unique
@@ -1037,7 +1046,15 @@
 template <class J> inline double ClusterSequence::_bj_dist(
                 const J * const jetA, const J * const jetB) const {
+  //#define FASTJET_NEW_DELTA_PHI
+#ifndef FASTJET_NEW_DELTA_PHI
+  //GPS+MC old version of Delta phi calculation
   double dphi = std::abs(jetA->phi - jetB->phi);
   double deta = (jetA->eta - jetB->eta);
   if (dphi > pi) {dphi = twopi - dphi;}
+#else 
+  //GPS+MC testing for 2015-02-faster-deltaR2
+  double dphi = pi-std::abs(pi-std::abs(jetA->phi - jetB->phi));
+  double deta = (jetA->eta - jetB->eta);
+#endif 
   return dphi*dphi + deta*deta;
 }
Index: external/fastjet/ClusterSequence1GhostPassiveArea.hh
===================================================================
--- external/fastjet/ClusterSequence1GhostPassiveArea.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequence1GhostPassiveArea.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequence1GhostPassiveArea.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClusterSequence1GhostPassiveArea.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -70,5 +70,5 @@
   /// AreaBase one rather than the ActiveArea one (which for which we
   /// do not have the information).
-  virtual double n_empty_jets(const Selector & selector) const {
+  virtual double n_empty_jets(const Selector & selector) const FASTJET_OVERRIDE {
     return ClusterSequenceAreaBase::n_empty_jets(selector);
   }
Index: external/fastjet/ClusterSequenceActiveArea.hh
===================================================================
--- external/fastjet/ClusterSequenceActiveArea.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequenceActiveArea.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequenceActiveArea.hh 3619 2014-08-13 14:17:19Z salam $
+// $Id: ClusterSequenceActiveArea.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -74,10 +74,10 @@
 	  const bool & writeout_combinations = false) ;
 
-  virtual double area (const PseudoJet & jet) const {
+  virtual double area (const PseudoJet & jet) const FASTJET_OVERRIDE {
                              return _average_area[jet.cluster_hist_index()];};
-  virtual double area_error (const PseudoJet & jet) const {
+  virtual double area_error (const PseudoJet & jet) const FASTJET_OVERRIDE {
                              return _average_area2[jet.cluster_hist_index()];};
 
-  virtual PseudoJet area_4vector (const PseudoJet & jet) const {
+  virtual PseudoJet area_4vector (const PseudoJet & jet) const FASTJET_OVERRIDE {
                     return _average_area_4vector[jet.cluster_hist_index()];};
 
@@ -110,9 +110,9 @@
   ///
   /// The selector passed as an argument needs to apply jet by jet.
-  virtual double empty_area(const Selector & selector) const;
+  virtual double empty_area(const Selector & selector) const FASTJET_OVERRIDE;
 
   /// return the true number of empty jets (replaces
   /// ClusterSequenceAreaBase::n_empty_jets(...))
-  virtual double n_empty_jets(const Selector & selector) const;
+  virtual double n_empty_jets(const Selector & selector) const FASTJET_OVERRIDE;
 
 protected:
Index: external/fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh
===================================================================
--- external/fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequenceActiveAreaExplicitGhosts.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClusterSequenceActiveAreaExplicitGhosts.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -97,5 +97,5 @@
 
   /// returns the area of a jet
-  virtual double area (const PseudoJet & jet) const;
+  virtual double area (const PseudoJet & jet) const FASTJET_OVERRIDE;
 
   /// returns a four vector corresponding to the sum (E-scheme) of the
@@ -103,8 +103,8 @@
   /// for a small contiguous area the p_t of the extended_area jet is
   /// equal to area of the jet.
-  virtual PseudoJet area_4vector (const PseudoJet & jet) const;
+  virtual PseudoJet area_4vector (const PseudoJet & jet) const FASTJET_OVERRIDE;
 
   /// true if a jet is made exclusively of ghosts
-  virtual bool is_pure_ghost(const PseudoJet & jet) const;
+  virtual bool is_pure_ghost(const PseudoJet & jet) const FASTJET_OVERRIDE;
 
   /// true if the entry in the history index corresponds to a
@@ -114,5 +114,5 @@
 
   /// this class does have explicit ghosts
-  virtual bool has_explicit_ghosts() const {return true;}
+  virtual bool has_explicit_ghosts() const FASTJET_OVERRIDE {return true;}
 
   /// return the total area, corresponding to a given Selector, that
@@ -120,5 +120,5 @@
   ///
   /// The selector needs to apply jet by jet
-  virtual double empty_area(const Selector & selector) const;
+  virtual double empty_area(const Selector & selector) const FASTJET_OVERRIDE;
 
   /// returns the total area under study
Index: external/fastjet/ClusterSequenceArea.hh
===================================================================
--- external/fastjet/ClusterSequenceArea.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequenceArea.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequenceArea.hh 3484 2014-07-29 21:39:39Z soyez $
+// $Id: ClusterSequenceArea.hh 4098 2016-03-15 16:38:22Z salam $
 //
 // Copyright (c) 2006-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -80,14 +80,14 @@
 
   /// return the area associated with the given jet
-  virtual double area       (const PseudoJet & jet) const {
+  virtual double area       (const PseudoJet & jet) const FASTJET_OVERRIDE {
     return _area_base->area(jet);}
 
   /// return the error (uncertainty) associated with the determination
   /// of the area of this jet
-  virtual double area_error (const PseudoJet & jet) const {
+  virtual double area_error (const PseudoJet & jet) const FASTJET_OVERRIDE {
     return _area_base->area_error(jet);}
 
   /// return the 4-vector area
-  virtual PseudoJet area_4vector(const PseudoJet & jet) const {
+  virtual PseudoJet area_4vector(const PseudoJet & jet) const FASTJET_OVERRIDE {
     return _area_base->area_4vector(jet);}
 
@@ -110,5 +110,5 @@
   /// jet (see the BackgroundEstimator and Subtractor tools for more
   /// advanced usage)
-  virtual double empty_area(const Selector & selector) const {
+  virtual double empty_area(const Selector & selector) const FASTJET_OVERRIDE {
     return _area_base->empty_area(selector);}
 
@@ -122,15 +122,15 @@
   /// jet (see the BackgroundEstimator and Subtractor tools for more
   /// advanced usage)
-  virtual double n_empty_jets(const Selector & selector) const {
+  virtual double n_empty_jets(const Selector & selector) const FASTJET_OVERRIDE {
     return _area_base->n_empty_jets(selector);
   }
 
   /// true if a jet is made exclusively of ghosts
-  virtual bool is_pure_ghost(const PseudoJet & jet) const {
+  virtual bool is_pure_ghost(const PseudoJet & jet) const FASTJET_OVERRIDE {
     return _area_base->is_pure_ghost(jet);
   }
 
   /// true if this ClusterSequence has explicit ghosts
-  virtual bool has_explicit_ghosts() const {
+  virtual bool has_explicit_ghosts() const FASTJET_OVERRIDE {
     return _area_base->has_explicit_ghosts();
   }
@@ -144,4 +144,5 @@
   /// jet (see the BackgroundEstimator and Subtractor tools for more
   /// advanced usage)
+  //FASTJET_DEPRECATED_MSG("ClusterSequenceArea::get_median_rho_and_sigma(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
   virtual void get_median_rho_and_sigma(const std::vector<PseudoJet> & all_jets,
 					const Selector & selector, 
@@ -149,7 +150,7 @@
                                         double & median, double & sigma,
                                         double & mean_area,
-					bool all_are_incl = false) const {
+					bool all_are_incl = false) const FASTJET_OVERRIDE {
     _warn_if_range_unsuitable(selector);
-    ClusterSequenceAreaBase::get_median_rho_and_sigma(
+    ClusterSequenceAreaBase::_get_median_rho_and_sigma(
                                  all_jets, selector, use_area_4vector,
 				 median, sigma, mean_area, all_are_incl);
@@ -160,9 +161,10 @@
   /// since we've overridden the 5-argument version above, we have to
   /// override the 4-argument version too.
+  //FASTJET_DEPRECATED_MSG("ClusterSequenceArea::get_median_rho_and_sigma(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
   virtual void get_median_rho_and_sigma(const Selector & selector, 
                                         bool use_area_4vector,
-                                        double & median, double & sigma) const {
-    ClusterSequenceAreaBase::get_median_rho_and_sigma(selector,use_area_4vector,
-                                                      median,sigma);
+                                        double & median, double & sigma) const FASTJET_OVERRIDE {
+    ClusterSequenceAreaBase::_get_median_rho_and_sigma(selector,use_area_4vector,
+                                                       median,sigma);
   }
 
@@ -171,10 +173,11 @@
   /// since we've overridden the multi-argument version above, we have to
   /// override the 5-argument version too.
+  //FASTJET_DEPRECATED_MSG("ClusterSequenceArea::get_median_rho_and_sigma(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
   virtual void get_median_rho_and_sigma(const Selector & selector, 
                                         bool use_area_4vector,
                                         double & median, double & sigma,
-					double & mean_area) const {
-    ClusterSequenceAreaBase::get_median_rho_and_sigma(selector,use_area_4vector,
-                                                      median,sigma, mean_area);
+					double & mean_area) const FASTJET_OVERRIDE {
+    ClusterSequenceAreaBase::_get_median_rho_and_sigma(selector,use_area_4vector,
+                                                       median,sigma, mean_area);
   }
 
@@ -183,11 +186,10 @@
   /// additionally checks compatibility between "range" and region in which
   /// ghosts are thrown.
+  //FASTJET_DEPRECATED_MSG("ClusterSequenceArea::parabolic_pt_per_unit_area(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")  
   virtual void parabolic_pt_per_unit_area(double & a, double & b, 
                                           const Selector & selector, 
                                           double exclude_above=-1.0, 
-                                          bool use_area_4vector=false) const {
-    _warn_if_range_unsuitable(selector);
-    ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
-                                a,b,selector, exclude_above, use_area_4vector);
+                                          bool use_area_4vector=false) const FASTJET_OVERRIDE {
+    return _parabolic_pt_per_unit_area(a,b,selector,exclude_above,use_area_4vector);
   }
 
@@ -204,8 +206,23 @@
                                  const JetDefinition & jet_def);
 
-  std::auto_ptr<ClusterSequenceAreaBase> _area_base;
+  SharedPtr<ClusterSequenceAreaBase> _area_base;
   AreaDefinition _area_def;
   static LimitedWarning _range_warnings;
   static LimitedWarning _explicit_ghosts_repeats_warnings;
+
+  // the following set of private methods are all deprecated. Their
+  // role is simply to hide the corresponding methods (without the
+  // first underscore) from the public interface so that they can be
+  // used internally until all the deprecated methods are removed.
+  // DO NOT USE ANY OF THESE METHODS: THEY ARE DEPRECATED AND WILL BE
+  // REMOVED.
+  virtual void _parabolic_pt_per_unit_area(double & a, double & b, 
+                                          const Selector & selector, 
+                                          double exclude_above=-1.0, 
+                                          bool use_area_4vector=false) const FASTJET_OVERRIDE {
+    _warn_if_range_unsuitable(selector);
+    ClusterSequenceAreaBase::_parabolic_pt_per_unit_area(
+                                a,b,selector, exclude_above, use_area_4vector);
+  }
 
 };
@@ -254,5 +271,5 @@
   }
   // now copy across the information from the area base class
-  _area_base = std::auto_ptr<ClusterSequenceAreaBase>(_area_base_ptr);
+  _area_base = SharedPtr<ClusterSequenceAreaBase>(_area_base_ptr);
   transfer_from_sequence(*_area_base);
 }
Index: external/fastjet/ClusterSequenceAreaBase.cc
===================================================================
--- external/fastjet/ClusterSequenceAreaBase.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequenceAreaBase.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,5 +1,5 @@
 
 //FJSTARTHEADER
-// $Id: ClusterSequenceAreaBase.cc 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClusterSequenceAreaBase.cc 4079 2016-03-09 12:20:03Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -79,32 +79,68 @@
 }
 
+// this is deprecated but used by other deprecated methods. So we hide
+// the implementation in a protected method so that (i) it can still
+// be used internally (without generating a compile-time warning when
+// building FastJet) and the interface can be marked as deprecated.
+// This can disappear once all the public interfaces have disappeared.
 double ClusterSequenceAreaBase::median_pt_per_unit_area(const Selector & selector) const {
-  return median_pt_per_unit_something(selector,false);
-}
-
+  return _median_pt_per_unit_area(selector);
+}
+
+// the hidden implementation
+double ClusterSequenceAreaBase::_median_pt_per_unit_area(const Selector & selector) const {
+  return _median_pt_per_unit_something(selector,false);
+}
+
+
+
+// this is deprecated but used by other deprecated methods. So we hide
+// the implementation in a protected method so that (i) it can still
+// be used internally (without generating a compile-time warning when
+// building FastJet) and the interface can be marked as deprecated.
+// This can disappear once all the public interfaces have disappeared.
 double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(const Selector & selector) const {
-  return median_pt_per_unit_something(selector,true);
-}
-
-
-//----------------------------------------------------------------------
-/// the median of (pt/area) for jets contained within range, counting
-/// the empty area as if it were made up of a collection of empty
-/// jets each of area (0.55 * pi R^2).
+  return _median_pt_per_unit_area_4vector(selector);
+}
+
+// the deprecated interface
+double ClusterSequenceAreaBase::_median_pt_per_unit_area_4vector(const Selector & selector) const {
+  return _median_pt_per_unit_something(selector,true);
+}
+
+
+//----------------------------------------------------------------------
+// this is deprecated but used by other deprecated methods. So we hide
+// the implementation in a protected method so that (i) it can still
+// be used internally (without generating a compile-time warning when
+// building FastJet) and the interface can be marked as deprecated.
+// This can disappear once all the public interfaces have disappeared.
 double ClusterSequenceAreaBase::median_pt_per_unit_something(
                 const Selector & selector, bool use_area_4vector) const {
-
+  return _median_pt_per_unit_something(selector, use_area_4vector);
+}
+
+// the median of (pt/area) for jets contained within range, counting
+// the empty area as if it were made up of a collection of empty
+// jets each of area (0.55 * pi R^2).
+double ClusterSequenceAreaBase::_median_pt_per_unit_something(
+                const Selector & selector, bool use_area_4vector) const {
   double median, sigma, mean_area;
-  get_median_rho_and_sigma(selector, use_area_4vector, median, sigma, mean_area);
+  _get_median_rho_and_sigma(selector, use_area_4vector, median, sigma, mean_area);
   return median;
-
-}
-
-
-//----------------------------------------------------------------------
-/// fits a form pt_per_unit_area(y) = a + b*y^2 for jets in range. 
-/// exclude_above allows one to exclude large values of pt/area from fit. 
-/// use_area_4vector = true uses the 4vector areas.
+}
+
+
+//----------------------------------------------------------------------
+/// fits a form pt_per_unit_area(y) = a + b*y^2 for jets in range.
+/// exclude_above allows one to exclude large values of pt/area from
+/// fit.  use_area_4vector = true uses the 4vector areas.
 void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
+       double & a, double & b, const Selector & selector, 
+       double exclude_above, bool use_area_4vector) const {
+  return _parabolic_pt_per_unit_area(a, b, selector, exclude_above, use_area_4vector);
+}
+
+void ClusterSequenceAreaBase::_parabolic_pt_per_unit_area(
        double & a, double & b, const Selector & selector, 
        double exclude_above, bool use_area_4vector) const {
@@ -161,16 +197,30 @@
 
 
-
+//----------------------------------------------------------------------
 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
             const Selector & selector, bool use_area_4vector,
             double & median, double & sigma, double & mean_area) const {
+  _get_median_rho_and_sigma(selector, use_area_4vector, median, sigma, mean_area);
+}
+
+void ClusterSequenceAreaBase::_get_median_rho_and_sigma(
+            const Selector & selector, bool use_area_4vector,
+            double & median, double & sigma, double & mean_area) const {
 
   vector<PseudoJet> incl_jets = inclusive_jets();
-  get_median_rho_and_sigma(incl_jets, selector, use_area_4vector,
-			   median, sigma, mean_area, true);
-}
-
+  _get_median_rho_and_sigma(incl_jets, selector, use_area_4vector,
+                            median, sigma, mean_area, true);
+}
 
 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
+            const vector<PseudoJet> & all_jets,
+            const Selector & selector, bool use_area_4vector,
+            double & median, double & sigma, double & mean_area,
+	    bool all_are_incl) const {
+  _get_median_rho_and_sigma(all_jets, selector, use_area_4vector,
+                            median, sigma, mean_area, all_are_incl);
+}
+
+void ClusterSequenceAreaBase::_get_median_rho_and_sigma(
             const vector<PseudoJet> & all_jets,
             const Selector & selector, bool use_area_4vector,
@@ -288,15 +338,21 @@
 
 
-/// return a vector of all subtracted jets, using area_4vector, given rho.
-/// Only inclusive_jets above ptmin are subtracted and returned.
-/// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
-/// i.e. not necessarily ordered in pt once subtracted
+// return a vector of all subtracted jets, using area_4vector, given rho.
+// Only inclusive_jets above ptmin are subtracted and returned.
+// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
+// i.e. not necessarily ordered in pt once subtracted
 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(const double rho,
                                                            const double ptmin) 
                                                            const {
+  return _subtracted_jets(rho,ptmin);
+}
+
+vector<PseudoJet> ClusterSequenceAreaBase::_subtracted_jets(const double rho,
+                                                            const double ptmin) 
+  const {
   vector<PseudoJet> sub_jets;
   vector<PseudoJet> jets_local = sorted_by_pt(inclusive_jets(ptmin));
   for (unsigned i=0; i<jets_local.size(); i++) {
-     PseudoJet sub_jet = subtracted_jet(jets_local[i],rho);
+     PseudoJet sub_jet = _subtracted_jet(jets_local[i],rho);
      sub_jets.push_back(sub_jet);
   }
@@ -304,14 +360,14 @@
 }
 
-/// return a vector of subtracted jets, using area_4vector.
-/// Only inclusive_jets above ptmin are subtracted and returned.
-/// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
-/// i.e. not necessarily ordered in pt once subtracted
+// return a vector of subtracted jets, using area_4vector.
+// Only inclusive_jets above ptmin are subtracted and returned.
+// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
+// i.e. not necessarily ordered in pt once subtracted
 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
                                                  const Selector & selector, 
 						 const double ptmin)
 						 const {
-  double rho = median_pt_per_unit_area_4vector(selector);
-  return subtracted_jets(rho,ptmin);
+  double rho = _median_pt_per_unit_area_4vector(selector);
+  return _subtracted_jets(rho,ptmin);
 }
 
@@ -320,4 +376,9 @@
 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
                                                   const double rho) const {
+  return _subtracted_jet(jet, rho);
+}
+
+PseudoJet ClusterSequenceAreaBase::_subtracted_jet(const PseudoJet & jet,
+                                                   const double rho) const {
   PseudoJet area4vect = area_4vector(jet);
   PseudoJet sub_jet;
@@ -343,6 +404,11 @@
 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
                                        const Selector & selector) const {
-  double rho = median_pt_per_unit_area_4vector(selector);
-  PseudoJet sub_jet = subtracted_jet(jet, rho);
+  return _subtracted_jet(jet, selector);
+}
+
+PseudoJet ClusterSequenceAreaBase::_subtracted_jet(const PseudoJet & jet,
+                                       const Selector & selector) const {
+  double rho = _median_pt_per_unit_area_4vector(selector);
+  PseudoJet sub_jet = _subtracted_jet(jet, rho);
   return sub_jet;
 }
@@ -353,6 +419,12 @@
                                               const double rho,
                                               bool use_area_4vector) const {
+  return _subtracted_pt(jet, rho, use_area_4vector);
+}
+
+double ClusterSequenceAreaBase::_subtracted_pt(const PseudoJet & jet,
+                                              const double rho,
+                                              bool use_area_4vector) const {
   if ( use_area_4vector ) { 
-     PseudoJet sub_jet = subtracted_jet(jet,rho);
+     PseudoJet sub_jet = _subtracted_jet(jet,rho);
      return sub_jet.perp();
   } else {
@@ -369,9 +441,9 @@
                                               bool use_area_4vector) const {
   if ( use_area_4vector ) { 
-     PseudoJet sub_jet = subtracted_jet(jet,selector);
+     PseudoJet sub_jet = _subtracted_jet(jet,selector);
      return sub_jet.perp();
   } else {
-     double rho = median_pt_per_unit_area(selector);
-     return subtracted_pt(jet,rho,false);
+     double rho = _median_pt_per_unit_area(selector);
+     return _subtracted_pt(jet,rho,false);
   }
 }  
Index: external/fastjet/ClusterSequenceAreaBase.hh
===================================================================
--- external/fastjet/ClusterSequenceAreaBase.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequenceAreaBase.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequenceAreaBase.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClusterSequenceAreaBase.hh 4098 2016-03-15 16:38:22Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -35,4 +35,5 @@
 #include "fastjet/LimitedWarning.hh"
 #include "fastjet/Selector.hh"
+#include "fastjet/internal/deprecated.hh"
 
 FASTJET_BEGIN_NAMESPACE
@@ -144,4 +145,5 @@
   /// apply jet-by-jet (see the BackgroundEstimator and Subtractor
   /// tools for more generic usages)
+  FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::median_pt_per_unit_area(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
   double median_pt_per_unit_area(const Selector & selector) const;
 
@@ -151,4 +153,5 @@
   /// The selector passed as an argument has to have a finite area and
   /// apply jet-by-jet
+  FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
   double median_pt_per_unit_area_4vector(const Selector & selector) const;
   
@@ -157,6 +160,7 @@
   /// - something_is_area_4vect = false -> use plain area
   /// - something_is_area_4vect = true  -> use 4-vector area
-  double median_pt_per_unit_something(
-                    const Selector & selector, bool use_area_4vector) const;
+  FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::median_pt_per_unit_something(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
+  double median_pt_per_unit_something(const Selector & selector,
+                                      bool use_area_4vector) const;
 
   /// using jets withing the selector range (and with 4-vector areas if
@@ -184,4 +188,5 @@
   /// using a scalar area causes one to neglect terms of relative
   /// order $R^2/8$ in the jet $p_t$.
+  //FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::get_median_rho_and_sigma(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
   virtual void get_median_rho_and_sigma(const Selector & selector, 
                                         bool use_area_4vector,
@@ -208,4 +213,5 @@
   /// apply jet-by-jet (see the BackgroundEstimator and Subtractor
   /// tools for more generic usages)
+  //FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::get_median_rho_and_sigma(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
   virtual void get_median_rho_and_sigma(const std::vector<PseudoJet> & all_jets,
 					const Selector & selector, 
@@ -221,10 +227,9 @@
   /// apply jet-by-jet (see the BackgroundEstimator and Subtractor
   /// tools for more generic usages)
+  //FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::get_median_rho_and_sigma(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
   virtual void get_median_rho_and_sigma(const Selector & selector, 
-                                bool use_area_4vector,
-                                double & median, double & sigma) const {
-    double mean_area;
-    get_median_rho_and_sigma(selector,  use_area_4vector,
-                             median,  sigma, mean_area);
+                                        bool use_area_4vector,
+                                        double & median, double & sigma) const{
+    return _get_median_rho_and_sigma(selector, use_area_4vector, median, sigma);
   }
   
@@ -238,4 +243,5 @@
   /// apply jet-by-jet (see the BackgroundEstimator and Subtractor
   /// tools for more generic usages)
+  //FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::parabolic_pt_per_unit_area(...) is deprecated since FastJet 3.0. Use the BackgroundEstimator series of tools instead")
   virtual void parabolic_pt_per_unit_area(double & a, double & b, 
                                           const Selector & selector, 
@@ -247,4 +253,5 @@
   /// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
   /// i.e. not necessarily ordered in pt once subtracted
+  FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::subtracted_jets(...) is deprecated since FastJet 3.0. Use the Subtractor tool (with the BackgroundEstimator series of tools) instead")
   std::vector<PseudoJet> subtracted_jets(const double rho,
                                          const double ptmin=0.0) const;
@@ -258,8 +265,10 @@
   /// apply jet-by-jet (see the BackgroundEstimator and Subtractor
   /// tools for more generic usages)
+  FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::subtracted_jets(...) is deprecated since FastJet 3.0. Use the Subtractor tool (with the BackgroundEstimator series of tools) instead")
   std::vector<PseudoJet> subtracted_jets(const Selector & selector, 
                                          const double ptmin=0.0) const;
 
   /// return a subtracted jet, using area_4vector, given rho
+  FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::subtracted_jet(...) is deprecated since FastJet 3.0. Use the Subtractor tool (with the BackgroundEstimator series of tools) instead")
   PseudoJet subtracted_jet(const PseudoJet & jet,
                            const double rho) const;
@@ -273,8 +282,10 @@
   /// apply jet-by-jet (see the BackgroundEstimator and Subtractor
   /// tools for more generic usages)
+  FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::subtracted_jet(...) is deprecated since FastJet 3.0. Use the Subtractor tool (with the BackgroundEstimator series of tools) instead")
   PseudoJet subtracted_jet(const PseudoJet & jet,
                            const Selector & selector) const;
 
   /// return the subtracted pt, given rho
+  FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::subtracted_pt(...) is deprecated since FastJet 3.0. Use the Subtractor tool (with the BackgroundEstimator series of tools) instead")
   double subtracted_pt(const PseudoJet & jet,
                        const double rho,
@@ -288,4 +299,5 @@
   /// apply jet-by-jet (see the BackgroundEstimator and Subtractor
   /// tools for more generic usages)
+  FASTJET_DEPRECATED_MSG("ClusterSequenceAreaBase::subtracted_pt(...) is deprecated since FastJet 3.0. Use the Subtractor tool (with the BackgroundEstimator series of tools) instead")
   double subtracted_pt(const PseudoJet & jet,
                        const Selector & selector,
@@ -296,4 +308,31 @@
   void _check_selector_good_for_median(const Selector &selector) const;
 
+  // the following set of protected methods are all deprecated. Their
+  // role is simply to hide the corresponding methods (without the
+  // first underscore) from the public interface so that they can be
+  // used internally until all the deprecated methods are removed.
+  // DO NOT USE ANY OF THESE METHODS: THEY ARE DEPRECATED AND WILL BE
+  // REMOVED.
+  virtual void _get_median_rho_and_sigma(const Selector & selector, 
+                                         bool use_area_4vector,
+                                         double & median, double & sigma,
+                                         double & mean_area) const;
+  virtual void _get_median_rho_and_sigma(const std::vector<PseudoJet> & all_jets,
+                                         const Selector & selector, 
+                                         bool use_area_4vector,
+                                         double & median, double & sigma,
+                                         double & mean_area,
+                                         bool all_are_inclusive = false) const;
+  virtual void _get_median_rho_and_sigma(const Selector & selector, 
+                                         bool use_area_4vector,
+                                         double & median, double & sigma) const {
+    double mean_area;
+    _get_median_rho_and_sigma(selector,  use_area_4vector,
+                              median,  sigma, mean_area);
+  }
+  virtual void _parabolic_pt_per_unit_area(double & a, double & b, 
+                                           const Selector & selector, 
+                                           double exclude_above=-1.0, 
+                                           bool use_area_4vector=false) const;
 
 private:
@@ -305,5 +344,24 @@
   /// check the jet algorithm is suitable (and if not issue a warning)
   void _check_jet_alg_good_for_median() const;
-  
+
+  // the following set of private methods are all deprecated. Their
+  // role is simply to hide the corresponding methods (without the
+  // first underscore) from the public interface so that they can be
+  // used internally until all the deprecated methods are removed.
+  // DO NOT USE ANY OF THESE METHODS: THEY ARE DEPRECATED AND WILL BE
+  // REMOVED.
+  double _median_pt_per_unit_area(const Selector & selector) const;
+  double _median_pt_per_unit_area_4vector(const Selector & selector) const;
+  double _median_pt_per_unit_something(const Selector & selector,
+                                       bool use_area_4vector) const;
+  std::vector<PseudoJet> _subtracted_jets(const double rho,
+                                          const double ptmin=0.0) const;
+  PseudoJet _subtracted_jet(const PseudoJet & jet,
+                            const double rho) const;
+  PseudoJet _subtracted_jet(const PseudoJet & jet,
+                            const Selector & selector) const;
+  double _subtracted_pt(const PseudoJet & jet,
+                        const double rho,
+                        bool use_area_4vector=false) const;
 };
 
Index: external/fastjet/ClusterSequencePassiveArea.hh
===================================================================
--- external/fastjet/ClusterSequencePassiveArea.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequencePassiveArea.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequencePassiveArea.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClusterSequencePassiveArea.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -65,5 +65,5 @@
   /// return an empty area that's appropriate to the passive area
   /// determination carried out
-  virtual double empty_area(const Selector & selector) const;
+  virtual double empty_area(const Selector & selector) const FASTJET_OVERRIDE;
 
 private:
Index: external/fastjet/ClusterSequenceStructure.hh
===================================================================
--- external/fastjet/ClusterSequenceStructure.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequenceStructure.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequenceStructure.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClusterSequenceStructure.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -76,5 +76,7 @@
 
   /// description
-  virtual std::string description() const{ return "PseudoJet with an associated ClusterSequence"; }
+  virtual std::string description() const FASTJET_OVERRIDE{
+    return "PseudoJet with an associated ClusterSequence";
+  }
 
   //-------------------------------------------------------------
@@ -85,21 +87,21 @@
   //-------------------------------------------------------------
   /// returns true if there is an associated ClusterSequence
-  virtual bool has_associated_cluster_sequence() const{ return true;}
+  virtual bool has_associated_cluster_sequence() const FASTJET_OVERRIDE{ return true;}
 
   /// get a (const) pointer to the parent ClusterSequence (NULL if
   /// inexistent)
-  virtual const ClusterSequence* associated_cluster_sequence() const;
+  virtual const ClusterSequence* associated_cluster_sequence() const FASTJET_OVERRIDE;
   
   /// returns true if there is a valid associated ClusterSequence
-  virtual bool has_valid_cluster_sequence() const;
+  virtual bool has_valid_cluster_sequence() const FASTJET_OVERRIDE;
 
   /// if the jet has a valid associated cluster sequence then return a
   /// pointer to it; otherwise throw an error
-  virtual const ClusterSequence * validated_cs() const;
+  virtual const ClusterSequence * validated_cs() const FASTJET_OVERRIDE;
 
 #ifndef __FJCORE__
   /// if the jet has valid area information then return a pointer to
   /// the associated ClusterSequenceAreaBase object; otherwise throw an error
-  virtual const ClusterSequenceAreaBase * validated_csab() const;
+  virtual const ClusterSequenceAreaBase * validated_csab() const FASTJET_OVERRIDE;
 #endif  // __FJCORE__
 
@@ -125,5 +127,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const;
+  virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const FASTJET_OVERRIDE;
 
   /// check if it has been recombined with another PseudoJet in which
@@ -133,5 +135,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const;
+  virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const FASTJET_OVERRIDE;
 
   /// check if it is the product of a recombination, in which case
@@ -141,5 +143,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual bool has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const;
+  virtual bool has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const FASTJET_OVERRIDE;
 
   /// check if the reference PseudoJet is contained in the second one
@@ -151,5 +153,5 @@
   /// false is returned if the 2 PseudoJet do not belong the same
   /// ClusterSequence
-  virtual bool object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const;
+  virtual bool object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const FASTJET_OVERRIDE;
 
   /// return true if the structure supports constituents. 
@@ -157,5 +159,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual bool has_constituents() const;
+  virtual bool has_constituents() const FASTJET_OVERRIDE;
 
   /// retrieve the constituents. 
@@ -163,5 +165,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const;
+  virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
 
@@ -170,5 +172,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual bool has_exclusive_subjets() const;
+  virtual bool has_exclusive_subjets() const FASTJET_OVERRIDE;
 
   /// return a vector of all subjets of the current jet (in the sense
@@ -183,5 +185,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference, const double & dcut) const;
+  virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference, const double & dcut) const FASTJET_OVERRIDE;
 
   /// return the size of exclusive_subjets(...); still n ln n with same
@@ -191,5 +193,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual int n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const;
+  virtual int n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const FASTJET_OVERRIDE;
 
   /// return the list of subjets obtained by unclustering the supplied
@@ -201,5 +203,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const;
+  virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const FASTJET_OVERRIDE;
 
   /// return the dij that was present in the merging nsub+1 -> nsub 
@@ -208,5 +210,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const;
+  virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const FASTJET_OVERRIDE;
 
   /// return the maximum dij that occurred in the whole event at the
@@ -216,5 +218,5 @@
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
-  virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const;
+  virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const FASTJET_OVERRIDE;
 
 
@@ -224,5 +226,5 @@
   /// by convention, a jet associated with a ClusterSequence will have
   /// its parents as pieces
-  virtual bool has_pieces(const PseudoJet &reference) const;
+  virtual bool has_pieces(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
   /// by convention, a jet associated with a ClusterSequence will have
@@ -235,5 +237,5 @@
   /// sequence. If the cluster sequence has gone out of scope, an
   /// error will be thrown
-  virtual std::vector<PseudoJet> pieces(const PseudoJet &reference) const;
+  virtual std::vector<PseudoJet> pieces(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
 
@@ -244,22 +246,22 @@
 
   /// check if it has a defined area
-  virtual bool has_area() const;
+  virtual bool has_area() const FASTJET_OVERRIDE;
 
   /// return the jet (scalar) area.
   /// throws an Error if there is no support for area in the parent CS
-  virtual double area(const PseudoJet &reference) const;
+  virtual double area(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
   /// return the error (uncertainty) associated with the determination
   /// of the area of this jet.
   /// throws an Error if there is no support for area in the parent CS
-  virtual double area_error(const PseudoJet &reference) const;
+  virtual double area_error(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
   /// return the jet 4-vector area.
   /// throws an Error if there is no support for area in the parent CS
-  virtual PseudoJet area_4vector(const PseudoJet &reference) const;
+  virtual PseudoJet area_4vector(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
   /// true if this jet is made exclusively of ghosts.
   /// throws an Error if there is no support for area in the parent CS
-  virtual bool is_pure_ghost(const PseudoJet &reference) const;
+  virtual bool is_pure_ghost(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
 #endif  // __FJCORE__
Index: external/fastjet/ClusterSequenceVoronoiArea.hh
===================================================================
--- external/fastjet/ClusterSequenceVoronoiArea.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequenceVoronoiArea.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequenceVoronoiArea.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClusterSequenceVoronoiArea.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -65,5 +65,5 @@
 
   /// return the area associated with the given jet
-  virtual inline double area(const PseudoJet & jet) const {
+  virtual inline double area(const PseudoJet & jet) const FASTJET_OVERRIDE {
     return _voronoi_area[jet.cluster_hist_index()];}
 
@@ -72,10 +72,10 @@
   /// made of sums of centres of all Voronoi cells in jet, each
   /// contributing with a normalisation equal to the area of the cell
-  virtual inline PseudoJet area_4vector(const PseudoJet & jet) const {
+  virtual inline PseudoJet area_4vector(const PseudoJet & jet) const FASTJET_OVERRIDE {
     return _voronoi_area_4vector[jet.cluster_hist_index()];}
 
   /// return the error of the area associated with the given jet
   /// (0 by definition for a voronoi area)
-  virtual inline double area_error(const PseudoJet & /*jet*/) const {
+  virtual inline double area_error(const PseudoJet & /*jet*/) const FASTJET_OVERRIDE {
     return 0.0;}
 
Index: external/fastjet/ClusterSequence_CP2DChan.cc
===================================================================
--- external/fastjet/ClusterSequence_CP2DChan.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequence_CP2DChan.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequence_CP2DChan.cc 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClusterSequence_CP2DChan.cc 4045 2016-03-03 10:01:55Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -46,6 +46,6 @@
   public:
     int orig, mirror;
-    MirrorInfo(int a, int b) : orig(a), mirror(b) {};
-    MirrorInfo() {};
+    MirrorInfo(int a, int b) : orig(a), mirror(b) {}
+    MirrorInfo() : orig(0), mirror(0) {} // set dummy values to keep static code checkers happy
   };
 
Index: external/fastjet/ClusterSequence_Delaunay.cc
===================================================================
--- external/fastjet/ClusterSequence_Delaunay.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/ClusterSequence_Delaunay.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClusterSequence_Delaunay.cc 3918 2015-07-03 14:19:13Z salam $
+// $Id: ClusterSequence_Delaunay.cc 4059 2016-03-03 20:49:48Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -69,6 +69,6 @@
 
   // initialise our DNN structure with the set of points
-  auto_ptr<DynamicNearestNeighbours> DNN;
-  bool verbose = false;
+  SharedPtr<DynamicNearestNeighbours> DNN;
+  const bool verbose = false;
 #ifndef DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
   bool ignore_nearest_is_mirror = (_Rparam < twopi);
Index: external/fastjet/CompositeJetStructure.hh
===================================================================
--- external/fastjet/CompositeJetStructure.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/CompositeJetStructure.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: CompositeJetStructure.hh 3652 2014-09-03 13:31:13Z salam $
+// $Id: CompositeJetStructure.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -64,5 +64,5 @@
 
   /// description
-  virtual std::string description() const;
+  virtual std::string description() const FASTJET_OVERRIDE;
 
   // things reimplemented from the base structure
@@ -70,5 +70,5 @@
   /// true unless the jet has no pieces (see also the description of
   /// constituents() below)
-  virtual bool has_constituents() const;
+  virtual bool has_constituents() const FASTJET_OVERRIDE;
 
   /// return the constituents (i.e. the union of the constituents of each piece)
@@ -78,5 +78,5 @@
   /// Note that as a consequence, a composite jet with no pieces will
   /// have an empty vector as constituents
-  virtual std::vector<PseudoJet> constituents(const PseudoJet &jet) const;
+  virtual std::vector<PseudoJet> constituents(const PseudoJet &jet) const FASTJET_OVERRIDE;
 
   //-------------------------------------------------------------------
@@ -84,8 +84,8 @@
   //-------------------------------------------------------------------
   /// true if it has pieces (always the case)
-  virtual bool has_pieces(const PseudoJet & /*jet*/) const {return true;}
+  virtual bool has_pieces(const PseudoJet & /*jet*/) const FASTJET_OVERRIDE {return true;}
 
   /// returns the pieces
-  virtual std::vector<PseudoJet> pieces(const PseudoJet &jet) const;
+  virtual std::vector<PseudoJet> pieces(const PseudoJet &jet) const FASTJET_OVERRIDE;
 
   // area-related material
@@ -93,8 +93,8 @@
 
   /// check if it has a well-defined area
-  virtual bool has_area() const;
+  virtual bool has_area() const FASTJET_OVERRIDE;
 
   /// return the jet (scalar) area.
-  virtual double area(const PseudoJet &reference) const;
+  virtual double area(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
   /// return the error (uncertainty) associated with the determination
@@ -102,13 +102,13 @@
   ///
   /// Be conservative: return the sum of the errors
-  virtual double area_error(const PseudoJet &reference) const;
+  virtual double area_error(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
   /// return the jet 4-vector area.
-  virtual PseudoJet area_4vector(const PseudoJet &reference) const;
+  virtual PseudoJet area_4vector(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
   /// true if this jet is made exclusively of ghosts.
   ///
   /// In this case, it will be true if all pieces are pure ghost
-  virtual bool is_pure_ghost(const PseudoJet &reference) const;
+  virtual bool is_pure_ghost(const PseudoJet &reference) const FASTJET_OVERRIDE;
 
   //unused: // allows one to modify the area information
Index: external/fastjet/DnnPlane.cc
===================================================================
--- external/fastjet/DnnPlane.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/DnnPlane.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: DnnPlane.cc 3918 2015-07-03 14:19:13Z salam $
+// $Id: DnnPlane.cc 3917 2015-07-03 14:07:50Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
Index: external/fastjet/Error.hh
===================================================================
--- external/fastjet/Error.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/Error.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: Error.hh 3809 2015-02-20 13:05:13Z soyez $
+// $Id: Error.hh 3807 2015-02-20 11:16:55Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -114,5 +114,4 @@
 };
 
-
 FASTJET_END_NAMESPACE
 
Index: external/fastjet/GhostedAreaSpec.hh
===================================================================
--- external/fastjet/GhostedAreaSpec.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/GhostedAreaSpec.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: GhostedAreaSpec.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: GhostedAreaSpec.hh 4074 2016-03-08 09:09:25Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -39,4 +39,5 @@
 #include "fastjet/Selector.hh"
 #include "fastjet/LimitedWarning.hh"
+#include "fastjet/internal/deprecated.hh"
 
 // 
@@ -173,4 +174,5 @@
   ///
   /// FJ2 placement is now deprecated.
+  FASTJET_DEPRECATED_MSG("This is deprecated since we strongly recomment to use the new ghost placement instead")
   void set_fj2_placement(bool  val);
 
Index: external/fastjet/JetDefinition.cc
===================================================================
--- external/fastjet/JetDefinition.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/JetDefinition.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: JetDefinition.cc 3677 2014-09-09 22:45:25Z soyez $
+// $Id: JetDefinition.cc 4074 2016-03-08 09:09:25Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -45,6 +45,6 @@
 JetDefinition::JetDefinition(JetAlgorithm jet_algorithm_in, 
 			     double R_in, 
+			     RecombinationScheme recomb_scheme_in,
 			     Strategy strategy_in,
-			     RecombinationScheme recomb_scheme_in,
                              int nparameters) :
   _jet_algorithm(jet_algorithm_in), _Rparam(R_in), _strategy(strategy_in) {
@@ -185,5 +185,5 @@
 
   // do not forget to delete the existing recombiner if needed
-  if (_shared_recombiner()) _shared_recombiner.reset();
+  if (_shared_recombiner) _shared_recombiner.reset();
 
   _recombiner = 0;
Index: external/fastjet/JetDefinition.hh
===================================================================
--- external/fastjet/JetDefinition.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/JetDefinition.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: JetDefinition.hh 3677 2014-09-09 22:45:25Z soyez $
+// $Id: JetDefinition.hh 4074 2016-03-08 09:09:25Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -35,4 +35,5 @@
 #include "fastjet/internal/numconsts.hh"
 #include "fastjet/PseudoJet.hh"
+#include "fastjet/internal/deprecated.hh"
 #include<string>
 #include<memory>
@@ -268,5 +269,5 @@
                 RecombinationScheme recomb_scheme_in = E_scheme,
                 Strategy strategy_in = Best) {
-    *this = JetDefinition(jet_algorithm_in, R_in, strategy_in, recomb_scheme_in, 1);
+    *this = JetDefinition(jet_algorithm_in, R_in, recomb_scheme_in, strategy_in, 1);
   }
 
@@ -277,5 +278,5 @@
                 Strategy strategy_in = Best) {
     double dummyR = 0.0;
-    *this = JetDefinition(jet_algorithm_in, dummyR, strategy_in, recomb_scheme_in, 0);
+    *this = JetDefinition(jet_algorithm_in, dummyR, recomb_scheme_in, strategy_in, 0);
   }
 
@@ -287,5 +288,5 @@
                 RecombinationScheme recomb_scheme_in = E_scheme,
                 Strategy strategy_in = Best) {
-    *this = JetDefinition(jet_algorithm_in, R_in, strategy_in, recomb_scheme_in, 2);
+    *this = JetDefinition(jet_algorithm_in, R_in, recomb_scheme_in, strategy_in, 2);
     set_extra_param(xtra_param_in);
   }
@@ -344,19 +345,31 @@
     _strategy = plugin_strategy;
     _Rparam = _plugin->R();
+    _extra_param = 0.0; // a dummy value to keep static code checkers happy
     _jet_algorithm = plugin_algorithm;
     set_recombination_scheme(E_scheme);
   }
 
-
   /// constructor to fully specify a jet-definition (together with
   /// information about how algorithically to run it).
+  JetDefinition(JetAlgorithm jet_algorithm_in, 
+                double R_in, 
+                RecombinationScheme recomb_scheme_in,
+                Strategy strategy_in,
+                int nparameters_in);
+
+  /// constructor to fully specify a jet-definition (together with
+  /// information about how algorithically to run it).
   ///
   /// the ordering of arguments here is old and deprecated (except
   /// as the common constructor for internal use)
+  FASTJET_DEPRECATED_MSG("This argument ordering is deprecated. Use JetDefinition(alg, R, strategy, scheme[, n_parameters]) instead")
   JetDefinition(JetAlgorithm jet_algorithm_in, 
                 double R_in, 
                 Strategy strategy_in,
                 RecombinationScheme recomb_scheme_in = E_scheme,
-                int nparameters_in = 1);
+                int nparameters_in = 1){
+    (*this) = JetDefinition(jet_algorithm_in,R_in,recomb_scheme_in,strategy_in,nparameters_in);
+  }
+
 
   /// cluster the supplied particles and returns a vector of resulting
@@ -388,5 +401,5 @@
   /// Recombiner *) may lead to memory corruption.
   void set_recombiner(const Recombiner * recomb) {
-    if (_shared_recombiner()) _shared_recombiner.reset(recomb);
+    if (_shared_recombiner) _shared_recombiner.reset(recomb);
     _recombiner = recomb;
     _default_recombiner = DefaultRecombiner(external_scheme);
@@ -521,11 +534,11 @@
       _recomb_scheme(recomb_scheme) {}
     
-    virtual std::string description() const;
+    virtual std::string description() const FASTJET_OVERRIDE;
     
     /// recombine pa and pb and put result into pab
     virtual void recombine(const PseudoJet & pa, const PseudoJet & pb, 
-                           PseudoJet & pab) const;
-
-    virtual void preprocess(PseudoJet & p) const;
+                           PseudoJet & pab) const FASTJET_OVERRIDE;
+
+    virtual void preprocess(PseudoJet & p) const FASTJET_OVERRIDE;
 
     /// return the index of the recombination scheme
Index: external/fastjet/LazyTiling25.cc
===================================================================
--- external/fastjet/LazyTiling25.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/LazyTiling25.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: LazyTiling25.cc 3808 2015-02-20 11:24:53Z soyez $
+// $Id: LazyTiling25.cc 3807 2015-02-20 11:16:55Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -369,6 +369,6 @@
     double dist = _distance_to_tile(jet, *near_tile) - tile_edge_security_margin;
     // cout << "      max info looked at tile " << *near_tile - &_tiles[0] 
-    // 	 << ", dist = " << dist << " " << (*near_tile)->max_NN_dist
-    // 	 << endl;
+    //      << ", dist = " << dist << " " << (*near_tile)->max_NN_dist
+    //      << " -> diff = " << dist-(*near_tile)->max_NN_dist << endl;
     if (dist > (*near_tile)->max_NN_dist) continue;
 
@@ -665,4 +665,31 @@
     int n_near_tiles = 0;
 
+    // GS comment:
+    //
+    // At this stage, we have perforned the clustering in
+    // ClusterSequence and we need to update the NNs. The objects we
+    // deal with are jetA and oldB (the once that have been clustered)
+    // as well as jetB (the result of the clustering)
+    //
+    // There are two types of objects we need to update:
+    //  - find jetB NN
+    //  - update the NN of points which had jetA or jetB as their NN
+    //
+    // Wile we find jetB's NN, browsing relevant tiles near jetB, we
+    // also search for points which had jetA or jetB as their
+    // NN. These are tagged. Then we list the relevant tiles where we
+    // can potentially have points to update (i.e. points which had
+    // jetA and oldB as NN) in the yet untagged relevant tiles near
+    // jetA and oldB.
+    //
+    // DEBUG:
+    // if (jetB != NULL) {
+    //   cout << "jetA = " << jetA->_jets_index << " (tile " << jetA->tile_index << "), "
+    //        << "oldB = " << oldB._jets_index  << " (tile " << oldB. tile_index << "), "
+    //        << "jetB = " << jetB->_jets_index << " (tile " << jetB->tile_index << ")" << endl;
+    // } else {
+    //   cout << "jetA = " << jetA->_jets_index << " (tile " << jetA->tile_index << ")" << endl;
+    // }
+    
     // Initialise jetB's NN distance as well as updating it for other
     // particles. While doing so, examine whether jetA or old jetB was
@@ -679,4 +706,8 @@
     	bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
         bool relevant = relevant_for_jetB || relevant_for_near_tile;
+
+        // cout << "  Relevance of tile " << *near_tile - & _tiles[0]
+        //      << " wrt jetB is " << relevant << endl;
+
         if (! relevant) continue;
         // now label this tile as having been considered (so that we 
@@ -702,7 +733,9 @@
     // and one new jet.
     int n_done_tiles = n_near_tiles;
+    //cout << "Looking at relevant tiles to update for jetA" << endl;
     _add_untagged_neighbours_to_tile_union_using_max_info(jetA, 
        					   tile_union, n_near_tiles);
     if (jetB != NULL) {
+      // cout << "Looking at relevant tiles to update for oldB" << endl;
 	_add_untagged_neighbours_to_tile_union_using_max_info(&oldB,
 							      tile_union,n_near_tiles);
Index: external/fastjet/LazyTiling9.cc
===================================================================
--- external/fastjet/LazyTiling9.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/LazyTiling9.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: LazyTiling9.cc 3808 2015-02-20 11:24:53Z soyez $
+// $Id: LazyTiling9.cc 3807 2015-02-20 11:16:55Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
Index: external/fastjet/LazyTiling9Alt.cc
===================================================================
--- external/fastjet/LazyTiling9Alt.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/LazyTiling9Alt.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: LazyTiling9Alt.cc 3808 2015-02-20 11:24:53Z soyez $
+// $Id: LazyTiling9Alt.cc 3807 2015-02-20 11:16:55Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
Index: external/fastjet/LazyTiling9SeparateGhosts.cc
===================================================================
--- external/fastjet/LazyTiling9SeparateGhosts.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/LazyTiling9SeparateGhosts.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: LazyTiling9SeparateGhosts.cc 3808 2015-02-20 11:24:53Z soyez $
+// $Id: LazyTiling9SeparateGhosts.cc 3807 2015-02-20 11:16:55Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
Index: external/fastjet/NNBase.hh
===================================================================
--- external/fastjet/NNBase.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
+++ external/fastjet/NNBase.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -0,0 +1,192 @@
+#ifndef __FASTJET_NNBASE_HH__
+#define __FASTJET_NNBASE_HH__
+
+//FJSTARTHEADER
+// $Id: NNBase.hh 4060 2016-03-03 21:46:49Z salam $
+//
+// Copyright (c) 2016, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet.
+//
+//  FastJet 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.
+//
+//  The algorithms that underlie FastJet have required considerable
+//  development. They are described in the original FastJet paper,
+//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
+//  FastJet as part of work towards a scientific publication, please
+//  quote the version you use and include a citation to the manual and
+//  optionally also to hep-ph/0512210.
+//
+//  FastJet 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 FastJet. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+//FJENDHEADER
+
+#include<fastjet/ClusterSequence.hh>
+
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+/// @ingroup advanced_usage
+/// \class _NoInfo
+/// internal dummy class, used as a default template argument
+class _NoInfo {};
+
+/// @ingroup advanced_usage
+/// \class NNInfo
+///
+/// internal helper template class to facilitate initialisation of a
+/// BJ with a PseudoJet and extra information. Implementations of
+/// NN-based clustering do not need to explicitly use or refer to
+/// this class!
+template<class I> class NNInfo {
+public:
+  NNInfo()         : _info(NULL) {}
+  NNInfo(I * info) : _info(info) {}
+  template<class BJ> void init_jet(BJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index, _info);}
+private:
+  I * _info;
+};
+
+/// @ingroup advanced_usage Internal helper specialisation of NNInfo
+/// for cases where there is no extra info
+template<> class NNInfo<_NoInfo>  {
+public:
+  NNInfo()           {}
+  NNInfo(_NoInfo * ) {}
+  template<class BJ> void init_jet(BJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index);}
+};
+
+
+//----------------------------------------------------------------------
+/// @ingroup advanced_usage
+/// \class NNBase
+/// Helps solve closest pair problems with generic interparticle and
+/// particle-beam distances.
+///
+/// \section Description Description and derived classes:
+///
+///   This is an abstract base class which defines the interface for
+///   several classes that help carry out nearest-neighbour
+///   clustering:
+///  
+///    - NNH          provides an implementation for generic measures,
+///  
+///    - NNFJN2Plain  provides an implementation for distances
+///                   satisfying the FastJet lemma i.e. distances for
+///                   which the minimum dij has the property that i is
+///                   the geometrical nearest neighbour of j, or vice
+///                   versa. I.e. the distance can be factorised in a
+///                   momentum factor and a geometric piece. This is
+///                   based on the fastjet N2Plain clustering strategy
+///  
+///    - NNFJN2Tiled  is a tiled version of NNFJN2Plain (based on the
+///                   N2Tiled FastJet clustering strategy). Like
+///                   NNPlain2 it applies to distance measures that
+///                   satisfy the FastJet lemma, with the additional
+///                   restriction that: (a) the underlying geometry
+///                   should be cylindrical (e.g. rapidity--azimuth)
+///                   and (b) the search for the geometric nearest
+///                   neighbour of each particle can be limited to
+///                   that particle's tile and its neighbouring tiles.
+///
+/// If you can use NNFJN2Plain it will usually be faster than
+/// NNH. NNFJN2Tiled, where it can be used, will be faster for
+/// multiplicities above a few tens of particles.
+///
+///   NOTE: IN ALL CASES, THE DISTANCE MUST BE SYMMETRIC (dij=dji)!!!
+///
+/// \section BJ Underlying BriefJet (BJ) class:
+/// 
+///   All derived classes must be templated with a BriefJet (BJ)
+///   class --- BJ should basically cache the minimal amount of
+///   information that is needed to efficiently calculate
+///   interparticle distances and particle-beam distances.
+///   
+///   This class can be used with or without an extra "Information"
+///   template, i.e. `NN*<BJ>` or `NN*<BJ,I>`. Accordingly BJ must provide
+///   one of the two following init functions:
+///
+///   \code
+///     void  BJ::init(const PseudoJet & jet);            // initialise with a PseudoJet
+///     void  BJ::init(const PseudoJet & jet, I * info);  // initialise with a PseudoJet + info
+///   \endcode
+///
+///   where info might be a pointer to a class that contains, e.g.,
+///   information about R, or other parameters of the jet algorithm
+///   
+///   The BJ then provides information about interparticle and
+///   particle-beam distances. The exact requirements depend on
+///   whether you use NNH, NNFJN2Plain or NNFJN2Tiled. (See the
+///   corresponding classes for details).
+///
+///
+/// \section Workflow Workflow:
+///
+///   In all cases, the usage of NNBase classes works as follows:
+///
+///   First, from the list of particles, create an `NN*<BJ>`
+///   object of the appropriate type with the appropriate BJ class
+///   (and optional extra info).
+///
+///   Then, cluster using a loop like this (assuming a FastJet plugin)
+///
+///   \code
+///     while (njets > 0) {
+///       int i, j, k;
+///       // get the i and j that minimize the distance
+///       double dij = nn.dij_min(i, j);  
+///     
+///       // do the appropriate recombination and update the nn
+///       if (j >= 0) {    // interparticle recombination
+///         cs.plugin_record_ij_recombination(i, j, dij, k);
+///         nn.merge_jets(i, j, cs.jets()[k], k); 
+///       } else {         // bbeam recombination
+///         double diB = cs.jets()[i].E()*cs.jets()[i].E(); // get new diB
+///         cs.plugin_record_iB_recombination(i, diB);
+///         nn.remove_jet(i);
+///       }
+///       njets--;
+///     }
+///   \endcode
+///
+/// For an example of how the NNH<BJ> class is used, see the
+/// JadePlugin or EECambridgePlugin.
+template<class I = _NoInfo> class NNBase : public NNInfo<I> {
+public:
+  /// Default constructor
+  NNBase() {}
+  /// Constuctor with additional Info 
+  NNBase(I * info) : NNInfo<I>(info) {}
+
+  /// initialisation from a given list of particles
+  virtual void start(const std::vector<PseudoJet> & jets) = 0;
+  
+  /// returns the dij_min and indices iA, iB, for the corresponding jets.
+  /// If iB < 0 then iA recombines with the beam
+  virtual double dij_min(int & iA, int & iB) = 0;
+
+  /// removes the jet pointed to by index iA
+  virtual void remove_jet(int iA) = 0;
+
+  /// merges the jets pointed to by indices A and B and replaces them with
+  /// jet, assigning it an index jet_index.
+  virtual void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index) =  0;
+
+  virtual ~NNBase() {};
+};
+
+
+FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
+
+
+#endif // __FASTJET_NNBASE_HH__
Index: external/fastjet/NNFJN2Plain.hh
===================================================================
--- external/fastjet/NNFJN2Plain.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
+++ external/fastjet/NNFJN2Plain.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -0,0 +1,401 @@
+#ifndef __FASTJET_NNFJN2PLAIN_HH__
+#define __FASTJET_NNFJN2PLAIN_HH__
+
+//FJSTARTHEADER
+// $Id: NNFJN2Plain.hh 4058 2016-03-03 15:39:40Z salam $
+//
+// Copyright (c) 2016, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet.
+//
+//  FastJet 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.
+//
+//  The algorithms that underlie FastJet have required considerable
+//  development. They are described in the original FastJet paper,
+//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
+//  FastJet as part of work towards a scientific publication, please
+//  quote the version you use and include a citation to the manual and
+//  optionally also to hep-ph/0512210.
+//
+//  FastJet 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 FastJet. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+//FJENDHEADER
+
+#include <fastjet/NNBase.hh>
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+//----------------------------------------------------------------------
+/// @ingroup advanced_usage
+/// \class NNFJN2Plain
+///
+/// Helps solve closest pair problems with factorised interparticle
+/// and beam distances (ie satisfying the FastJet lemma)
+///
+/// (see NNBase.hh for an introductory description)
+///
+/// This variant provides an implementation based on the N2Plain
+/// clustering strategy in FastJet. The interparticle and beam
+/// distances should be of the form
+///
+/// \code
+///   dij = min(mom_factor(i), mom_factor(j)) * geometrical_distance(i,j)
+///   diB = mom_factor(i) * geometrical_beam_distance(i)
+/// \endcode
+///
+/// The class is templated with a BJ (brief jet) class and can be used
+/// with or without an extra "Information" template,
+/// i.e. NNFJN2Plain<BJ> or NNFJN2Plain<BJ,I>
+///
+/// For the NNH_N2Plain<BJ> version of the class to function, BJ must
+/// provide four member functions
+///  
+/// \code
+///   void   BJ::init(const PseudoJet & jet);                   // initialise with a PseudoJet
+///   double BJ::geometrical_distance(const BJ * other_bj_jet); // distance between this and other_bj_jet (geometrical part)
+///   double BJ::geometrical_beam_distance();                   // distance to the beam (geometrical part)
+///   double BJ::momentum_factor();                             // extra momentum factor
+/// \endcode
+///
+/// For the NNH_N2Plain<BJ,I> version to function, the BJ::init(...)
+/// member must accept an extra argument
+///
+/// \code
+///   void  BJ::init(const PseudoJet & jet, I * info);   // initialise with a PseudoJet + info
+/// \endcode
+///
+/// NOTE: THE DISTANCE MUST BE SYMMETRIC I.E. SATISFY
+/// \code
+///     a.geometrical_distance(b) == b.geometrical_distance(a)
+/// \endcode
+///
+/// Note that you are strongly advised to add the following lines to
+/// your BJ class to allow it to be used also with NNH:
+///
+/// \code
+///   /// make this BJ class compatible with the use of NNH
+///   double BJ::distance(const BJ * other_bj_jet){
+///     double mom1 = momentum_factor();
+///     double mom2 = other_bj_jet->momentum_factor();
+///     return (mom1<mom2 ? mom1 : mom2) * geometrical_distance(other_bj_jet);
+///   }
+///   double BJ::beam_distance(){
+///     return momentum_factor() * geometrical_beam_distance();
+///   }
+/// \endcode
+///
+template<class BJ, class I = _NoInfo> class NNFJN2Plain : public NNBase<I> {
+public:
+
+  /// constructor with an initial set of jets (which will be assigned indices
+  /// `0...jets.size()-1`)
+  NNFJN2Plain(const std::vector<PseudoJet> & jets)           : NNBase<I>()     {start(jets);}
+  NNFJN2Plain(const std::vector<PseudoJet> & jets, I * info) : NNBase<I>(info) {start(jets);}
+  
+  /// initialisation from a given list of particles
+  virtual void start(const std::vector<PseudoJet> & jets);
+
+  /// returns the dij_min and indices iA, iB, for the corresponding jets.
+  /// If iB < 0 then iA recombines with the beam
+  double dij_min(int & iA, int & iB);
+
+  /// removes the jet pointed to by index iA
+  void remove_jet(int iA);
+
+  /// merges the jets pointed to by indices A and B and replace them with
+  /// jet, assigning it an index jet_index.
+  void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
+
+  /// a destructor
+  ~NNFJN2Plain() {
+    delete[] briefjets;
+    delete[] diJ;
+  }
+    
+private:
+  class NNBJ; // forward declaration
+
+  // return the full distance of a particle to its NN
+  inline double compute_diJ(const NNBJ * const jet) const {
+    double mom_fact = jet->momentum_factor();
+    if (jet->NN != NULL) {
+      double other_mom_fact = jet->NN->momentum_factor();
+      if (other_mom_fact < mom_fact) {mom_fact = other_mom_fact;}
+    }
+    return jet->NN_dist * mom_fact;
+  }
+
+  /// establish the nearest neighbour for jet, and cross check consistency
+  /// of distances for the other jets that are encountered. Assumes
+  /// jet not contained within begin...end
+  void set_NN_crosscheck(NNBJ * jet, NNBJ * begin, NNBJ * end);
+
+  /// establish the nearest neighbour for jet; don't cross check other jets'
+  /// distances and allow jet to be contained within begin...end
+  void set_NN_nocross   (NNBJ * jet, NNBJ * begin, NNBJ * end);
+
+  /// contains the briefjets
+  NNBJ * briefjets;
+
+  /// semaphores for the current extent of our structure
+  NNBJ * head, * tail;
+
+  /// currently active number of jets
+  int n;
+
+  /// where_is[i] contains a pointer to the jet with index i
+  std::vector<NNBJ *> where_is;
+
+  /// a table containing all the (full) distances to each object's NN
+  double * diJ;
+
+  /// a class that wraps around the BJ, supplementing it with extra information
+  /// such as pointers to neighbours, etc.
+  class NNBJ : public BJ {
+  public:
+    void init(const PseudoJet & jet, int index_in) {
+      BJ::init(jet);
+      other_init(index_in);
+    }
+    void init(const PseudoJet & jet, int index_in, I * info) {
+      BJ::init(jet, info);
+      other_init(index_in);
+    }
+    void other_init(int index_in) {
+      _index = index_in;
+      NN_dist = BJ::geometrical_beam_distance();
+      NN = NULL;
+    }
+    int index() const {return _index;}
+    
+    double NN_dist;
+    NNBJ * NN;
+    
+  private:
+    int _index;
+  };
+
+};
+
+
+
+//----------------------------------------------------------------------
+template<class BJ, class I> void NNFJN2Plain<BJ,I>::start(const std::vector<PseudoJet> & jets) {
+  n = jets.size();
+  briefjets = new NNBJ[n];
+  where_is.resize(2*n);
+
+  NNBJ * jetA = briefjets;
+  
+  // initialise the basic jet info 
+  for (int i = 0; i< n; i++) {
+    // the this-> in the next line is required by standard compiler
+    // see e.g. http://stackoverflow.com/questions/10639053/name-lookups-in-c-templates
+    this->init_jet(jetA, jets[i], i);
+    where_is[i] = jetA;
+    jetA++; // move on to next entry of briefjets
+  }
+  tail = jetA; // a semaphore for the end of briefjets
+  head = briefjets; // a nicer way of naming start
+
+  // now initialise the NN distances: jetA will run from 1..n-1; and
+  // jetB from 0..jetA-1
+  for (jetA = head + 1; jetA != tail; jetA++) {
+    // set NN info for jetA based on jets running from head..jetA-1,
+    // checking in the process whether jetA itself is an undiscovered
+    // NN of one of those jets.
+    set_NN_crosscheck(jetA, head, jetA);
+  }
+
+  // now create the diJ (where J is i's NN) table -- remember that 
+  // we differ from standard normalisation here by a factor of R2
+  diJ = new double[n];
+  jetA = head;
+  for (int i = 0; i < n; i++) {
+    diJ[i] = compute_diJ(jetA);
+    jetA++; // have jetA follow i
+  }
+}
+
+
+//----------------------------------------------------------------------
+template<class BJ, class I> double NNFJN2Plain<BJ,I>::dij_min(int & iA, int & iB) {
+  // find the minimum of the diJ on this round
+  double diJ_min = diJ[0];
+  int diJ_min_jet = 0;
+  for (int i = 1; i < n; i++) {
+    if (diJ[i] < diJ_min) {
+      diJ_min_jet = i;
+      diJ_min  = diJ[i];
+    }
+  }
+
+  // return information to user about recombination
+  NNBJ * jetA = & briefjets[diJ_min_jet];
+  //std::cout << jetA - briefjets << std::endl; 
+  iA = jetA->index();
+  iB = jetA->NN ? jetA->NN->index() : -1;
+  return diJ_min;
+}
+
+
+//----------------------------------------------------------------------
+// remove jetA from the list
+template<class BJ, class I> void NNFJN2Plain<BJ,I>::remove_jet(int iA) {
+  NNBJ * jetA = where_is[iA];
+  // now update our nearest neighbour info and diJ table
+  // first reduce size of table
+  tail--; n--;
+  // Copy last jet contents and diJ info into position of jetA
+  *jetA = *tail;
+  // update the info on where the given index is stored
+  where_is[jetA->index()] = jetA;
+  diJ[jetA - head] = diJ[tail-head];
+
+  // updating NN infos. 2 cases to watch for (see below)
+  for (NNBJ * jetI = head; jetI != tail; jetI++) {
+    // see if jetI had jetA as a NN -- if so recalculate the NN
+    if (jetI->NN == jetA) {
+      set_NN_nocross(jetI, head, tail);
+      diJ[jetI-head] = compute_diJ(jetI); // update diJ 
+    } 
+    // if jetI's NN is the new tail then relabel it so that it becomes jetA
+    if (jetI->NN == tail) {jetI->NN = jetA;}
+  }
+}
+
+
+//----------------------------------------------------------------------
+template<class BJ, class I> void NNFJN2Plain<BJ,I>::merge_jets(int iA, int iB, 
+					const PseudoJet & jet, int index) {
+
+  NNBJ * jetA = where_is[iA];
+  NNBJ * jetB = where_is[iB];
+
+  // If necessary relabel A & B to ensure jetB < jetA, that way if
+  // the larger of them == newtail then that ends up being jetA and 
+  // the new jet that is added as jetB is inserted in a position that
+  // has a future!
+  if (jetA < jetB) std::swap(jetA,jetB);
+
+  // initialise jetB based on the new jet
+  this->init_jet(jetB, jet, index);
+  // and record its position (making sure we have the space)
+  if (index >= int(where_is.size())) where_is.resize(2*index);
+  where_is[jetB->index()] = jetB;
+
+  // now update our nearest neighbour info
+  // first reduce size of table
+  tail--; n--;
+  // Copy last jet contents into position of jetA
+  *jetA = *tail;
+  // update the info on where the tail's index is stored
+  where_is[jetA->index()] = jetA;
+  diJ[jetA - head] = diJ[tail-head];
+
+  // initialise jetB NN's distance and update NN infos
+  for (NNBJ * jetI = head; jetI != tail; jetI++) {
+    // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
+    if (jetI->NN == jetA || jetI->NN == jetB) {
+      set_NN_nocross(jetI, head, tail);
+      diJ[jetI-head] = compute_diJ(jetI); // update diJ 
+    } 
+
+    // check whether new jetB is closer than jetI's current NN and
+    // if need be update things
+    double dist = jetI->geometrical_distance(jetB);
+    if (dist < jetI->NN_dist) {  // we need to update I
+      if (jetI != jetB) {
+        jetI->NN_dist = dist;
+        jetI->NN = jetB;
+        diJ[jetI-head] = compute_diJ(jetI); // update diJ...
+      }
+    }
+    if (dist < jetB->NN_dist) { // we need to update B
+      if (jetI != jetB) {
+        jetB->NN_dist = dist;
+        jetB->NN      = jetI;
+      }
+    }
+
+    // if jetI's NN is the new tail then relabel it so that it becomes jetA
+    if (jetI->NN == tail) {jetI->NN = jetA;}
+  }
+
+  // update info for jetB
+  diJ[jetB-head] = compute_diJ(jetB);
+}
+
+
+//----------------------------------------------------------------------
+// this function assumes that jet is not contained within begin...end
+template <class BJ, class I> void NNFJN2Plain<BJ,I>::set_NN_crosscheck(NNBJ * jet, 
+		    NNBJ * begin, NNBJ * end) {
+  double NN_dist = jet->geometrical_beam_distance();
+  NNBJ * NN      = NULL;
+  for (NNBJ * jetB = begin; jetB != end; jetB++) {
+    double dist = jet->geometrical_distance(jetB);
+    if (dist < NN_dist) {
+      NN_dist = dist;
+      NN = jetB;
+    }
+    if (dist < jetB->NN_dist) {
+      jetB->NN_dist = dist;
+      jetB->NN = jet;
+    }
+  }
+  jet->NN = NN;
+  jet->NN_dist = NN_dist;
+}
+
+
+//----------------------------------------------------------------------
+// set the NN for jet without checking whether in the process you might
+// have discovered a new nearest neighbour for another jet
+template <class BJ, class I>  void NNFJN2Plain<BJ,I>::set_NN_nocross(
+                 NNBJ * jet, NNBJ * begin, NNBJ * end) {
+  double NN_dist = jet->geometrical_beam_distance();
+  NNBJ * NN      = NULL;
+  // if (head < jet) {
+  //   for (NNBJ * jetB = head; jetB != jet; jetB++) {
+  if (begin < jet) {
+    for (NNBJ * jetB = begin; jetB != jet; jetB++) {
+      double dist = jet->geometrical_distance(jetB);
+      if (dist < NN_dist) {
+	NN_dist = dist;
+	NN = jetB;
+      }
+    }
+  }
+  // if (tail > jet) {
+  //   for (NNBJ * jetB = jet+1; jetB != tail; jetB++) {
+  if (end > jet) {
+    for (NNBJ * jetB = jet+1; jetB != end; jetB++) {
+      double dist = jet->geometrical_distance(jetB);
+      if (dist < NN_dist) {
+	NN_dist = dist;
+	NN = jetB;
+      }
+    }
+  }
+  jet->NN = NN;
+  jet->NN_dist = NN_dist;
+}
+
+
+
+
+FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
+
+
+#endif // __FASTJET_NNFJN2PLAIN_HH__
Index: external/fastjet/NNFJN2Tiled.hh
===================================================================
--- external/fastjet/NNFJN2Tiled.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
+++ external/fastjet/NNFJN2Tiled.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -0,0 +1,715 @@
+#ifndef __FASTJET_NNFJN2TILED_HH__
+#define __FASTJET_NNFJN2TILED_HH__
+
+//FJSTARTHEADER
+// $Id: NNFJN2Tiled.hh 4056 2016-03-03 15:27:35Z soyez $
+//
+// Copyright (c) 2016, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet.
+//
+//  FastJet 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.
+//
+//  The algorithms that underlie FastJet have required considerable
+//  development. They are described in the original FastJet paper,
+//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
+//  FastJet as part of work towards a scientific publication, please
+//  quote the version you use and include a citation to the manual and
+//  optionally also to hep-ph/0512210.
+//
+//  FastJet 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 FastJet. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+//FJENDHEADER
+
+#include <fastjet/NNBase.hh>
+#include <fastjet/internal/TilingExtent.hh>
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+//----------------------------------------------------------------------
+/// @ingroup advanced_usage
+/// \class NNFJN2Tiled
+///
+/// Helps solve closest pair problems with factorised interparticle
+/// and beam distances (ie satisfying the FastJet lemma) that are on
+/// a cylindrical geometry and allow tiling.
+///
+/// (see NNBase.hh for an introductory description)
+///
+/// This variant provides an implementation based on the N2Tiled
+/// clustering strategy in FastJet. As for the NNFJN2Plain case, the
+/// interparticle and beam distances should be of the form
+///
+/// \code
+///   dij = min(mom_factor(i), mom_factor(j)) * geometrical_distance(i,j)
+///   diB = mom_factor(i) * geometrical_beam_distance(i)
+/// \endcode
+///
+/// Additionally, the NNFJN2Tiled class takes a tile_size parameter
+/// that controls the size of the tiles. It must be such that, for any
+/// two points in non-neighbouring (and non-identical) tiles, the
+/// geometrical distance between the 2 points is larger than the
+/// geometrical beam distance of each of the 2 points.
+///
+/// It is templated with a BJ (brief jet) class and can be used with or
+/// without an extra "Information" template, i.e. NNFJN2Tiled<BJ> or
+/// NNFJN2Tiled<BJ,I>
+///
+/// For the NNFJN2Tiled<BJ> version of the class to function, BJ must provide 
+/// three member functions
+///  
+/// \code
+///   void   BJ::init(const PseudoJet & jet);                   // initialise with a PseudoJet
+///   double BJ::geometrical_distance(const BJ * other_bj_jet); // distance between this and other_bj_jet (geometrical part)
+///   double BJ::geometrical_beam_distance();                   // distance to the beam (geometrical part)
+///   double BJ::momentum_factor();                             // extra momentum factor
+/// \endcode
+///
+/// For the NNFJN2Tiled<BJ,I> version to function, the BJ::init(...) member
+/// must accept an extra argument
+///
+/// \code
+///   void BJ::init(const PseudoJet & jet, I * info);   // initialise with a PseudoJet + info
+/// \endcode
+///
+/// NOTE: THE DISTANCE MUST BE SYMMETRIC I.E. SATISFY
+/// \code
+///     a.geometrical_distance(b) == b.geometrical_distance(a)
+/// \endcode
+///
+/// Finally, the BJ class needs to provide access to the variables used
+/// for the rectangular tiling:
+///
+/// \code
+///   double BJ::rap(); // rapidity-like variable
+///   double BJ::phi(); // azimutal-angle-like variable (should be > -2pi)
+/// \endcode
+///
+/// Note that you are strongly advised to add the following lines
+/// to your BJ class to allow it to be used also with NNH:
+///
+/// \code
+///   /// make this BJ class compatible with the use of NNH
+///   double BJ::distance(const BJ * other_bj_jet){
+///     double mom1 = momentum_factor();
+///     double mom2 = other_bj_jet->momentum_factor();
+///     return (mom1<mom2 ? mom1 : mom2) * geometrical_distance(other_bj_jet);
+///   }
+///   double BJ::beam_distance(){
+///     return momentum_factor() * geometrical_beam_distance();
+///   }
+/// \endcode
+///
+template<class BJ, class I = _NoInfo> class NNFJN2Tiled : public NNBase<I> {
+public:
+
+  /// constructor with an initial set of jets (which will be assigned indices
+  /// `0...jets.size()-1`)
+  NNFJN2Tiled(const std::vector<PseudoJet> & jets, double requested_tile_size)
+    : NNBase<I>(),     _requested_tile_size(requested_tile_size) {start(jets);}
+  NNFJN2Tiled(const std::vector<PseudoJet> & jets, double requested_tile_size, I * info)
+    : NNBase<I>(info), _requested_tile_size(requested_tile_size) {start(jets);}
+  
+  void start(const std::vector<PseudoJet> & jets);
+
+  /// return the dij_min and indices iA, iB, for the corresponding jets.
+  /// If iB < 0 then iA recombines with the beam
+  double dij_min(int & iA, int & iB);
+
+  /// remove the jet pointed to by index iA
+  void remove_jet(int iA);
+
+  /// merge the jets pointed to by indices A and B and replace them with
+  /// jet, assigning it an index jet_index.
+  void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
+
+  /// a destructor
+  ~NNFJN2Tiled() {
+    delete[] briefjets;
+    delete[] diJ;
+  }
+    
+private:
+  class TiledJet; // forward declaration
+  class Tile;
+  class diJ_plus_link;
+
+  // Set up the tiles:
+  void _initialise_tiles(const std::vector<PseudoJet> & particles);
+
+  // return the full distance of a particle to its NN
+  inline double _compute_diJ(const TiledJet * const jet) const {
+    double mom_fact = jet->momentum_factor();
+    if (jet->NN != NULL) {
+      double other_mom_fact = jet->NN->momentum_factor();
+      if (other_mom_fact < mom_fact) {mom_fact = other_mom_fact;}
+    }
+    return jet->NN_dist * mom_fact;
+  }
+
+  // reasonably robust return of tile index given irap and iphi, in particular
+  // it works even if iphi is negative
+  inline int _tile_index (int irap, int iphi) const {
+    // note that (-1)%n = -1 so that we have to add _n_tiles_phi
+    // before performing modulo operation
+    return (irap-_tiles_irap_min)*_n_tiles_phi
+                  + (iphi+_n_tiles_phi) % _n_tiles_phi;
+  }
+
+  int  _tile_index(const double rap, const double phi) const;
+  void _tiledjet_set_jetinfo ( TiledJet * const tiled_jet, const PseudoJet &jet, int index);
+  void _bj_remove_from_tiles(TiledJet * const jet);
+  void _initialise_tiles();
+  void _print_tiles(TiledJet * briefjets ) const;
+  void _add_neighbours_to_tile_union(const int tile_index, int & n_near_tiles) const;
+  void _add_untagged_neighbours_to_tile_union(const int tile_index, int & n_near_tiles);
+
+
+  /// contains the briefjets
+  TiledJet * briefjets;
+
+  /// semaphores for the current extent of our structure
+  TiledJet * head;
+
+  /// currently active number of jets
+  int n;
+
+  /// where_is[i] contains a pointer to the jet with index i
+  std::vector<TiledJet *> where_is;
+
+  /// helper to keep tracks of tiles to be checked for updates
+  std::vector<int> tile_union;
+
+  /// a table containing all the (full) distances to each object's NN
+  diJ_plus_link * diJ;
+
+  /// tiling information
+  std::vector<Tile> _tiles;
+  double _requested_tile_size;
+  double _tiles_rap_min, _tiles_rap_max;
+  double _tile_size_rap, _tile_size_phi;
+  int    _n_tiles_phi,_tiles_irap_min,_tiles_irap_max;
+  
+  /// a class that wraps around the BJ, supplementing it with extra information
+  /// such as pointers to neighbours, etc.
+  class TiledJet : public BJ {
+  public:
+    void init(const PseudoJet & jet, int index_in) {
+      BJ::init(jet);
+      other_init(index_in);
+    }
+    void init(const PseudoJet & jet, int index_in, I * info) {
+      BJ::init(jet, info);
+      other_init(index_in);
+    }
+    void other_init(int index_in) {
+      _index = index_in;
+      NN_dist = BJ::geometrical_beam_distance();
+      NN = NULL;
+    }
+    int jet_index() const {return _index;}
+    
+    double NN_dist;
+    TiledJet * NN, *previous, * next;
+    int tile_index, diJ_posn;
+    // routines that are useful in the minheap version of tiled
+    // clustering ("misuse" the otherwise unused diJ_posn, so as
+    // to indicate whether jets need to have their minheap entries
+    // updated).
+    inline void label_minheap_update_needed() {diJ_posn = 1;}
+    inline void label_minheap_update_done()   {diJ_posn = 0;}
+    inline bool minheap_update_needed() const {return diJ_posn==1;}
+    
+  private:
+    int _index;
+  };
+
+  /// number of neighbours that a tile will have (rectangular geometry
+  /// gives 9 neighbours).
+  static const int n_tile_neighbours = 9;
+  //----------------------------------------------------------------------
+  /// The fundamental structures to be used for the tiled N^2 algorithm
+  /// (see CCN27-44 for some discussion of pattern of tiling)
+  class Tile {
+  public:
+    /// pointers to neighbouring tiles, including self
+    Tile *   begin_tiles[n_tile_neighbours]; 
+    /// neighbouring tiles, excluding self
+    Tile **  surrounding_tiles; 
+    /// half of neighbouring tiles, no self
+    Tile **  RH_tiles;  
+    /// just beyond end of tiles
+    Tile **  end_tiles; 
+    /// start of list of BriefJets contained in this tile
+    TiledJet * head;    
+    /// sometimes useful to be able to tag a tile
+    bool     tagged;    
+  };
+
+  // structure that holds the real, full, distance (as well as a pointer to the corresponding TiledJet)
+  struct diJ_plus_link {
+    double     diJ; // the distance
+    TiledJet * jet; // the jet (i) for which we've found this distance
+                    // (whose NN will the J).
+  };
+
+};
+
+
+
+//----------------------------------------------------------------------
+template<class BJ, class I> void NNFJN2Tiled<BJ,I>::start(const std::vector<PseudoJet> & jets) {
+
+  _initialise_tiles(jets);
+
+  n = jets.size();
+
+  briefjets = new TiledJet[n];
+  where_is.resize(2*n);
+
+  TiledJet * jetA = briefjets, * jetB;
+
+  // will be used quite deep inside loops, but declare it here so that
+  // memory (de)allocation gets done only once
+  tile_union.resize(3*n_tile_neighbours);
+  
+  // initialise the basic jet info 
+  for (int i = 0; i< n; i++) {
+    _tiledjet_set_jetinfo(jetA, jets[i], i);
+    where_is[i] = jetA;
+    jetA++; // move on to next entry of briefjets
+  }
+
+  head = briefjets; // a nicer way of naming start
+
+  // set up the initial nearest neighbour information
+  typename std::vector<Tile>::const_iterator tile;
+  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
+    // first do it on this tile
+    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
+      for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
+	double dist = jetA->geometrical_distance(jetB);
+	if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
+	if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
+      }
+    }
+    // then do it for RH tiles
+    for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
+      for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
+	for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
+	  double dist = jetA->geometrical_distance(jetB);
+	  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
+	  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
+	}
+      }
+    }
+    // no need to do it for LH tiles, since they are implicitly done
+    // when we set NN for both jetA and jetB on the RH tiles.
+  }
+  
+  diJ = new diJ_plus_link[n];
+  jetA = head;
+  for (int i = 0; i < n; i++) {
+    diJ[i].diJ = _compute_diJ(jetA); // kt distance * R^2
+    diJ[i].jet = jetA;  // our compact diJ table will not be in	     
+    jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
+                        // so set up bi-directional correspondence here.
+    jetA++; // have jetA follow i 
+  }
+}
+
+
+//----------------------------------------------------------------------
+template<class BJ, class I> double NNFJN2Tiled<BJ,I>::dij_min(int & iA, int & iB) {
+  // find the minimum of the diJ on this round
+  diJ_plus_link * best, *stop; // pointers a bit faster than indices
+                               // could use best to keep track of diJ
+                               // min, but it turns out to be
+                               // marginally faster to have a separate
+                               // variable (avoids n dereferences at
+                               // the expense of n/2 assignments).
+  double diJ_min = diJ[0].diJ; // initialise the best one here.
+  best = diJ;                  // and here
+  stop = diJ+n;
+  for (diJ_plus_link * here = diJ+1; here != stop; here++) {
+    if (here->diJ < diJ_min) {best = here; diJ_min  = here->diJ;}
+  }
+
+  // return information to user about recombination
+  TiledJet * jetA = best->jet;
+  iA = jetA->jet_index();
+  iB = jetA->NN ? jetA->NN->jet_index() : -1;
+  return diJ_min;
+}
+
+
+//----------------------------------------------------------------------
+// remove jetA from the list
+template<class BJ, class I> void NNFJN2Tiled<BJ,I>::remove_jet(int iA) {
+  TiledJet * jetA = where_is[iA];
+
+  _bj_remove_from_tiles(jetA);
+
+  // first establish the set of tiles over which we are going to
+  // have to run searches for updated and new nearest-neighbours --
+  // basically a combination of vicinity of the tiles of the two old
+  // and one new jet.
+  int n_near_tiles = 0;
+  _add_untagged_neighbours_to_tile_union(jetA->tile_index, n_near_tiles);
+
+  // now update our nearest neighbour info and diJ table
+  // first reduce size of table
+  n--;
+  // then compactify the diJ by taking the last of the diJ and copying
+  // it to the position occupied by the diJ for jetA
+  diJ[n].jet->diJ_posn = jetA->diJ_posn;
+  diJ[jetA->diJ_posn] = diJ[n];
+
+  // updating other particles' NN.
+  // Run over all tiles in our union 
+  for (int itile = 0; itile < n_near_tiles; itile++) {
+    Tile * tile_ptr = &_tiles[tile_union[itile]];
+    tile_ptr->tagged = false; // reset tag, since we're done with unions
+    // run over all jets in the current tile
+    for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
+      // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
+      if (jetI->NN == jetA) {
+        jetI->NN_dist = jetI->geometrical_beam_distance();
+        jetI->NN      = NULL;
+        // now go over tiles that are neighbours of I (include own tile)
+        for (Tile ** near_tile  = tile_ptr->begin_tiles; 
+             near_tile != tile_ptr->end_tiles; near_tile++) {
+          // and then over the contents of that tile
+          for (TiledJet * jetJ  = (*near_tile)->head; jetJ != NULL; jetJ = jetJ->next) {
+            double dist = jetI->geometrical_distance(jetJ);
+            if (dist < jetI->NN_dist && jetJ != jetI) {
+              jetI->NN_dist = dist; jetI->NN = jetJ;
+            }
+          }
+        }
+        diJ[jetI->diJ_posn].diJ = _compute_diJ(jetI); // update diJ kt-dist
+      }
+    }
+  }
+
+}
+
+
+//----------------------------------------------------------------------
+template<class BJ, class I> void NNFJN2Tiled<BJ,I>::merge_jets(int iA, int iB, 
+					const PseudoJet & jet, int index) {
+
+  TiledJet * jetA = where_is[iA];
+  TiledJet * jetB = where_is[iB];
+
+  // jet-jet recombination
+  // If necessary relabel A & B to ensure jetB < jetA, that way if
+  // the larger of them == newtail then that ends up being jetA and 
+  // the new jet that is added as jetB is inserted in a position that
+  // has a future!
+  if (jetA < jetB) {std::swap(jetA,jetB);}
+
+  // what was jetB will now become the new jet
+  _bj_remove_from_tiles(jetA);
+  TiledJet oldB = * jetB;  // take a copy because we will need it...
+  _bj_remove_from_tiles(jetB);
+  _tiledjet_set_jetinfo(jetB, jet, index); // cause jetB to become _jets[nn]
+                                           // (also registers the jet in the tiling)
+  where_is[index] = jetB;
+  
+  // first establish the set of tiles over which we are going to
+  // have to run searches for updated and new nearest-neighbours --
+  // basically a combination of vicinity of the tiles of the two old
+  // and one new jet.
+  int n_near_tiles = 0;
+  _add_untagged_neighbours_to_tile_union(jetA->tile_index, n_near_tiles);
+  if (jetB->tile_index != jetA->tile_index) {
+    _add_untagged_neighbours_to_tile_union(jetB->tile_index, n_near_tiles);
+  }
+  if (oldB.tile_index != jetA->tile_index && 
+      oldB.tile_index != jetB->tile_index) {
+    _add_untagged_neighbours_to_tile_union(oldB.tile_index, n_near_tiles);
+  }
+
+  // now update our nearest neighbour info and diJ table
+  // first reduce size of table
+  n--;
+  // then compactify the diJ by taking the last of the diJ and copying
+  // it to the position occupied by the diJ for jetA
+  diJ[n].jet->diJ_posn = jetA->diJ_posn;
+  diJ[jetA->diJ_posn] = diJ[n];
+
+  // Initialise jetB's NN distance as well as updating it for 
+  // other particles.
+  // Run over all tiles in our union 
+  for (int itile = 0; itile < n_near_tiles; itile++) {
+    Tile * tile_ptr = &_tiles[tile_union[itile]];
+    tile_ptr->tagged = false; // reset tag, since we're done with unions
+    // run over all jets in the current tile
+    for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
+      // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
+      if ((jetI->NN == jetA) || (jetI->NN == jetB)) {
+        jetI->NN_dist = jetI->geometrical_beam_distance();
+        jetI->NN      = NULL;
+        // now go over tiles that are neighbours of I (include own tile)
+        for (Tile ** near_tile  = tile_ptr->begin_tiles; near_tile != tile_ptr->end_tiles; near_tile++) {
+          // and then over the contents of that tile
+          for (TiledJet * jetJ  = (*near_tile)->head; jetJ != NULL; jetJ = jetJ->next) {
+            double dist = jetI->geometrical_distance(jetJ);
+            if (dist < jetI->NN_dist && jetJ != jetI) {
+              jetI->NN_dist = dist; jetI->NN = jetJ;
+            }
+          }
+        }
+        diJ[jetI->diJ_posn].diJ = _compute_diJ(jetI); // update diJ kt-dist
+      }
+      // check whether new jetB is closer than jetI's current NN and
+      // if jetI is closer than jetB's current (evolving) nearest
+      // neighbour. Where relevant update things
+      double dist = jetI->geometrical_distance(jetB);
+      if (dist < jetI->NN_dist) {
+        if (jetI != jetB) {
+          jetI->NN_dist = dist;
+          jetI->NN = jetB;
+          diJ[jetI->diJ_posn].diJ = _compute_diJ(jetI); // update diJ...
+        }
+      }
+      if (dist < jetB->NN_dist) {
+        if (jetI != jetB) {
+          jetB->NN_dist = dist;
+          jetB->NN      = jetI;}
+      }
+    }
+  }
+
+  // finally, register the updated kt distance for B
+  diJ[jetB->diJ_posn].diJ = _compute_diJ(jetB);
+}
+
+
+//----------------------------------------------------------------------
+/// Set up the tiles:
+///  - decide the range in eta
+///  - allocate the tiles
+///  - set up the cross-referencing info between tiles
+///
+/// The neighbourhood of a tile is set up as follows
+///
+/// 	      LRR
+///           LXR
+///           LLR
+///
+/// such that tiles is an array containing XLLLLRRRR with pointers
+///                                         |   \ RH_tiles
+///                                         \ surrounding_tiles
+///
+/// with appropriate precautions when close to the edge of the tiled
+/// region.
+///
+template <class BJ, class I>
+void NNFJN2Tiled<BJ,I>::_initialise_tiles(const std::vector<PseudoJet> &particles) {
+
+  // first decide tile sizes (with a lower bound to avoid huge memory use with
+  // very small R)
+  double default_size = _requested_tile_size>0.1 ? _requested_tile_size : 0.1;
+  _tile_size_rap = default_size;
+  // it makes no sense to go below 3 tiles in phi -- 3 tiles is
+  // sufficient to make sure all pair-wise combinations up to pi in
+  // phi are possible
+  _n_tiles_phi   = int(floor(twopi/default_size));
+  if (_n_tiles_phi<3) _n_tiles_phi = 3;
+  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
+
+  TilingExtent tiling_analysis(particles);
+  _tiles_rap_min = tiling_analysis.minrap();
+  _tiles_rap_max = tiling_analysis.maxrap();
+
+  // now adjust the values
+  _tiles_irap_min = int(floor(_tiles_rap_min/_tile_size_rap));
+  _tiles_irap_max = int(floor( _tiles_rap_max/_tile_size_rap));
+  _tiles_rap_min = _tiles_irap_min * _tile_size_rap;
+  _tiles_rap_max = _tiles_irap_max * _tile_size_rap;
+
+  // allocate the tiles
+  _tiles.resize((_tiles_irap_max-_tiles_irap_min+1)*_n_tiles_phi);
+
+  // now set up the cross-referencing between tiles
+  for (int irap = _tiles_irap_min; irap <= _tiles_irap_max; irap++) {
+    for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
+      Tile * tile = & _tiles[_tile_index(irap,iphi)];
+      // no jets in this tile yet
+      tile->head = NULL; // first element of tiles points to itself
+      tile->begin_tiles[0] =  tile;
+      Tile ** pptile = & (tile->begin_tiles[0]);
+      pptile++;
+      //
+      // set up L's in column to the left of X
+      tile->surrounding_tiles = pptile;
+      if (irap > _tiles_irap_min) {
+	// with the itile subroutine, we can safely run tiles from
+	// idphi=-1 to idphi=+1, because it takes care of
+	// negative and positive boundaries
+	for (int idphi = -1; idphi <=+1; idphi++) {
+	  *pptile = & _tiles[_tile_index(irap-1,iphi+idphi)];
+	  pptile++;
+	}	
+      }
+      // now set up last L (below X)
+      *pptile = & _tiles[_tile_index(irap,iphi-1)];
+      pptile++;
+      // set up first R (above X)
+      tile->RH_tiles = pptile;
+      *pptile = & _tiles[_tile_index(irap,iphi+1)];
+      pptile++;
+      // set up remaining R's, to the right of X
+      if (irap < _tiles_irap_max) {
+	for (int idphi = -1; idphi <= +1; idphi++) {
+	  *pptile = & _tiles[_tile_index(irap+1,iphi+idphi)];
+	  pptile++;
+	}	
+      }
+      // now put semaphore for end tile
+      tile->end_tiles = pptile;
+      // finally make sure tiles are untagged
+      tile->tagged = false;
+    }
+  }
+
+}
+
+//----------------------------------------------------------------------
+/// return the tile index corresponding to the given rap,phi point
+template <class BJ, class I>
+int NNFJN2Tiled<BJ,I>::_tile_index(const double rap, const double phi) const {
+  int irap, iphi;
+  if      (rap <= _tiles_rap_min) {irap = 0;}
+  else if (rap >= _tiles_rap_max) {irap = _tiles_irap_max-_tiles_irap_min;}
+  else {
+    //irap = int(floor((rap - _tiles_rap_min) / _tile_size_rap));
+    irap = int(((rap - _tiles_rap_min) / _tile_size_rap));
+    // following needed in case of rare but nasty rounding errors
+    if (irap > _tiles_irap_max-_tiles_irap_min) {
+      irap = _tiles_irap_max-_tiles_irap_min;} 
+  }
+  // allow for some extent of being beyond range in calculation of phi
+  // as well
+  //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
+  // with just int and no floor, things run faster but beware
+  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
+  return (iphi + irap * _n_tiles_phi);
+}
+
+//----------------------------------------------------------------------
+template <class BJ, class I>
+void NNFJN2Tiled<BJ,I>::_bj_remove_from_tiles(TiledJet * const jet) {
+  Tile * tile = & _tiles[jet->tile_index];
+
+  if (jet->previous == NULL) {
+    // we are at head of the tile, so reset it.
+    // If this was the only jet on the tile then tile->head will now be NULL
+    tile->head = jet->next;
+  } else {
+    // adjust link from previous jet in this tile
+    jet->previous->next = jet->next;
+  }
+  if (jet->next != NULL) {
+    // adjust backwards-link from next jet in this tile
+    jet->next->previous = jet->previous;
+  }
+}
+
+
+//----------------------------------------------------------------------
+// overloaded version which additionally sets up information regarding the
+// tiling
+template <class BJ, class I>
+inline void NNFJN2Tiled<BJ,I>::_tiledjet_set_jetinfo(TiledJet * const tile_jet,
+                                                   const PseudoJet &jet, 
+                                                   int index) {
+  // the this-> in the next line is required by standard compiler
+  // see e.g. http://stackoverflow.com/questions/10639053/name-lookups-in-c-templates
+  this->init_jet(tile_jet, jet, index); 
+
+  // Then do the setup specific to the tiled case.
+
+  // Find out which tile it belonds to
+  tile_jet->tile_index = _tile_index(tile_jet->rap(), tile_jet->phi());
+
+  // Insert it into the tile's linked list of jets
+  Tile * tile = &_tiles[tile_jet->tile_index];
+  tile_jet->previous   = NULL;
+  tile_jet->next       = tile->head;
+  if (tile_jet->next != NULL) {tile_jet->next->previous = tile_jet;}
+  tile->head      = tile_jet;
+}
+
+//----------------------------------------------------------------------
+/// Add to the vector tile_union the tiles that are in the neighbourhood
+/// of the specified tile_index, including itself -- start adding
+/// from position n_near_tiles-1, and increase n_near_tiles as
+/// you go along (could have done it more C++ like with vector with reserved
+/// space, but fear is that it would have been slower, e.g. checking
+/// for end of vector at each stage to decide whether to resize it)
+template <class BJ, class I>
+void NNFJN2Tiled<BJ,I>::_add_neighbours_to_tile_union(const int tile_index, 
+                                                    int & n_near_tiles) const {
+  for (Tile * const * near_tile = _tiles[tile_index].begin_tiles; 
+       near_tile != _tiles[tile_index].end_tiles; near_tile++){
+    // get the tile number
+    tile_union[n_near_tiles] = *near_tile - & _tiles[0];
+    n_near_tiles++;
+  }
+}
+
+//----------------------------------------------------------------------
+/// Like _add_neighbours_to_tile_union, but only adds neighbours if 
+/// their "tagged" status is false; when a neighbour is added its
+/// tagged status is set to true.
+///
+/// Note that with a high level of warnings (-pedantic -Wextra -ansi,
+/// gcc complains about tile_index maybe being used uninitialised for
+/// oldB in ClusterSequence::_minheap_faster_tiled_N2_cluster(). We
+/// have explicitly checked that it was harmless so we could disable
+/// the gcc warning by hand using the construct below
+///
+///  #pragma GCC diagnostic push
+///  #pragma GCC diagnostic ignored "-Wpragmas"
+///  #pragma GCC diagnostic ignored "-Wuninitialized"
+///  #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
+///    ...
+///  #pragma GCC diagnostic pop
+///
+/// the @GCC diagnostic push/pop directive was only introduced in
+/// gcc-4.6, so for broader usage, we'd need to insert #pragma GCC
+/// diagnostic ignored "-Wpragmas" at the top of this file
+template <class BJ, class I>
+inline void NNFJN2Tiled<BJ,I>::_add_untagged_neighbours_to_tile_union(
+               const int tile_index, 
+	       int & n_near_tiles)  {
+  for (Tile ** near_tile = _tiles[tile_index].begin_tiles; 
+       near_tile != _tiles[tile_index].end_tiles; near_tile++){
+    if (! (*near_tile)->tagged) {
+      (*near_tile)->tagged = true;
+      // get the tile number
+      tile_union[n_near_tiles] = *near_tile - & _tiles[0];
+      n_near_tiles++;
+    }
+  }
+}
+
+
+
+FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
+
+
+#endif // __FASTJET_NNFJN2TILED_HH__
Index: external/fastjet/NNH.hh
===================================================================
--- external/fastjet/NNH.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/NNH.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: NNH.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: NNH.hh 4034 2016-03-02 00:20:27Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -32,35 +32,7 @@
 //FJENDHEADER
 
-#include<fastjet/ClusterSequence.hh>
-
+#include <fastjet/NNBase.hh>
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
-
-/// @ingroup advanced_usage
-/// \class _NoInfo
-/// dummy class, used as a default template argument
-class _NoInfo {};
-
-/// @ingroup advanced_usage
-/// \class NNHInfo
-/// template that will help initialise a BJ with a PseudoJet and extra information
-template<class I> class NNHInfo {
-public:
-  NNHInfo()         : _info(NULL) {}
-  NNHInfo(I * info) : _info(info) {}
-  template<class NNBJ> void init_jet(NNBJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index, _info);}
-private:
-  I * _info;
-};
-
-/// @ingroup advanced_usage
-/// Specialisation of NNHInfo for cases where there is no extra info
-template<> class NNHInfo<_NoInfo>  {
-public:
-  NNHInfo()           {}
-  NNHInfo(_NoInfo * ) {}
-  template<class NNBJ> void init_jet(NNBJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index);}
-};
-
 
 //----------------------------------------------------------------------
@@ -68,17 +40,11 @@
 /// \class NNH
 /// Help solve closest pair problems with generic interparticle and
-/// beam distance.
-///
-/// Class to help solve closest pair problems with generic interparticle
-/// distances and a beam distance, using Anderberg's Nearest Neighbour
-/// Heuristic.
-///
-/// It is templated with a BJ (brief jet) class --- BJ should
-/// basically cache the minimal amount of information that is needed
-/// to efficiently calculate interparticle distances and particle-beam
-/// distances.
-///
-/// This class can be used with or without an extra "Information" template, 
-/// i.e. NNB<BJ> or NNH<BJ,I>
+/// beam distance (generic case)
+///
+/// (see NNBase.hh for an introductory description)
+///
+/// This variant provides an implementation for any distance measure.
+/// It is templated with a BJ (brief jet) classand can be used with or
+/// without an extra "Information" template, i.e. NNH<BJ> or NNH<BJ,I>
 /// 
 /// For the NNH<BJ> version of the class to function, BJ must provide 
@@ -88,5 +54,5 @@
 ///  - double BJ::distance(const BJ * other_bj_jet); // distance between this and other_bj_jet
 ///  - double BJ::beam_distance()                  ; // distance to the beam
-///
+/// 
 /// For the NNH<BJ,I> version to function, the BJ::init(...) member
 /// must accept an extra argument
@@ -94,6 +60,6 @@
 ///  - void   BJ::init(const PseudoJet & jet, I * info);   // initialise with a PseudoJet + info
 ///
-/// where info might be a pointer to a class that contains, e.g., information
-/// about R, or other parameters of the jet algorithm 
+/// NOTE: THE DISTANCE MUST BE SYMMETRIC I.E. SATISFY
+///     a.distance(b) == b.distance(a)
 ///
 /// For an example of how the NNH<BJ> class is used, see the Jade (and
@@ -107,16 +73,13 @@
 /// implementations.
 ///
-///
-/// Implementation note: this class derives from NNHInfo, which deals
-/// with storing any global information that is needed during the clustering
-
-template<class BJ, class I = _NoInfo> class NNH : public NNHInfo<I> {
+template<class BJ, class I = _NoInfo> class NNH : public NNBase<I> {
 public:
 
   /// constructor with an initial set of jets (which will be assigned indices
   /// 0 ... jets.size()-1
-  NNH(const std::vector<PseudoJet> & jets) {start(jets);}
-  NNH(const std::vector<PseudoJet> & jets, I * info) : NNHInfo<I>(info) {start(jets);}
-  
+  NNH(const std::vector<PseudoJet> & jets)           : NNBase<I>()     {start(jets);}
+  NNH(const std::vector<PseudoJet> & jets, I * info) : NNBase<I>(info) {start(jets);}
+
+  // initialisation from a given list of particles
   void start(const std::vector<PseudoJet> & jets);
 
Index: external/fastjet/PseudoJet.cc
===================================================================
--- external/fastjet/PseudoJet.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/PseudoJet.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: PseudoJet.cc 3652 2014-09-03 13:31:13Z salam $
+// $Id: PseudoJet.cc 4100 2016-03-15 20:50:22Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -409,9 +409,9 @@
 string PseudoJet::description() const{
   // the "default" case of a PJ which does not belong to any cluster sequence
-  if (!_structure())
+  if (!_structure)
     return "standard PseudoJet (with no associated clustering information)";
   
   // for all the other cases, the description comes from the structure
-  return _structure()->description();
+  return _structure->description();
 }
 
@@ -429,5 +429,5 @@
 // ClusterSequence
 bool PseudoJet::has_associated_cluster_sequence() const{
-  return (_structure()) && (_structure->has_associated_cluster_sequence());
+  return (_structure) && (_structure->has_associated_cluster_sequence());
 }
 
@@ -446,5 +446,5 @@
 // ClusterSequence that is still valid
 bool PseudoJet::has_valid_cluster_sequence() const{
-  return (_structure()) && (_structure->has_valid_cluster_sequence());
+  return (_structure) && (_structure->has_valid_cluster_sequence());
 }
 
@@ -467,7 +467,7 @@
 
 //----------------------------------------------------------------------
-// return true if there is some strusture associated with this PseudoJet
+// return true if there is some structure associated with this PseudoJet
 bool PseudoJet::has_structure() const{
-  return _structure();
+  return bool(_structure);
 }
 
@@ -478,6 +478,6 @@
 // return NULL if there is no associated structure
 const PseudoJetStructureBase* PseudoJet::structure_ptr() const {
-  if (!_structure()) return NULL;
-  return _structure();
+  //if (!_structure) return NULL;
+  return _structure.get();
 }
   
@@ -493,6 +493,6 @@
 // underlying structure.
 PseudoJetStructureBase* PseudoJet::structure_non_const_ptr(){
-  if (!_structure()) return NULL;
-  return _structure();
+  //if (!_structure) return NULL;
+  return _structure.get();
 }
   
@@ -503,7 +503,7 @@
 // throw an error if there is no associated structure
 const PseudoJetStructureBase* PseudoJet::validated_structure_ptr() const {
-  if (!_structure()) 
+  if (!_structure) 
     throw Error("Trying to access the structure of a PseudoJet which has no associated structure");
-  return _structure();
+  return _structure.get();
 }
   
@@ -573,5 +573,5 @@
 // returns true if the PseudoJet has constituents
 bool PseudoJet::has_constituents() const{
-  return (_structure()) && (_structure->has_constituents());
+  return (_structure) && (_structure->has_constituents());
 }
 
@@ -586,5 +586,5 @@
 // returns true if the PseudoJet has support for exclusive subjets
 bool PseudoJet::has_exclusive_subjets() const{
-  return (_structure()) && (_structure->has_exclusive_subjets());
+  return (_structure) && (_structure->has_exclusive_subjets());
 }
 
@@ -670,5 +670,5 @@
 // ClusterSequence have no pieces and this methos will return false.
 bool PseudoJet::has_pieces() const{
-  return ((_structure()) && (_structure->has_pieces(*this)));
+  return ((_structure) && (_structure->has_pieces(*this)));
 }
 
@@ -766,33 +766,4 @@
 
 
-
-//----------------------------------------------------------------------
-/// given a vector of values with a one-to-one correspondence with the
-/// vector of objects, sort objects into an order such that the
-/// associated values would be in increasing order
-template<class T> vector<T>  objects_sorted_by_values(
-                       const vector<T> & objects, 
-		       const vector<double> & values) {
-
-  assert(objects.size() == values.size());
-
-  // get a vector of indices
-  vector<int> indices(values.size());
-  for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
-  
-  // sort the indices
-  sort_indices(indices, values);
-  
-  // copy the objects 
-  vector<T> objects_sorted(objects.size());
-  
-  // place the objects in the correct order
-  for (size_t i = 0; i < indices.size(); i++) {
-    objects_sorted[i] = objects[indices[i]];
-  }
-
-  return objects_sorted;
-}
-
 //----------------------------------------------------------------------
 /// return a vector of jets sorted into decreasing kt2
Index: external/fastjet/PseudoJet.hh
===================================================================
--- external/fastjet/PseudoJet.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/PseudoJet.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: PseudoJet.hh 3566 2014-08-11 15:36:34Z salam $
+// $Id: PseudoJet.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -437,5 +437,7 @@
   /// retrieve a pointer to the (const) user information
   const UserInfoBase * user_info_ptr() const{
-    if (!_user_info()) return NULL;
+    // the line below is not needed since the next line would anyway
+    // return NULL in that case
+    //if (!_user_info) return NULL;
     return _user_info.get();
   }
@@ -841,4 +843,5 @@
 inline bool operator!=( const double val, const PseudoJet & a) {return !(a==val);}
 
+/// returns the 4-vector dot product of a and b
 inline double dot_product(const PseudoJet & a, const PseudoJet & b) {
   return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
@@ -880,5 +883,27 @@
 /// touch the values vector in the process).
 template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects, 
-					      const std::vector<double> & values);
+					      const std::vector<double> & values) {
+  //assert(objects.size() == values.size());
+  if (objects.size() != values.size()){
+    throw Error("fastjet::objects_sorted_by_values(...): the size of the 'objects' vector must match the size of the 'values' vector");
+  }
+  
+  // get a vector of indices
+  std::vector<int> indices(values.size());
+  for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
+  
+  // sort the indices
+  sort_indices(indices, values);
+  
+  // copy the objects 
+  std::vector<T> objects_sorted(objects.size());
+  
+  // place the objects in the correct order
+  for (size_t i = 0; i < indices.size(); i++) {
+    objects_sorted[i] = objects[indices[i]];
+  }
+
+  return objects_sorted;
+}
 
 /// \if internal_doc
@@ -970,5 +995,5 @@
 template<typename TransformerType>
 bool PseudoJet::has_structure_of() const{
-  if (!_structure()) return false;
+  if (!_structure) return false;
 
   return dynamic_cast<const typename TransformerType::StructureType *>(_structure.get()) != 0;
@@ -980,5 +1005,5 @@
 template<typename TransformerType>
 const typename TransformerType::StructureType & PseudoJet::structure_of() const{
-  if (!_structure()) 
+  if (!_structure) 
     throw Error("Trying to access the structure of a PseudoJet without an associated structure");
 
Index: external/fastjet/RangeDefinition.hh
===================================================================
--- external/fastjet/RangeDefinition.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/RangeDefinition.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: RangeDefinition.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: RangeDefinition.hh 4074 2016-03-08 09:09:25Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -35,4 +35,5 @@
 #include "fastjet/Error.hh"
 #include "fastjet/LimitedWarning.hh"
+#include "fastjet/internal/deprecated.hh"
 #include<sstream>
 #include<iostream>
@@ -51,7 +52,9 @@
 public:
   /// default constructor
+  FASTJET_DEPRECATED_MSG("RangeDefinition is deprecated since FastJet 3.0. Use the Selector mechanism instead")
   RangeDefinition() { _warn_deprecated(); }
 
   /// constructor for a range definition given by |y|<rapmax
+  FASTJET_DEPRECATED_MSG("RangeDefinition is deprecated since FastJet 3.0. Use the Selector mechanism instead")
   RangeDefinition(double rapmax) {  _warn_deprecated(); 
                      assert ( rapmax > 0.0 );
Index: external/fastjet/RectangularGrid.cc
===================================================================
--- external/fastjet/RectangularGrid.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/RectangularGrid.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: RectangularGrid.cc 3771 2014-12-22 21:13:22Z salam $
+// $Id: RectangularGrid.cc 3773 2014-12-22 22:44:46Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
Index: external/fastjet/RectangularGrid.hh
===================================================================
--- external/fastjet/RectangularGrid.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/RectangularGrid.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: RectangularGrid.hh 3771 2014-12-22 21:13:22Z salam $
+// $Id: RectangularGrid.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -118,24 +118,28 @@
   RectangularGrid();
 
-  virtual int n_tiles() const {return _ntotal;}
+  virtual int n_tiles() const FASTJET_OVERRIDE {return _ntotal;}
 
-  virtual int n_good_tiles() const {return _ngood;}
+  virtual int n_good_tiles() const FASTJET_OVERRIDE {return _ngood;}
 
   // this was being kept inline, but it seems to make little
   // difference whether it is or not (at least on Gavin's mac)
-  virtual int tile_index(const PseudoJet & p) const;
+  virtual int tile_index(const PseudoJet & p) const FASTJET_OVERRIDE;
 
   /// returns whether a given tile is good
   // tested in "issue" 2014-08-08-testing-rect-grid
-  virtual bool tile_is_good(int itile) const {return _tile_selector.worker() ? _is_good[itile] : true;}
+  virtual bool tile_is_good(int itile) const FASTJET_OVERRIDE {
+    return _tile_selector.worker() ? _is_good[itile] : true;
+  }
 
   /// returns the area of tile itile.
-  virtual double tile_area(int /* itile */) const {return mean_tile_area();}
+  virtual double tile_area(int /* itile */) const FASTJET_OVERRIDE {
+    return mean_tile_area();
+  }
 
   /// returns the mean area of tiles.
-  virtual double mean_tile_area() const {return _dphi*_dy;};
+  virtual double mean_tile_area() const FASTJET_OVERRIDE {return _dphi*_dy;};
 
   /// returns a textual description of the grid
-  virtual std::string description() const;
+  virtual std::string description() const FASTJET_OVERRIDE;
   
   /// returns the minimum rapidity extent of the grid
@@ -149,5 +153,5 @@
 
   /// returns true if the grid is in a suitably initialised state
-  virtual bool is_initialised() const {return _ntotal > 0;}
+  virtual bool is_initialised() const FASTJET_OVERRIDE {return _ntotal > 0;}
 
 private:
Index: external/fastjet/SharedPtr.hh
===================================================================
--- external/fastjet/SharedPtr.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/SharedPtr.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: SharedPtr.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: SharedPtr.hh 4051 2016-03-03 14:33:38Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -38,4 +38,6 @@
 // for our SharedPtr simply to be derived from the STL TR1 one.
 // #define __FASTJET_USETR1SHAREDPTR
+
+#include "fastjet/internal/deprecated.hh"
 
 #ifdef __FASTJET_USETR1SHAREDPTR
@@ -72,4 +74,7 @@
   SharedPtr(const SharedPtr<T> & t) : std::tr1::shared_ptr<T>(t) {}
   // for some reason operator() doesn't get inherited
+  #ifdef FASTJET_HAVE_EXPLICIT_FOR_OPERATORS
+  explicit
+  #endif
   inline operator bool() const {return (this->get()!=NULL);}
   /// return the pointer we're pointing to  
@@ -211,5 +216,9 @@
   }
   
-  /// return the pointer we're pointing to  
+  /// return the pointer we're pointing to
+  ///
+  /// Since FastJet 3.2.0, this is depracated since it is no longer
+  /// part of std::shared_ptr<T>. Use SharedPtr<T>::get() instead
+  FASTJET_DEPRECATED_MSG("Use SharedPtr<T>::get() instead")
   T* operator ()() const{
     if (_ptr==NULL) return NULL;
@@ -257,4 +266,7 @@
   /// conversion to bool
   /// This will allow you to use the indirection nicely
+  #ifdef FASTJET_HAVE_EXPLICIT_FOR_OPERATORS
+  explicit
+  #endif
   inline operator bool() const{
     return (get()!=NULL);
Index: external/fastjet/TilingExtent.cc
===================================================================
--- external/fastjet/TilingExtent.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/TilingExtent.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: TilingExtent.cc 3433 2014-07-23 08:17:03Z salam $
+// $Id: TilingExtent.cc 4034 2016-03-02 00:20:27Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -41,5 +41,9 @@
   _determine_rapidity_extent(cs.jets());
 }
-  
+
+TilingExtent::TilingExtent(const vector<PseudoJet> &particles) {
+  _determine_rapidity_extent(particles);
+}
+
 void TilingExtent::_determine_rapidity_extent(const vector<PseudoJet> & particles) {
   // have a binning of rapidity that goes from -nrap to nrap
Index: external/fastjet/VERSION
===================================================================
--- external/fastjet/VERSION	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/VERSION	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,1 +1,1 @@
-3.1.3
+3.2.1
Index: external/fastjet/WrappedStructure.hh
===================================================================
--- external/fastjet/WrappedStructure.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/WrappedStructure.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: WrappedStructure.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: WrappedStructure.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -57,5 +57,5 @@
   WrappedStructure(const SharedPtr<PseudoJetStructureBase> & to_be_shared)
     : _structure(to_be_shared){
-    if (!_structure())
+    if (!_structure)
       throw Error("Trying to construct a wrapped structure around an empty (NULL) structure");
   }
@@ -65,5 +65,5 @@
 
   /// description
-  virtual std::string description() const{ 
+  virtual std::string description() const FASTJET_OVERRIDE{ 
     return "PseudoJet wrapping the structure ("+_structure->description()+")"; 
   }
@@ -76,5 +76,5 @@
   //-------------------------------------------------------------
   /// returns true if there is an associated ClusterSequence
-  virtual bool has_associated_cluster_sequence() const {
+  virtual bool has_associated_cluster_sequence() const FASTJET_OVERRIDE {
     return _structure->has_associated_cluster_sequence();
   }
@@ -82,5 +82,5 @@
   /// get a (const) pointer to the parent ClusterSequence (NULL if
   /// inexistent)
-  virtual const ClusterSequence* associated_cluster_sequence() const{
+  virtual const ClusterSequence* associated_cluster_sequence() const FASTJET_OVERRIDE{
     return _structure->associated_cluster_sequence();
   }
@@ -88,5 +88,5 @@
   /// returns true if this PseudoJet has an associated and still
   /// valid ClusterSequence.
-  virtual bool has_valid_cluster_sequence() const {
+  virtual bool has_valid_cluster_sequence() const FASTJET_OVERRIDE {
     return _structure->has_valid_cluster_sequence();
   }
@@ -94,5 +94,5 @@
   /// if the jet has a valid associated cluster sequence then return a
   /// pointer to it; otherwise throw an error
-  virtual const ClusterSequence * validated_cs() const{
+  virtual const ClusterSequence * validated_cs() const FASTJET_OVERRIDE{
     return _structure->validated_cs();
   }
@@ -100,5 +100,5 @@
   /// if the jet has valid area information then return a pointer to
   /// the associated ClusterSequenceAreaBase object; otherwise throw an error
-  virtual const ClusterSequenceAreaBase * validated_csab() const{
+  virtual const ClusterSequenceAreaBase * validated_csab() const FASTJET_OVERRIDE{
     return _structure->validated_csab();
   }
@@ -120,5 +120,5 @@
   ///
   /// By default, throws an Error
-  virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const{
+  virtual bool has_partner(const PseudoJet &reference, PseudoJet &partner) const FASTJET_OVERRIDE{
     return _structure->has_partner(reference, partner);
   }
@@ -129,5 +129,5 @@
   /// 
   /// By default, throws an Error
-  virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const{
+  virtual bool has_child(const PseudoJet &reference, PseudoJet &child) const FASTJET_OVERRIDE{
     return _structure->has_child(reference, child);
   }
@@ -138,5 +138,6 @@
   ///
   /// By default, throws an Error
-  virtual bool has_parents(const PseudoJet &reference, PseudoJet &parent1, PseudoJet &parent2) const{
+  virtual bool has_parents(const PseudoJet &reference,
+                           PseudoJet &parent1, PseudoJet &parent2) const FASTJET_OVERRIDE{
     return _structure->has_parents(reference, parent1, parent2);
   }
@@ -146,5 +147,6 @@
   ///
   /// By default, throws an Error
-  virtual bool object_in_jet(const PseudoJet &reference, const PseudoJet &jet) const{
+  virtual bool object_in_jet(const PseudoJet &reference,
+                             const PseudoJet &jet) const FASTJET_OVERRIDE{
     return _structure->object_in_jet(reference, jet);
   }
@@ -154,5 +156,5 @@
   ///
   /// false by default
-  virtual bool has_constituents() const {
+  virtual bool has_constituents() const  FASTJET_OVERRIDE{
     return _structure->has_constituents();
   }
@@ -161,10 +163,10 @@
   ///
   /// By default, throws an Error
-  virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const{
+  virtual std::vector<PseudoJet> constituents(const PseudoJet &reference) const FASTJET_OVERRIDE{
     return _structure->constituents(reference);
   }
 
   /// return true if the structure supports exclusive_subjets. 
-  virtual bool has_exclusive_subjets() const {
+  virtual bool has_exclusive_subjets() const  FASTJET_OVERRIDE{
     return _structure->has_exclusive_subjets();
   }
@@ -180,5 +182,6 @@
   ///
   /// By default, throws an Error
-  virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference, const double & dcut) const{
+  virtual std::vector<PseudoJet> exclusive_subjets(const PseudoJet &reference,
+                                                   const double & dcut) const FASTJET_OVERRIDE{
     return _structure->exclusive_subjets(reference, dcut);
   }
@@ -189,5 +192,6 @@
   ///
   /// By default, throws an Error
-  virtual int n_exclusive_subjets(const PseudoJet &reference, const double & dcut) const{
+  virtual int n_exclusive_subjets(const PseudoJet &reference,
+                                  const double & dcut) const FASTJET_OVERRIDE{
     return _structure->n_exclusive_subjets(reference, dcut);
   }
@@ -198,5 +202,6 @@
   ///
   /// By default, throws an Error
-  virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference, int nsub) const{
+  virtual std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet &reference,
+                                                          int nsub) const FASTJET_OVERRIDE{
     return _structure->exclusive_subjets_up_to (reference, nsub);
   }
@@ -206,5 +211,5 @@
   ///
   /// By default, throws an Error
-  virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const{
+  virtual double exclusive_subdmerge(const PseudoJet &reference, int nsub) const FASTJET_OVERRIDE{
     return _structure->exclusive_subdmerge(reference, nsub);
   }
@@ -215,5 +220,5 @@
   ///
   /// By default, throws an Error
-  virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const{
+  virtual double exclusive_subdmerge_max(const PseudoJet &reference, int nsub) const FASTJET_OVERRIDE{
     return _structure->exclusive_subdmerge_max(reference, nsub);
   }
@@ -226,5 +231,5 @@
   ///
   /// false by default
-  virtual bool has_pieces(const PseudoJet &reference) const {
+  virtual bool has_pieces(const PseudoJet &reference) const FASTJET_OVERRIDE {
     return _structure->has_pieces(reference);
   }
@@ -233,5 +238,5 @@
   ///
   /// By default, throws an Error
-  virtual std::vector<PseudoJet> pieces(const PseudoJet &reference) const{
+  virtual std::vector<PseudoJet> pieces(const PseudoJet &reference) const FASTJET_OVERRIDE {
     return _structure->pieces(reference);
   }
@@ -244,5 +249,5 @@
   ///
   /// false by default
-  virtual bool has_area() const {
+  virtual bool has_area() const FASTJET_OVERRIDE  {
     return _structure->has_area();
   }
@@ -251,5 +256,5 @@
   ///
   /// By default, throws an Error
-  virtual double area(const PseudoJet &reference) const{
+  virtual double area(const PseudoJet &reference) const FASTJET_OVERRIDE{
     return _structure->area(reference);
   }
@@ -259,5 +264,5 @@
   ///
   /// By default, throws an Error
-  virtual double area_error(const PseudoJet &reference) const{
+  virtual double area_error(const PseudoJet &reference) const FASTJET_OVERRIDE{
     return _structure->area_error(reference);
   }
@@ -266,5 +271,5 @@
   ///
   /// By default, throws an Error
-  virtual PseudoJet area_4vector(const PseudoJet &reference) const{
+  virtual PseudoJet area_4vector(const PseudoJet &reference) const FASTJET_OVERRIDE{
     return _structure->area_4vector(reference);
   }
@@ -273,5 +278,5 @@
   ///
   /// By default, throws an Error
-  virtual bool is_pure_ghost(const PseudoJet &reference) const{
+  virtual bool is_pure_ghost(const PseudoJet &reference) const FASTJET_OVERRIDE{
     return _structure->is_pure_ghost(reference);
   }
Index: external/fastjet/config_auto.h
===================================================================
--- external/fastjet/config_auto.h	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/config_auto.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -64,6 +64,11 @@
 #endif
 
-/* defined if C++11 features have been enabled */
-/* #undef HAVE_CXX11_FEATURES */
+/* compile the deprecated parts of the interface using auto-ptr */
+#ifndef FASTJET_HAVE_AUTO_PTR_INTERFACE 
+#define FASTJET_HAVE_AUTO_PTR_INTERFACE  /**/ 
+#endif
+
+/* compiler supports c++14 deprecated keyword */
+/* #undef HAVE_CXX14_DEPRECATED */
 
 /* defined if demangling is enabled at configure time and is supported through
@@ -81,4 +86,12 @@
 #endif
 
+/* compiler supports the "explicit" keyword for operators */
+/* #undef HAVE_EXPLICIT_FOR_OPERATORS */
+
+/* compiler supports GNU c++ deprecated attribute */
+#ifndef FASTJET_HAVE_GNUCXX_DEPRECATED 
+#define FASTJET_HAVE_GNUCXX_DEPRECATED  /**/ 
+#endif
+
 /* Define to 1 if you have the <inttypes.h> header file. */
 #ifndef FASTJET_HAVE_INTTYPES_H 
@@ -96,4 +109,7 @@
 #endif
 
+/* compiler supports the "override" keyword */
+/* #undef HAVE_OVERRIDE */
+
 /* Define to 1 if you have the <stdint.h> header file. */
 #ifndef FASTJET_HAVE_STDINT_H 
@@ -131,6 +147,5 @@
 #endif
 
-/* Define to the sub-directory in which libtool stores uninstalled libraries.
-   */
+/* Define to the sub-directory where libtool stores uninstalled libraries. */
 #ifndef FASTJET_LT_OBJDIR 
 #define FASTJET_LT_OBJDIR  ".libs/" 
@@ -154,5 +169,5 @@
 /* Define to the full name and version of this package. */
 #ifndef FASTJET_PACKAGE_STRING 
-#define FASTJET_PACKAGE_STRING  "FastJet 3.1.3" 
+#define FASTJET_PACKAGE_STRING  "FastJet 3.2.1" 
 #endif
 
@@ -162,7 +177,12 @@
 #endif
 
+/* Define to the home page for this package. */
+#ifndef FASTJET_PACKAGE_URL 
+#define FASTJET_PACKAGE_URL  "" 
+#endif
+
 /* Define to the version of this package. */
 #ifndef FASTJET_PACKAGE_VERSION 
-#define FASTJET_PACKAGE_VERSION  "3.1.3" 
+#define FASTJET_PACKAGE_VERSION  "3.2.1" 
 #endif
 
@@ -174,5 +194,5 @@
 /* Version number of package */
 #ifndef FASTJET_VERSION 
-#define FASTJET_VERSION  "3.1.3" 
+#define FASTJET_VERSION  "3.2.1" 
 #endif
 
@@ -184,15 +204,15 @@
 /* Minor version of this package */
 #ifndef FASTJET_VERSION_MINOR 
-#define FASTJET_VERSION_MINOR  1 
+#define FASTJET_VERSION_MINOR  2 
 #endif
 
 /* Version of the package under the form XYYZZ (instead of X.Y.Z) */
 #ifndef FASTJET_VERSION_NUMBER 
-#define FASTJET_VERSION_NUMBER  30103 
+#define FASTJET_VERSION_NUMBER  30201 
 #endif
 
 /* Patch version of this package */
 #ifndef FASTJET_VERSION_PATCHLEVEL 
-#define FASTJET_VERSION_PATCHLEVEL  3 
+#define FASTJET_VERSION_PATCHLEVEL  1 
 #endif
 
Index: external/fastjet/config_raw.h
===================================================================
--- external/fastjet/config_raw.h	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/config_raw.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -38,6 +38,9 @@
 #define ENABLE_PLUGIN_TRACKJET /**/
 
-/* defined if C++11 features have been enabled */
-/* #undef HAVE_CXX11_FEATURES */
+/* compile the deprecated parts of the interface using auto-ptr */
+#define HAVE_AUTO_PTR_INTERFACE /**/
+
+/* compiler supports c++14 deprecated keyword */
+/* #undef HAVE_CXX14_DEPRECATED */
 
 /* defined if demangling is enabled at configure time and is supported through
@@ -51,4 +54,10 @@
 #define HAVE_EXECINFO_H 1
 
+/* compiler supports the "explicit" keyword for operators */
+/* #undef HAVE_EXPLICIT_FOR_OPERATORS */
+
+/* compiler supports GNU c++ deprecated attribute */
+#define HAVE_GNUCXX_DEPRECATED /**/
+
 /* Define to 1 if you have the <inttypes.h> header file. */
 #define HAVE_INTTYPES_H 1
@@ -59,4 +68,7 @@
 /* Define to 1 if you have the <memory.h> header file. */
 #define HAVE_MEMORY_H 1
+
+/* compiler supports the "override" keyword */
+/* #undef HAVE_OVERRIDE */
 
 /* Define to 1 if you have the <stdint.h> header file. */
@@ -81,6 +93,5 @@
 #define HAVE_UNISTD_H 1
 
-/* Define to the sub-directory in which libtool stores uninstalled libraries.
-   */
+/* Define to the sub-directory where libtool stores uninstalled libraries. */
 #define LT_OBJDIR ".libs/"
 
@@ -95,11 +106,14 @@
 
 /* Define to the full name and version of this package. */
-#define PACKAGE_STRING "FastJet 3.1.3"
+#define PACKAGE_STRING "FastJet 3.2.1"
 
 /* Define to the one symbol short name of this package. */
 #define PACKAGE_TARNAME "fastjet"
 
+/* Define to the home page for this package. */
+#define PACKAGE_URL ""
+
 /* Define to the version of this package. */
-#define PACKAGE_VERSION "3.1.3"
+#define PACKAGE_VERSION "3.2.1"
 
 /* Define to 1 if you have the ANSI C header files. */
@@ -107,5 +121,5 @@
 
 /* Version number of package */
-#define VERSION "3.1.3"
+#define VERSION "3.2.1"
 
 /* Major version of this package */
@@ -113,11 +127,11 @@
 
 /* Minor version of this package */
-#define VERSION_MINOR 1
+#define VERSION_MINOR 2
 
 /* Version of the package under the form XYYZZ (instead of X.Y.Z) */
-#define VERSION_NUMBER 30103
+#define VERSION_NUMBER 30201
 
 /* Patch version of this package */
-#define VERSION_PATCHLEVEL 3
+#define VERSION_PATCHLEVEL 1
 
 /* Pre-release version of this package */
Index: external/fastjet/config_win.h
===================================================================
--- external/fastjet/config_win.h	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/config_win.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,10 +1,9 @@
-#define FASTJET_PACKAGE_STRING  "FastJet 3.1.3"
-#define FASTJET_PACKAGE_VERSION  "3.1.3"
-#define FASTJET_VERSION  "3.1.3"
-#define FASTJET_VERSION_MAJOR       3.1.3
-#define FASTJET_VERSION_MINOR       3.1.3
-#define FASTJET_VERSION_PATCHLEVEL  3.1.3
-#define FASTJET_VERSION_PRERELEASE  "3.1.3"
-#define FASTJET_VERSION_NUMBER      00000
+#define FASTJET_PACKAGE_STRING  "FastJet 3.2.1"
+#define FASTJET_PACKAGE_VERSION  "3.2.1"
+#define FASTJET_VERSION  "3.2.1"
+#define FASTJET_VERSION_MAJOR       3
+#define FASTJET_VERSION_MINOR       2
+#define FASTJET_VERSION_PATCHLEVEL  1
+#define FASTJET_VERSION_NUMBER      30201
 
 /* The ATLASCone plugin is disabled by default*/
Index: external/fastjet/contribs/Nsubjettiness/ChangeLog
===================================================================
--- external/fastjet/contribs/Nsubjettiness/ChangeLog	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/contribs/Nsubjettiness/ChangeLog	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,2 +1,11 @@
+2016-06-08 <jthaler>
+   Fixed bug in MeasureDefinition.cc where axes were not completely defined,
+      leading to problems with multi-pass axes 
+2016-04-04 <jthaler>
+   Fixed Njettiness.cc to give value of _current_tau_components even if less
+      than N constituents
+   Delete extraneous code in example_advanced_usage.cc
+2016-03-29 <jthaler>
+   Update for FJ 3.2.0 to deal with SharedPtr () deprecation
 2015-09-28 <jthaler>
    Updated NEWS for 2.2.1 release.
Index: external/fastjet/contribs/Nsubjettiness/MeasureDefinition.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/MeasureDefinition.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/contribs/Nsubjettiness/MeasureDefinition.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -5,5 +5,5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
-//  $Id: MeasureDefinition.cc 819 2015-06-12 21:23:24Z jthaler $
+//  $Id: MeasureDefinition.cc 946 2016-06-14 19:11:27Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -524,4 +524,5 @@
       old_axes[k].set_rap( seedAxes[k].rap() );
       old_axes[k].set_phi( seedAxes[k].phi() );
+      old_axes[k].set_mom( seedAxes[k].modp()  );
    }
    
Index: external/fastjet/contribs/Nsubjettiness/NEWS
===================================================================
--- external/fastjet/contribs/Nsubjettiness/NEWS	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/contribs/Nsubjettiness/NEWS	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -32,4 +32,9 @@
     N-jettiness as a jet finder using the new ConicalGeometric measure.
 
+-- 2.2.4:  (Jun 14, 2016) Fixed bug where multi-pass minimization could yield
+           pathological axes (thanks Gregory Soyez) 
+-- 2.2.3:  (Apr 4, 2016) Fixed bug where a jet with fewer than N constituents
+           could give random value for tau_N (thanks Nathan Hartland)
+-- 2.2.2:  (Mar 29, 2016)  Updating SharedPtr interface for FJ 3.2
 -- 2.2.1:  (Sept 28, 2015)  Fix of small Makefile bug
 -- 2.2.0:  (Sept 7, 2015)  Inclusion of the XCone jet algorithm, as well as a
Index: external/fastjet/contribs/Nsubjettiness/Njettiness.cc
===================================================================
--- external/fastjet/contribs/Nsubjettiness/Njettiness.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/contribs/Nsubjettiness/Njettiness.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -5,5 +5,5 @@
 //  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
 //
-//  $Id: Njettiness.cc 821 2015-06-15 18:50:53Z jthaler $
+//  $Id: Njettiness.cc 933 2016-04-04 22:23:32Z jthaler $
 //----------------------------------------------------------------------
 // This file is part of FastJet contrib.
@@ -46,5 +46,5 @@
 // setAxes for Manual mode
 void Njettiness::setAxes(const std::vector<fastjet::PseudoJet> & myAxes) {
-   if (_axes_def()->needsManualAxes()) {
+   if (_axes_def->needsManualAxes()) {
       _currentAxes = myAxes;
    } else {
@@ -59,23 +59,32 @@
       _currentAxes = inputJets;
       _currentAxes.resize(n_jets,fastjet::PseudoJet(0.0,0.0,0.0,0.0));
-      _current_tau_components = TauComponents();
+     
+     // Put in empty tau components
+     std::vector<double> dummy_jet_pieces;
+     _current_tau_components = TauComponents(UNDEFINED_SHAPE,
+                                             dummy_jet_pieces,
+                                             0.0,
+                                             1.0,
+                                             _currentAxes,
+                                             _currentAxes
+                                             );
       _seedAxes = _currentAxes;
       _currentPartition = TauPartition(n_jets); // empty partition
    } else {
-      assert(_axes_def()); // this should never fail.
+      assert(_axes_def); // this should never fail.
       
-      if (_axes_def()->needsManualAxes()) { // if manual mode
+      if (_axes_def->needsManualAxes()) { // if manual mode
          // take current axes as seeds
          _seedAxes = _currentAxes;
          
          // refine axes if requested
-         _currentAxes = _axes_def->get_refined_axes(n_jets,inputJets,_seedAxes, _measure_def());
+         _currentAxes = _axes_def->get_refined_axes(n_jets,inputJets,_seedAxes, _measure_def.get());
       } else { // non-manual axes
          
           //set starting point for minimization
-         _seedAxes = _axes_def->get_starting_axes(n_jets,inputJets,_measure_def());
+         _seedAxes = _axes_def->get_starting_axes(n_jets,inputJets,_measure_def.get());
          
          // refine axes as needed
-         _currentAxes = _axes_def->get_refined_axes(n_jets,inputJets,_seedAxes, _measure_def());
+         _currentAxes = _axes_def->get_refined_axes(n_jets,inputJets,_seedAxes, _measure_def.get());
          
          // NOTE:  The above two function calls are combined in "AxesDefinition::get_axes"
Index: external/fastjet/contribs/Nsubjettiness/VERSION
===================================================================
--- external/fastjet/contribs/Nsubjettiness/VERSION	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/contribs/Nsubjettiness/VERSION	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,1 +1,1 @@
-2.2.1
+2.2.4
Index: external/fastjet/contribs/SoftKiller/Makefile
===================================================================
--- external/fastjet/contribs/SoftKiller/Makefile	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ 	(revision )
@@ -1,81 +1,0 @@
-# If you are using this Makefile standalone and fastjet-config is not
-# in your path, edit this line to specify the full path
-FASTJETCONFIG=fastjet-config
-PREFIX=`$(FASTJETCONFIG) --prefix`
-CXX=g++
-CXXFLAGS= -O3 -Wall -g
-install_script = $(SHELL) ../utils/install-sh
-check_script = ../utils/check.sh
-
-# global contrib-wide Makefile include may override some of the above
-# variables (leading "-" means don't give an error if you can't find
-# the file)
--include ../.Makefile.inc
-
-#------------------------------------------------------------------------
-# things that are specific to this contrib
-NAME=SoftKiller
-SRCS=SoftKiller.cc
-EXAMPLES=example
-INSTALLED_HEADERS=SoftKiller.hh
-#------------------------------------------------------------------------
-
-CXXFLAGS+= $(shell $(FASTJETCONFIG) --cxxflags)
-LDFLAGS += -lm $(shell $(FASTJETCONFIG) --libs)
-
-OBJS  = $(SRCS:.cc=.o)
-EXAMPLES_SRCS  = $(EXAMPLES:=.cc)
-
-install_HEADER  = $(install_script) -c -m 644
-install_LIB     = $(install_script) -c -m 644
-install_DIR     = $(install_script) -d
-install_DATA    = $(install_script) -c -m 644
-install_PROGRAM = $(install_script) -c -s
-install_SCRIPT  = $(install_script) -c
-
-.PHONY: clean distclean examples check install
-
-# compilation of the code (default target)
-all: lib$(NAME).a
-
-lib$(NAME).a: $(OBJS) 
-	ar cru lib$(NAME).a $(OBJS)
-	ranlib lib$(NAME).a
-
-# building the examples
-examples: $(EXAMPLES)
-
-# the following construct makes it possible to automatically build
-# each of the examples listed in $EXAMPLES
-$(EXAMPLES): % : %.o all
-	$(CXX) -o $@ $< -L. -l$(NAME) $(LDFLAGS)
-
-# check that everything went fine
-check: examples
-	@for prog in $(EXAMPLES); do\
-	  $(check_script) $${prog} ../data/Pythia-Zp2jets-lhc-pileup-1ev.dat || exit 1; \
-	done
-	@echo "All tests successful"
-
-# cleaning the directory
-clean:
-	rm -f *~ *.o
-
-distclean: clean
-	rm -f lib$(NAME).a $(EXAMPLES)
-
-# install things in PREFIX/...
-install: all
-	$(install_DIR) $(PREFIX)/include/fastjet/contrib
-	for header in $(INSTALLED_HEADERS); do\
-	  $(install_HEADER) $$header $(PREFIX)/include/fastjet/contrib/;\
-	done
-	$(install_DIR) $(PREFIX)/lib
-	$(install_LIB) lib$(NAME).a $(PREFIX)/lib
-
-depend:
-	makedepend -Y --   -- $(SRCS) $(EXAMPLES_SRCS)
-# DO NOT DELETE
-
-SoftKiller.o: SoftKiller.hh
-example.o: SoftKiller.hh
Index: external/fastjet/internal/ClosestPair2D.hh
===================================================================
--- external/fastjet/internal/ClosestPair2D.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/internal/ClosestPair2D.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClosestPair2D.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClosestPair2D.hh 4059 2016-03-03 20:49:48Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -38,4 +38,5 @@
 #include "fastjet/internal/SearchTree.hh"
 #include "fastjet/internal/MinHeap.hh"
+#include "fastjet/SharedPtr.hh"
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
@@ -134,6 +135,6 @@
 
 
-  triplet<std::auto_ptr<Tree> >  _trees;
-  std::auto_ptr<MinHeap> _heap;
+  triplet<SharedPtr<Tree> >  _trees;
+  SharedPtr<MinHeap>     _heap;
   std::vector<Point>     _points;
   std::stack<Point *>    _available_points;
Index: external/fastjet/internal/ClosestPair2DBase.hh
===================================================================
--- external/fastjet/internal/ClosestPair2DBase.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/internal/ClosestPair2DBase.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: ClosestPair2DBase.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: ClosestPair2DBase.hh 4049 2016-03-03 13:45:47Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -48,5 +48,5 @@
   double x, y;
 
-  Coord2D() {};
+  Coord2D() : x(0.0), y(0.0) {};
 
   Coord2D(double a, double b): x(a), y(b) {};
Index: external/fastjet/internal/DnnPlane.hh
===================================================================
--- external/fastjet/internal/DnnPlane.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/internal/DnnPlane.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: DnnPlane.hh 3918 2015-07-03 14:19:13Z salam $
+// $Id: DnnPlane.hh 3917 2015-07-03 14:07:50Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
Index: external/fastjet/internal/LazyTiling9Alt.hh
===================================================================
--- external/fastjet/internal/LazyTiling9Alt.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/internal/LazyTiling9Alt.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: LazyTiling9Alt.hh 3808 2015-02-20 11:24:53Z soyez $
+// $Id: LazyTiling9Alt.hh 3807 2015-02-20 11:16:55Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
Index: external/fastjet/internal/TilingExtent.hh
===================================================================
--- external/fastjet/internal/TilingExtent.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/internal/TilingExtent.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -5,5 +5,5 @@
 
 //FJSTARTHEADER
-// $Id: TilingExtent.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: TilingExtent.hh 4034 2016-03-02 00:20:27Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -46,4 +46,7 @@
   TilingExtent(ClusterSequence & cs);
 
+  /// constructor that takes a list of PseudoJets
+  TilingExtent(const std::vector<PseudoJet> &particles);
+
   /// returns the suggested minimum rapidity for the tiling
   double minrap() const {return _minrap;}
Index: external/fastjet/internal/base.hh
===================================================================
--- external/fastjet/internal/base.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/internal/base.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,5 +1,5 @@
 
 //FJSTARTHEADER
-// $Id: base.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: base.hh 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -33,4 +33,6 @@
 #define __FASTJET_FASTJET_BASE_HH__
 
+#include "fastjet/config.h"
+
 /// \namespace fastjet
 /// the FastJet namespace
@@ -43,3 +45,11 @@
 #define FASTJET_END_NAMESPACE   }
 
+// define a macro to mark virtual function in derived classes as
+// overriding the base-class definition
+#ifdef FASTJET_HAVE_OVERRIDE
+#define FASTJET_OVERRIDE  override
+#else
+#define FASTJET_OVERRIDE  
+#endif
+
 #endif // __FASTJET_FASTJET_BASE_HH__
Index: external/fastjet/internal/deprecated.hh
===================================================================
--- external/fastjet/internal/deprecated.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
+++ external/fastjet/internal/deprecated.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -0,0 +1,50 @@
+
+//FJSTARTHEADER
+// $Id: deprecated.hh 4098 2016-03-15 16:38:22Z salam $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet.
+//
+//  FastJet 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.
+//
+//  The algorithms that underlie FastJet have required considerable
+//  development. They are described in the original FastJet paper,
+//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
+//  FastJet as part of work towards a scientific publication, please
+//  quote the version you use and include a citation to the manual and
+//  optionally also to hep-ph/0512210.
+//
+//  FastJet 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 FastJet. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+//FJENDHEADER
+
+#ifndef __FASTJET_FASTJET_DEPRECATED_HH__
+#define __FASTJET_FASTJET_DEPRECATED_HH__
+
+#include "fastjet/config.h"
+
+// define a deprecation macro based on the capabilities of the compiler
+// (as determined at configure time).
+#if defined(FASTJET_HAVE_CXX14_DEPRECATED)
+#define FASTJET_DEPRECATED               [[deprecated]]
+#define FASTJET_DEPRECATED_MSG(message)  [[deprecated(message)]]
+#elif defined(FASTJET_HAVE_GNUCXX_DEPRECATED)
+#define FASTJET_DEPRECATED               __attribute__((__deprecated__))
+#define FASTJET_DEPRECATED_MSG(message)  __attribute__((__deprecated__))
+#else
+#define FASTJET_DEPRECATED               
+#define FASTJET_DEPRECATED_MSG(message) 
+#endif
+
+#endif // __FASTJET_FASTJET_DEPRECATED_HH__
Index: external/fastjet/plugins/Jade/JadePlugin.cc
===================================================================
--- external/fastjet/plugins/Jade/JadePlugin.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/Jade/JadePlugin.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: JadePlugin.cc 3433 2014-07-23 08:17:03Z salam $
+// $Id: JadePlugin.cc 4063 2016-03-04 10:31:40Z salam $
 //
 // Copyright (c) 2007-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -35,4 +35,5 @@
 //#include "fastjet/internal/ClusterSequence_N2.icc"
 #include "fastjet/NNH.hh"
+#include "fastjet/NNFJN2Plain.hh"
 
 // other stuff
@@ -51,4 +52,23 @@
 //----------------------------------------------------------------------
 /// class to help run a JADE algorithm
+///
+/// This class works both with NNH and NNFJN2Plain clustering
+/// helpers. They both use the same init(...) call, but for the
+/// clustering:
+///
+/// - NNH uses distance(...) and beam_distance()
+/// - NNFJPlainN2 uses geometrical_distance(...), momentum_factor()
+///   and geometrical_beam_distance()
+///
+/// For NNFJPlainN2 the 2 E_i E_j (1-cos theta_{ij}) factor
+/// gets broken up into
+///
+///     sqrt(2)*min(E_i,E_j) * [sqrt(2)*max(E_i,E_j) (1 - cos \theta_{ij})]
+///
+/// The second factor is what we call the "geometrical_distance" even
+/// though it isn't actually purely geometrical. But the fact that it
+/// gets multiplied by min(E_i,E_j) to get the full distance is
+/// sufficient for the validity of the FJ lemma, allowing for the use
+/// of NNFJN2Plain.
 class JadeBriefJet {
 public:
@@ -69,8 +89,28 @@
   }
 
+  double geometrical_distance(const JadeBriefJet * jet) const {
+    double dij = 1 - nx*jet->nx
+                   - ny*jet->ny
+                   - nz*jet->nz;
+    dij *= max(rt2E,jet->rt2E);
+    return dij;
+  }
+
+  double momentum_factor() const {
+    return rt2E;
+  }
+  
   double beam_distance() const {
     return numeric_limits<double>::max();
   }
 
+  double geometrical_beam_distance() const {
+    // get a number that is almost the same as max(), just a little
+    // smaller so as to ensure that when we divide it by rt2E and then
+    // multiply it again, we won't get an overflow
+    const double almost_max = numeric_limits<double>::max() * (1 - 1e-13);
+    return almost_max / rt2E;
+  }
+  
 private:
   double rt2E, nx, ny, nz;
@@ -82,11 +122,64 @@
   ostringstream desc;
   desc << "e+e- JADE algorithm plugin";
+  switch(_strategy) {
+  case strategy_NNH:
+    desc << ", using NNH strategy"; break;
+  case strategy_NNFJN2Plain:
+    desc << ", using NNFJN2Plain strategy"; break;
+  default:
+    throw Error("Unrecognized strategy in JadePlugin");
+  }
+
   return desc.str();
 }
 
-//----------------------------------------------------------------------
-void JadePlugin::run_clustering(ClusterSequence & cs) const {
+// //----------------------------------------------------------------------
+// void JadePlugin::run_clustering(ClusterSequence & cs) const {
+//   int njets = cs.jets().size();
+// 
+//   //SharedPtr<NNBase<> > nn;
+//   NNBase<> * nn;
+//   switch(_strategy) {
+//   case strategy_NNH:
+//     //nn.reset(new NNH<JadeBriefJet>(cs.jets()));
+//     nn = new NNH<JadeBriefJet>(cs.jets());
+//     break;
+//   case strategy_NNFJN2Plain:
+//     //nn.reset(new NNFJN2Plain<JadeBriefJet>(cs.jets()));
+//     nn = new NNFJN2Plain<JadeBriefJet>(cs.jets());
+//     break;
+//   default:
+//     throw Error("Unrecognized strategy in JadePlugin");
+//   }
+//   //NNH<JadeBriefJet> nnh(cs.jets());
+//   //NNFJN2Plain<JadeBriefJet> nnh(cs.jets());
+// 
+//   // if testing against Hoeth's implementation, need to rescale the
+//   // dij by Q^2.
+//   //double Q2 = cs.Q2(); 
+// 
+//   while (njets > 0) {
+//     int i, j, k;
+//     double dij = nn->dij_min(i, j);
+// 
+//     if (j >= 0) {
+//       cs.plugin_record_ij_recombination(i, j, dij, k);
+//       nn->merge_jets(i, j, cs.jets()[k], k);
+//     } else {
+//       double diB = cs.jets()[i].E()*cs.jets()[i].E(); // get new diB
+//       cs.plugin_record_iB_recombination(i, diB);
+//       nn->remove_jet(i);
+//     }
+//     njets--;
+//   }
+//   delete nn;
+// }
+
+
+template<class N> void JadePlugin::_actual_run_clustering(ClusterSequence & cs) const {
+
   int njets = cs.jets().size();
-  NNH<JadeBriefJet> nnh(cs.jets());
+
+  N nn(cs.jets());
 
   // if testing against Hoeth's implementation, need to rescale the
@@ -96,17 +189,34 @@
   while (njets > 0) {
     int i, j, k;
-    double dij = nnh.dij_min(i, j);
+    double dij = nn.dij_min(i, j);
 
     if (j >= 0) {
       cs.plugin_record_ij_recombination(i, j, dij, k);
-      nnh.merge_jets(i, j, cs.jets()[k], k);
+      nn.merge_jets(i, j, cs.jets()[k], k);
     } else {
       double diB = cs.jets()[i].E()*cs.jets()[i].E(); // get new diB
       cs.plugin_record_iB_recombination(i, diB);
-      nnh.remove_jet(i);
+      nn.remove_jet(i);
     }
     njets--;
   }
+
 }
 
+//----------------------------------------------------------------------
+void JadePlugin::run_clustering(ClusterSequence & cs) const {
+
+  switch(_strategy) {
+  case strategy_NNH:
+    _actual_run_clustering<NNH<JadeBriefJet> >(cs);
+    break;
+  case strategy_NNFJN2Plain:
+    _actual_run_clustering<NNFJN2Plain<JadeBriefJet> >(cs);
+    break;
+  default:
+    throw Error("Unrecognized strategy in JadePlugin");
+  }
+}
+
+
 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
Index: external/fastjet/plugins/Jade/fastjet/JadePlugin.hh
===================================================================
--- external/fastjet/plugins/Jade/fastjet/JadePlugin.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/Jade/fastjet/JadePlugin.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: JadePlugin.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: JadePlugin.hh 4061 2016-03-03 21:51:25Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -77,6 +77,11 @@
 class JadePlugin : public JetDefinition::Plugin {
 public:
+  /// enum that contains the two clustering strategy options; for
+  /// higher multiplicities, strategy_NNFJN2Plain is about a factor of
+  /// two faster.
+  enum Strategy { strategy_NNH = 0, strategy_NNFJN2Plain = 1};
+  
   /// Main constructor for the Jade Plugin class.  
-  JadePlugin (){}
+  JadePlugin (Strategy strategy = strategy_NNFJN2Plain) : _strategy(strategy) {}
 
   /// copy constructor
@@ -100,4 +105,7 @@
 private:
 
+  template<class N> void _actual_run_clustering(ClusterSequence &) const;
+  
+  Strategy _strategy;
 };
 
Index: external/fastjet/plugins/SISCone/SISConePlugin.cc
===================================================================
--- external/fastjet/plugins/SISCone/SISConePlugin.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/SISConePlugin.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -78,7 +78,7 @@
 // static members declaration              //
 /////////////////////////////////////////////
-std::auto_ptr<SISConePlugin>           SISConePlugin::stored_plugin;
-std::auto_ptr<std::vector<PseudoJet> > SISConePlugin::stored_particles;
-std::auto_ptr<Csiscone>                SISConePlugin::stored_siscone;
+SharedPtr<SISConePlugin>           SISConePlugin::stored_plugin;
+SharedPtr<std::vector<PseudoJet> > SISConePlugin::stored_particles;
+SharedPtr<Csiscone>                SISConePlugin::stored_siscone;
 
 
Index: external/fastjet/plugins/SISCone/config.h
===================================================================
--- external/fastjet/plugins/SISCone/config.h	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/config.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,63 +1,117 @@
-/* siscone/config.h.  Generated from config.h.in by configure.  */
+#ifndef _SISCONE_CONFIG_H
+#define _SISCONE_CONFIG_H 1
+ 
+/* siscone/config.h. Generated automatically at end of configure. */
+/* siscone/config_raw.h.  Generated from config.h.in by configure.  */
 /* config.h.in.  Generated from configure.ac by autoheader.  */
 
 /* Define to 1 if you have the <dlfcn.h> header file. */
-#define HAVE_DLFCN_H 1
+#ifndef SISCONE_HAVE_DLFCN_H 
+#define SISCONE_HAVE_DLFCN_H  1 
+#endif
 
 /* Define to 1 if you have the <inttypes.h> header file. */
-#define HAVE_INTTYPES_H 1
+#ifndef SISCONE_HAVE_INTTYPES_H 
+#define SISCONE_HAVE_INTTYPES_H  1 
+#endif
 
 /* Define to 1 if you have the `m' library (-lm). */
-#define HAVE_LIBM 1
+#ifndef SISCONE_HAVE_LIBM 
+#define SISCONE_HAVE_LIBM  1 
+#endif
 
 /* Define to 1 if you have the <memory.h> header file. */
-#define HAVE_MEMORY_H 1
+#ifndef SISCONE_HAVE_MEMORY_H 
+#define SISCONE_HAVE_MEMORY_H  1 
+#endif
 
 /* Define to 1 if you have the <stdint.h> header file. */
-#define HAVE_STDINT_H 1
+#ifndef SISCONE_HAVE_STDINT_H 
+#define SISCONE_HAVE_STDINT_H  1 
+#endif
 
 /* Define to 1 if you have the <stdlib.h> header file. */
-#define HAVE_STDLIB_H 1
+#ifndef SISCONE_HAVE_STDLIB_H 
+#define SISCONE_HAVE_STDLIB_H  1 
+#endif
 
 /* Define to 1 if you have the <strings.h> header file. */
-#define HAVE_STRINGS_H 1
+#ifndef SISCONE_HAVE_STRINGS_H 
+#define SISCONE_HAVE_STRINGS_H  1 
+#endif
 
 /* Define to 1 if you have the <string.h> header file. */
-#define HAVE_STRING_H 1
+#ifndef SISCONE_HAVE_STRING_H 
+#define SISCONE_HAVE_STRING_H  1 
+#endif
 
 /* Define to 1 if you have the <sys/stat.h> header file. */
-#define HAVE_SYS_STAT_H 1
+#ifndef SISCONE_HAVE_SYS_STAT_H 
+#define SISCONE_HAVE_SYS_STAT_H  1 
+#endif
 
 /* Define to 1 if you have the <sys/types.h> header file. */
-#define HAVE_SYS_TYPES_H 1
+#ifndef SISCONE_HAVE_SYS_TYPES_H 
+#define SISCONE_HAVE_SYS_TYPES_H  1 
+#endif
 
 /* Define to 1 if you have the <unistd.h> header file. */
-#define HAVE_UNISTD_H 1
+#ifndef SISCONE_HAVE_UNISTD_H 
+#define SISCONE_HAVE_UNISTD_H  1 
+#endif
 
-/* Define to the sub-directory in which libtool stores uninstalled libraries.
-   */
-#define LT_OBJDIR ".libs/"
+/* Define to the sub-directory where libtool stores uninstalled libraries. */
+#ifndef SISCONE_LT_OBJDIR 
+#define SISCONE_LT_OBJDIR  ".libs/" 
+#endif
 
 /* Name of package */
-#define PACKAGE "siscone"
+#ifndef SISCONE_PACKAGE 
+#define SISCONE_PACKAGE  "siscone" 
+#endif
 
 /* Define to the address where bug reports for this package should be sent. */
-#define PACKAGE_BUGREPORT ""
+#ifndef SISCONE_PACKAGE_BUGREPORT 
+#define SISCONE_PACKAGE_BUGREPORT  "" 
+#endif
 
 /* Define to the full name of this package. */
-#define PACKAGE_NAME "SISCone"
+#ifndef SISCONE_PACKAGE_NAME 
+#define SISCONE_PACKAGE_NAME  "SISCone" 
+#endif
 
 /* Define to the full name and version of this package. */
-#define PACKAGE_STRING "SISCone 3.0.0"
+#ifndef SISCONE_PACKAGE_STRING 
+#define SISCONE_PACKAGE_STRING  "SISCone 3.0.3" 
+#endif
 
 /* Define to the one symbol short name of this package. */
-#define PACKAGE_TARNAME "siscone"
+#ifndef SISCONE_PACKAGE_TARNAME 
+#define SISCONE_PACKAGE_TARNAME  "siscone" 
+#endif
+
+/* Define to the home page for this package. */
+#ifndef SISCONE_PACKAGE_URL 
+#define SISCONE_PACKAGE_URL  "" 
+#endif
 
 /* Define to the version of this package. */
-#define PACKAGE_VERSION "3.0.0"
+#ifndef SISCONE_PACKAGE_VERSION 
+#define SISCONE_PACKAGE_VERSION  "3.0.3" 
+#endif
 
 /* Define to 1 if you have the ANSI C header files. */
-#define STDC_HEADERS 1
+#ifndef SISCONE_STDC_HEADERS 
+#define SISCONE_STDC_HEADERS  1 
+#endif
+
+/* use unique_ptr instead of auto_ptr */
+/* #undef USES_UNIQUE_PTR_AS_AUTO_PTR */
 
 /* Version number of package */
-#define VERSION "3.0.0"
+#ifndef SISCONE_VERSION 
+#define SISCONE_VERSION  "3.0.3" 
+#endif
+ 
+/* once: _SISCONE_CONFIG_H */
+#endif
Index: external/fastjet/plugins/SISCone/config_raw.h
===================================================================
--- external/fastjet/plugins/SISCone/config_raw.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
+++ external/fastjet/plugins/SISCone/config_raw.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -0,0 +1,68 @@
+/* siscone/config_raw.h.  Generated from config.h.in by configure.  */
+/* config.h.in.  Generated from configure.ac by autoheader.  */
+
+/* Define to 1 if you have the <dlfcn.h> header file. */
+#define HAVE_DLFCN_H 1
+
+/* Define to 1 if you have the <inttypes.h> header file. */
+#define HAVE_INTTYPES_H 1
+
+/* Define to 1 if you have the `m' library (-lm). */
+#define HAVE_LIBM 1
+
+/* Define to 1 if you have the <memory.h> header file. */
+#define HAVE_MEMORY_H 1
+
+/* Define to 1 if you have the <stdint.h> header file. */
+#define HAVE_STDINT_H 1
+
+/* Define to 1 if you have the <stdlib.h> header file. */
+#define HAVE_STDLIB_H 1
+
+/* Define to 1 if you have the <strings.h> header file. */
+#define HAVE_STRINGS_H 1
+
+/* Define to 1 if you have the <string.h> header file. */
+#define HAVE_STRING_H 1
+
+/* Define to 1 if you have the <sys/stat.h> header file. */
+#define HAVE_SYS_STAT_H 1
+
+/* Define to 1 if you have the <sys/types.h> header file. */
+#define HAVE_SYS_TYPES_H 1
+
+/* Define to 1 if you have the <unistd.h> header file. */
+#define HAVE_UNISTD_H 1
+
+/* Define to the sub-directory where libtool stores uninstalled libraries. */
+#define LT_OBJDIR ".libs/"
+
+/* Name of package */
+#define PACKAGE "siscone"
+
+/* Define to the address where bug reports for this package should be sent. */
+#define PACKAGE_BUGREPORT ""
+
+/* Define to the full name of this package. */
+#define PACKAGE_NAME "SISCone"
+
+/* Define to the full name and version of this package. */
+#define PACKAGE_STRING "SISCone 3.0.3"
+
+/* Define to the one symbol short name of this package. */
+#define PACKAGE_TARNAME "siscone"
+
+/* Define to the home page for this package. */
+#define PACKAGE_URL ""
+
+/* Define to the version of this package. */
+#define PACKAGE_VERSION "3.0.3"
+
+/* Define to 1 if you have the ANSI C header files. */
+#define STDC_HEADERS 1
+
+/* use unique_ptr instead of auto_ptr */
+/* #undef USES_UNIQUE_PTR_AS_AUTO_PTR */
+
+/* Version number of package */
+#define VERSION "3.0.3"
Index: external/fastjet/plugins/SISCone/defines.h
===================================================================
--- external/fastjet/plugins/SISCone/defines.h	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/defines.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -22,6 +22,6 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 225                                                          $//
-// $Date:: 2008-05-20 16:59:47 +0200 (Tue, 20 May 2008)                     $//
+// $Revision:: 401                                                          $//
+// $Date:: 2016-05-19 16:44:37 +0200 (Thu, 19 May 2016)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -35,7 +35,7 @@
 // defined in siscone.h
 // Otherwise, config.h
-// It is also defined as "PACKAGE_NAME" in config.h but this method 
+// It is also defined as "SISCONE_PACKAGE_NAME" in config.h but this method 
 // might lead to conflicts
-//#define PROGRAM   PACKAGE_NAME 
+//#define PROGRAM   SISCONE_PACKAGE_NAME 
 
 // program version
@@ -43,5 +43,5 @@
 //   siscone::siscone_version
 // defined in siscone.h
-// It is also defined as "VERSION" in config.h but this method 
+// It is also defined as "SISCONE_VERSION" in config.h but this method 
 // might lead to conflicts
 
Index: external/fastjet/plugins/SISCone/fastjet/SISConePlugin.hh
===================================================================
--- external/fastjet/plugins/SISCone/fastjet/SISConePlugin.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/fastjet/SISConePlugin.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -191,7 +191,7 @@
   // part needed for the cache 
   // variables for caching the results and the input
-  static std::auto_ptr<SISConePlugin          > stored_plugin;
-  static std::auto_ptr<std::vector<PseudoJet> > stored_particles;
-  static std::auto_ptr<siscone::Csiscone      > stored_siscone;
+  static SharedPtr<SISConePlugin          > stored_plugin;
+  static SharedPtr<std::vector<PseudoJet> > stored_particles;
+  static SharedPtr<siscone::Csiscone      > stored_siscone;
 };
 
Index: external/fastjet/plugins/SISCone/fastjet/SISConeSphericalPlugin.hh
===================================================================
--- external/fastjet/plugins/SISCone/fastjet/SISConeSphericalPlugin.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/fastjet/SISConeSphericalPlugin.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -166,7 +166,7 @@
   // part needed for the cache 
   // variables for caching the results and the input
-  static std::auto_ptr<SISConeSphericalPlugin        > stored_plugin;
-  static std::auto_ptr<std::vector<PseudoJet>        > stored_particles;
-  static std::auto_ptr<siscone_spherical::CSphsiscone> stored_siscone;
+  static SharedPtr<SISConeSphericalPlugin        > stored_plugin;
+  static SharedPtr<std::vector<PseudoJet>        > stored_particles;
+  static SharedPtr<siscone_spherical::CSphsiscone> stored_siscone;
 };
 
Index: external/fastjet/plugins/SISCone/siscone.cc
===================================================================
--- external/fastjet/plugins/SISCone/siscone.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/siscone.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -21,16 +21,9 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 371                                                          $//
-// $Date:: 2014-09-09 10:05:32 +0200 (Tue, 09 Sep 2014)                     $//
+// $Revision:: 403                                                          $//
+// $Date:: 2016-05-19 16:52:05 +0200 (Thu, 19 May 2016)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
-//#ifdef HAVE_CONFIG_H
 #include "config.h"
-//#else
-//#define PACKAGE_NAME "SISCone"
-//#define VERSION "3.0.0"
-//#warning "No config.h file available, using preset values"
-//#endif
-
 #include "ranlux.h"
 #include "momentum.h"
@@ -254,4 +247,6 @@
   // print the banner
   if (_banner_ostr != 0){
+    ios::fmtflags flags_to_restore(_banner_ostr->flags());
+
     (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl;
     (*_banner_ostr) << "#                    SISCone   version " << setw(28) << left << siscone_version() << "o" << endl;
@@ -269,4 +264,5 @@
 
     _banner_ostr->flush();
+    _banner_ostr->flags(flags_to_restore);
   }
 }
@@ -279,10 +275,10 @@
  * return SISCone package name.
  * This is nothing but "SISCone", it is a replacement to the
- * PACKAGE_NAME string defined in config.h and which is not
- * public by default.
+ * SISCONE_PACKAGE_NAME string defined in config.h and which is not
+ * guaranteed to be public.
  * return the SISCone name as a string
  */
 string siscone_package_name(){
-  return PACKAGE_NAME;
+  return SISCONE_PACKAGE_NAME;
 }
 
@@ -293,6 +289,6 @@
  */
 string siscone_version(){
-  return VERSION;
-}
-
-}
+  return SISCONE_VERSION;
+}
+
+}
Index: external/fastjet/plugins/SISCone/siscone.h
===================================================================
--- external/fastjet/plugins/SISCone/siscone.h	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/siscone.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -22,6 +22,6 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 369                                                          $//
-// $Date:: 2014-09-04 16:57:55 +0200 (Thu, 04 Sep 2014)                     $//
+// $Revision:: 401                                                          $//
+// $Date:: 2016-05-19 16:44:37 +0200 (Thu, 19 May 2016)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -160,6 +160,6 @@
  * return SISCone package name.
  * This is nothing but "SISCone", it is a replacement to the
- * PACKAGE_NAME string defined in config.h and which is not
- * public by default.
+ * SISCONE_PACKAGE_NAME string defined in config.h and which is not
+ * guaranteed to be public.
  * \return the SISCone name as a string
  */
Index: external/fastjet/plugins/SISCone/split_merge.cc
===================================================================
--- external/fastjet/plugins/SISCone/split_merge.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/split_merge.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -21,6 +21,6 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 370                                                          $//
-// $Date:: 2014-09-04 17:03:15 +0200 (Thu, 04 Sep 2014)                     $//
+// $Revision:: 390                                                          $//
+// $Date:: 2016-03-03 11:06:52 +0100 (Thu, 03 Mar 2016)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -54,4 +54,6 @@
   pt_tilde = 0.0;
   sm_var2 = 0.0;
+  pass = CJET_INEXISTENT_PASS; // initialised to a value that should
+                               // notappear in the end (after clustering)
 }
 
@@ -659,5 +661,5 @@
 
 #ifdef DEBUG_SPLIT_MERGE
-  cout << "PR-Jet " << jets.size() << " [size " << next_jet.contents.size() << "]:";
+  cout << "PR-Jet " << jets.size() << " [size " << jet.contents.size() << "]:";
 #endif
     
@@ -1179,6 +1181,6 @@
   }
 
-  return 0.0;
-}
-
-}
+  //return 0.0;
+}
+
+}
Index: external/fastjet/plugins/SISCone/split_merge.h
===================================================================
--- external/fastjet/plugins/SISCone/split_merge.h	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/split_merge.h	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -22,6 +22,6 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 367                                                          $//
-// $Date:: 2014-09-04 15:57:37 +0200 (Thu, 04 Sep 2014)                     $//
+// $Revision:: 405                                                          $//
+// $Date:: 2016-05-23 20:15:02 +0200 (Mon, 23 May 2016)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -29,4 +29,5 @@
 #define __SPLIT_MERGE_H__
 
+#include "config.h"
 #include "defines.h"
 #include "geom_2d.h"
@@ -39,4 +40,6 @@
 
 namespace siscone{
+
+const int CJET_INEXISTENT_PASS = -2;
 
 /**
@@ -75,4 +78,6 @@
   /// pass at which the jet has been found
   /// It starts at 0 (first pass), -1 means infinite rapidity
+  /// (it will be initialised to "CJET_INEXISTENT_PASS" which should
+  /// never appear after clustering)
   int pass;
 };
@@ -441,6 +446,10 @@
   // jet information
   /// list of jet candidates
+#ifdef SISCONE_USES_UNIQUE_PTR_AS_AUTO_PTR
+  std::unique_ptr<std::multiset<Cjet,Csplit_merge_ptcomparison> > candidates;
+#else
   std::auto_ptr<std::multiset<Cjet,Csplit_merge_ptcomparison> > candidates;
-
+#endif
+  
   /// minimal pt2
   double pt_min2;
Index: external/fastjet/plugins/SISCone/vicinity.cc
===================================================================
--- external/fastjet/plugins/SISCone/vicinity.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/plugins/SISCone/vicinity.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -21,6 +21,6 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 123                                                          $//
-// $Date:: 2007-03-01 02:52:16 +0100 (Thu, 01 Mar 2007)                     $//
+// $Revision:: 388                                                          $//
+// $Date:: 2016-03-03 10:42:25 +0100 (Thu, 03 Mar 2016)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -80,4 +80,5 @@
   VR2 = VR = 0.0;
 
+  ve_list = NULL;
   set_particle_list(_particle_list);
 }
Index: external/fastjet/tools/Filter.cc
===================================================================
--- external/fastjet/tools/Filter.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/tools/Filter.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: Filter.cc 3760 2014-12-19 10:05:10Z soyez $
+// $Id: Filter.cc 4080 2016-03-09 15:01:57Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -31,4 +31,5 @@
 #include "fastjet/tools/Filter.hh"
 #include "fastjet/tools/Recluster.hh"
+#include "fastjet/tools/Subtractor.hh"
 #include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
 #include <cassert>
@@ -98,7 +99,8 @@
   } else if (_rho!=0){
     if (subjets.size()>0){
-      const ClusterSequenceAreaBase *csab = subjets[0].validated_csab();
+      //const ClusterSequenceAreaBase *csab = subjets[0].validated_csab();
       for (unsigned int i=0;i<subjets.size();i++){
-        subjets[i]=csab->subtracted_jet(subjets[i], _rho);
+        //subjets[i]=csab->subtracted_jet(subjets[i], _rho);
+        subjets[i]=Subtractor(_rho)(subjets[i]);
       }
     }
Index: external/fastjet/tools/Filter.hh
===================================================================
--- external/fastjet/tools/Filter.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/tools/Filter.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: Filter.hh 3845 2015-03-08 08:35:36Z soyez $
+// $Id: Filter.hh 3846 2015-03-08 08:35:43Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -143,4 +143,8 @@
   /// deciding which ones to keep. It takes precedence over a non-zero rho.
   void set_subtractor(const FunctionOfPseudoJet<PseudoJet> * subtractor_in) {_subtractor = subtractor_in;}
+
+  /// Set a subtractor that is applied to all individual subjets before
+  /// deciding which ones to keep. It takes precedence over a non-zero rho.
+  const FunctionOfPseudoJet<PseudoJet> * subtractor() const{ return _subtractor;}
 
   /// runs the filtering and sets kept and rejected to be the jets of interest
Index: external/fastjet/tools/GridMedianBackgroundEstimator.hh
===================================================================
--- external/fastjet/tools/GridMedianBackgroundEstimator.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/tools/GridMedianBackgroundEstimator.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -3,5 +3,5 @@
 
 //FJSTARTHEADER
-// $Id: GridMedianBackgroundEstimator.hh 3778 2014-12-24 09:28:09Z salam $
+// $Id: GridMedianBackgroundEstimator.hh 3969 2015-09-21 08:57:59Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -194,5 +194,5 @@
   /// determination of rho_m.
   ///
-  /// Note that support for sigma_m is automatic is one has sigma and
+  /// Note that support for sigma_m is automatic if one has sigma and
   /// rho_m support.
   bool has_rho_m() const {return _enable_rho_m;}
Index: external/fastjet/tools/JetMedianBackgroundEstimator.cc
===================================================================
--- external/fastjet/tools/JetMedianBackgroundEstimator.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/tools/JetMedianBackgroundEstimator.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: JetMedianBackgroundEstimator.cc 3517 2014-08-01 14:23:13Z soyez $
+// $Id: JetMedianBackgroundEstimator.cc 4047 2016-03-03 13:21:49Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -194,5 +194,5 @@
 
   _csi = jets[0].structure_shared_ptr();
-  ClusterSequenceStructure * csi = dynamic_cast<ClusterSequenceStructure*>(_csi());
+  ClusterSequenceStructure * csi = dynamic_cast<ClusterSequenceStructure*>(_csi.get());
   const ClusterSequenceAreaBase * csab = csi->validated_csab();
 
@@ -467,5 +467,5 @@
 
   // determine the number of empty jets
-  const ClusterSequenceAreaBase * csab = (dynamic_cast<ClusterSequenceStructure*>(_csi()))->validated_csab();
+  const ClusterSequenceAreaBase * csab = (dynamic_cast<ClusterSequenceStructure*>(_csi.get()))->validated_csab();
   if (csab->has_explicit_ghosts()) {
     _empty_area = 0.0;
@@ -503,9 +503,9 @@
 // throw an error otherwise
 void JetMedianBackgroundEstimator::_check_csa_alive() const{
-  ClusterSequenceStructure* csa = dynamic_cast<ClusterSequenceStructure*>(_csi());
+  ClusterSequenceStructure* csa = dynamic_cast<ClusterSequenceStructure*>(_csi.get());
   if (csa == 0) {
     throw Error("JetMedianBackgroundEstimator: there is no cluster sequence associated with the JetMedianBackgroundEstimator");
   }
-  if (! dynamic_cast<ClusterSequenceStructure*>(_csi())->has_associated_cluster_sequence())
+  if (! dynamic_cast<ClusterSequenceStructure*>(_csi.get())->has_associated_cluster_sequence())
     throw Error("JetMedianBackgroundEstimator: modifications are no longer possible as the underlying ClusterSequence has gone out of scope");
 }
@@ -521,5 +521,5 @@
   // cluster sequence
   if (_jet_def.jet_algorithm() == undefined_jet_algorithm){
-    const ClusterSequence * cs = dynamic_cast<ClusterSequenceStructure*>(_csi())->validated_cs();
+    const ClusterSequence * cs = dynamic_cast<ClusterSequenceStructure*>(_csi.get())->validated_cs();
     jet_def = &(cs->jet_def());
   }
Index: external/fastjet/tools/Recluster.hh
===================================================================
--- external/fastjet/tools/Recluster.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/tools/Recluster.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -2,5 +2,5 @@
 #define __FASTJET_TOOLS_RECLUSTER_HH__
 
-// $Id: Recluster.hh 3760 2014-12-19 10:05:10Z soyez $
+// $Id: Recluster.hh 3753 2014-12-17 15:19:55Z salam $
 //
 // Copyright (c) 2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
Index: external/fastjet/tools/Subtractor.cc
===================================================================
--- external/fastjet/tools/Subtractor.cc	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/tools/Subtractor.cc	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: Subtractor.cc 3670 2014-09-08 14:17:59Z soyez $
+// $Id: Subtractor.cc 3970 2015-09-21 10:31:17Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -39,4 +39,5 @@
 const double Subtractor::_invalid_rho = -numeric_limits<double>::infinity();
 
+LimitedWarning Subtractor::_unused_rho_m_warning;
 
 //----------------------------------------------------------------------
Index: external/fastjet/tools/Subtractor.hh
===================================================================
--- external/fastjet/tools/Subtractor.hh	(revision 973b92a3bf3d2d9073e3b7917710e4fac19714c7)
+++ external/fastjet/tools/Subtractor.hh	(revision 4154bbda8cb3b2dc05c3f4258ee0b6bf25132273)
@@ -1,4 +1,4 @@
 //FJSTARTHEADER
-// $Id: Subtractor.hh 3670 2014-09-08 14:17:59Z soyez $
+// $Id: Subtractor.hh 3970 2015-09-21 10:31:17Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -203,5 +203,5 @@
   static const double _invalid_rho;
 
-  mutable LimitedWarning _unused_rho_m_warning;
+  static LimitedWarning _unused_rho_m_warning;
 };
 
