Index: external/fastjet/tools/BackgroundEstimatorBase.cc
===================================================================
--- external/fastjet/tools/BackgroundEstimatorBase.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/BackgroundEstimatorBase.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: BackgroundEstimatorBase.cc 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 
Index: external/fastjet/tools/BackgroundEstimatorBase.hh
===================================================================
--- external/fastjet/tools/BackgroundEstimatorBase.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/BackgroundEstimatorBase.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -2,8 +2,8 @@
 #define __FASTJET_BACKGROUND_ESTIMATOR_BASE_HH__
 
-//STARTHEADER
-// $Id: BackgroundEstimatorBase.hh 2689 2011-11-14 14:51:06Z soyez $
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: BackgroundEstimatorBase.hh 3516 2014-08-01 14:07:58Z salam $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -16,7 +16,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -28,5 +30,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/ClusterSequenceAreaBase.hh>
@@ -43,7 +45,6 @@
 ///
 /// Abstract base class that provides the basic interface for classes
-/// that estimate levels of background radiation in hadrion and
+/// that estimate levels of background radiation in hadron and
 /// heavy-ion collider events.
-///
 ///
 class BackgroundEstimatorBase {
@@ -99,6 +100,40 @@
   /// determination of sigma
   virtual bool has_sigma() {return false;}
-  //\}
-  
+
+  //----------------------------------------------------------------
+  // now do the same thing for rho_m and sigma_m
+
+  /// returns rho_m, the purely longitudinal, particle-mass-induced
+  /// component of the background density per unit area
+  virtual double rho_m() const{
+    throw Error("rho_m() not supported for this Background Estimator");
+  }
+
+  /// returns sigma_m, a measure of the fluctuations in the purely
+  /// longitudinal, particle-mass-induced component of the background
+  /// density per unit area; must be multipled by sqrt(area) to get
+  /// fluctuations for a region of a given area.
+  virtual double sigma_m() const { 
+    throw Error("sigma_m() not supported for this Background Estimator");
+  }
+
+  /// Returns rho_m locally at the jet position. As for rho(jet), it is non-const.
+  virtual double rho_m(const PseudoJet & /*jet*/){
+    throw Error("rho_m(jet) not supported for this Background Estimator");
+  }
+
+  /// Returns sigma_m locally at the jet position. As for rho(jet), it is non-const.
+  virtual double sigma_m(const PseudoJet & /*jet*/) { 
+    throw Error("sigma_m(jet) not supported for this Background Estimator");
+  }
+
+  /// Returns true if this background estimator has support for
+  /// determination of rho_m.
+  ///
+  /// Note that support for sigma_m is automatic is one has sigma and
+  /// rho_m support.
+  virtual bool has_rho_m() const {return false;}
+  //\}
+
 
   /// @name configuring the behaviour
@@ -114,4 +149,7 @@
   /// The BackgroundRescalingYPolynomial class can be used to get a
   /// rescaling that depends just on rapidity.
+  ///
+  /// There is currently no support for different rescaling classes 
+  /// for rho and rho_m determinations.
   virtual void set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class_in) { _rescaling_class = rescaling_class_in; }
 
Index: external/fastjet/tools/Boost.hh
===================================================================
--- external/fastjet/tools/Boost.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/Boost.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id: Boost.hh 2689 2011-11-14 14:51:06Z soyez $
+//FJSTARTHEADER
+// $Id: Boost.hh 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #ifndef __FASTJET_TOOL_BOOST_HH__
Index: external/fastjet/tools/CASubJetTagger.cc
===================================================================
--- external/fastjet/tools/CASubJetTagger.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/CASubJetTagger.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: CASubJetTagger.cc 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/tools/CASubJetTagger.hh>
Index: external/fastjet/tools/CASubJetTagger.hh
===================================================================
--- external/fastjet/tools/CASubJetTagger.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/CASubJetTagger.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id: CASubJetTagger.hh 2616 2011-09-30 18:03:40Z salam $
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: CASubJetTagger.hh 3433 2014-07-23 08:17:03Z salam $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #ifndef __CASUBJET_TAGGER_HH__
Index: external/fastjet/tools/Filter.cc
===================================================================
--- external/fastjet/tools/Filter.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/Filter.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: Filter.cc 3633 2014-08-15 13:23:52Z soyez $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,7 +27,8 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include "fastjet/tools/Filter.hh"
+#include "fastjet/tools/Recluster.hh"
 #include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
 #include <cassert>
@@ -45,4 +48,8 @@
 // class description
 string Filter::description() const {
+  if (!_initialised){
+    return "uninitialised Filter";
+  }
+
   ostringstream ostr;
   ostr << "Filter with subjet_def = ";
@@ -73,4 +80,9 @@
 // by the filtering
 PseudoJet Filter::result(const PseudoJet &jet) const {
+  if (!_initialised){
+    //Q: do we throw or do we return an empty PJ?
+    throw Error("uninitialised Filter");
+  }
+
   // start by getting the list of subjets (including a list of sanity
   // checks)
@@ -78,10 +90,18 @@
   //     _set_filtered_elements_cafilt)
   vector<PseudoJet> subjets; 
-  JetDefinition subjet_def;
-  bool discard_area;
-  // NB: on return, subjet_def is set to the jet definition actually
-  //     used (so that we can make use of its recombination scheme 
-  //     when joining the jets to be kept).
-  _set_filtered_elements(jet, subjets, subjet_def, discard_area);
+  //JetDefinition subjet_def;
+  bool ca_optimised = _set_filtered_elements(jet, subjets);
+
+  // apply subtraction if needed:
+  if (_subtractor){
+    subjets = (*_subtractor)(subjets);
+  } else if (_rho!=0){
+    if (subjets.size()>0){
+      const ClusterSequenceAreaBase *csab = subjets[0].validated_csab();
+      for (unsigned int i=0;i<subjets.size();i++){
+        subjets[i]=csab->subtracted_jet(subjets[i], _rho);
+      }
+    }
+  }
 
   // now build the vector of kept and rejected subjets
@@ -94,199 +114,25 @@
 
   // gather the info under the form of a PseudoJet
-  return _finalise(jet, kept, rejected, subjet_def, discard_area);
+  return _finalise(jet, kept, rejected, ca_optimised);
 }
 
 
 // sets filtered_elements to be all the subjets on which filtering will work
-void Filter::_set_filtered_elements(const PseudoJet & jet,
-                                    vector<PseudoJet> & filtered_elements,
-                                    JetDefinition & subjet_def,
-                                    bool & discard_area) const {
-  // sanity checks
-  //-------------------------------------------------------------------
-  // make sure that the jet has constituents
-  if (! jet.has_constituents())
-    throw Error("Filter can only be applied on jets having constituents");
+//
+// return true when the subjets have been optained using teh optimised
+// method for C/A
+bool Filter::_set_filtered_elements(const PseudoJet & jet,
+                                    vector<PseudoJet> & filtered_elements) const {
+  // create the recluster instance
+  Recluster recluster;
+  if ((_Rfilt>=0) || (_Rfiltfunc))
+    recluster = Recluster(cambridge_algorithm, (_Rfiltfunc) ? (*_Rfiltfunc)(jet) : _Rfilt, Recluster::keep_all);
+  else
+    recluster = Recluster(_subjet_def, false, Recluster::keep_all);
 
-  // for a whole variety of checks, we shall need the "recursive"
-  // pieces of the jet (the jet itself or recursing down to its most
-  // fundamental pieces).
-  // So we do compute these once and for all
-  vector<PseudoJet> all_pieces; //.clear();
-  if ((!_get_all_pieces(jet, all_pieces)) || (all_pieces.size()==0))
-    throw Error("Attempt to filter a jet that has no associated ClusterSequence or is not a superposition of jets associated with a ClusterSequence");
-  
-  // if the filter uses subtraction, make sure we have a CS that supports area and has
-  // explicit ghosts 
-  if (_uses_subtraction()) {
-    if (!jet.has_area())   
-      throw Error("Attempt to filter and subtract (non-zero rho or subtractor) without area info for the original jet");
-
-    if (!_check_explicit_ghosts(all_pieces))
-      throw Error("Attempt to filter and subtract (non-zero rho or subtractor) without explicit ghosts");
-  }
-
-  // if we're dealing with a dynamic determination of the filtering
-  // radius, do it now
-  if ((_Rfilt>=0) || (_Rfiltfunc)){
-    double Rfilt = (_Rfiltfunc) ? (*_Rfiltfunc)(jet) : _Rfilt;
-    const JetDefinition::Recombiner * common_recombiner = _get_common_recombiner(all_pieces);
-    if (common_recombiner) {
-      if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
-        RecombinationScheme scheme = 
-          static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
-        subjet_def = JetDefinition(cambridge_algorithm, Rfilt, scheme);
-      } else {
-        subjet_def = JetDefinition(cambridge_algorithm, Rfilt, common_recombiner);
-      }
-    } else {
-      subjet_def = JetDefinition(cambridge_algorithm, Rfilt);
-    }
-  } else {
-    subjet_def = _subjet_def;
-  }
-
-  // get the jet definition to be use and whether we can apply our
-  // simplified C/A+C/A filter
-  //
-  // we apply C/A clustering iff
-  //  - the request subjet_def is C/A
-  //  - the jet is either directly coming from C/A or if it is a
-  //    superposition of C/A jets
-  //  - the pieces agree with the recombination scheme of subjet_def
-  //------------------------------------------------------------------
-  bool simple_cafilt = _check_ca(all_pieces);
-
-  // extract the subjets
-  //-------------------------------------------------------------------
-  discard_area = false;
-  if (simple_cafilt){
-    // first make sure that 'filtered_elements' is empty
-    filtered_elements.clear();
-    _set_filtered_elements_cafilt(jet, filtered_elements, subjet_def.R());
-    // in the following case, areas can be erroneous and will be discarded
-    discard_area = (!_uses_subtraction()) && (jet.has_area()) && (!_check_explicit_ghosts(all_pieces));
-  } else {
-    // here we'll simply recluster the jets.
-    //
-    // To include an area support we need
-    //  - the jet to have an area
-    //  - subtraction requested or explicit ghosts
-    bool do_areas = (jet.has_area()) && ((_uses_subtraction()) || (_check_explicit_ghosts(all_pieces))); 
-    _set_filtered_elements_generic(jet, filtered_elements, subjet_def, do_areas);
-  }
-
-  // order the filtered elements in pt
-  filtered_elements = sorted_by_pt(filtered_elements);
+  // get the subjets
+  //JetDefinition subjet_def;
+  return recluster.get_new_jets_and_def(jet, filtered_elements);
 }
-
-// set the filtered elements in the simple case of C/A+C/A
-//
-// WATCH OUT: this could be recursively called, so filtered elements
-//            of 'jet' are APPENDED to 'filtered_elements'
-void Filter::_set_filtered_elements_cafilt(const PseudoJet & jet, 
-                                           vector<PseudoJet> & filtered_elements, 
-                                           double Rfilt) const{
-  // we know that the jet is either a C/A jet or a superposition of
-  // such pieces
-  if (jet.has_associated_cluster_sequence()){
-    // just extract the exclusive subjets of 'jet'
-    const ClusterSequence *cs = jet.associated_cluster_sequence(); 
-    vector<PseudoJet> local_fe;
-
-    double dcut = Rfilt / cs->jet_def().R();
-    if (dcut>=1.0){
-      local_fe.push_back(jet);
-    } else {
-      dcut *= dcut;
-      local_fe = jet.exclusive_subjets(dcut);
-    }
-
-    // subtract the jets if needed
-    // Note that this one would work on pieces!!
-    //-----------------------------------------------------------------
-    if (_uses_subtraction()){
-      const ClusterSequenceAreaBase * csab = jet.validated_csab();
-      for (unsigned int i=0;i<local_fe.size();i++) {
-        if (_subtractor) {
-          local_fe[i] = (*_subtractor)(local_fe[i]);
-        } else {
-          local_fe[i] = csab->subtracted_jet(local_fe[i], _rho);
-        }
-      }
-    }
-
-    copy(local_fe.begin(), local_fe.end(), back_inserter(filtered_elements));
-    return;
-  }
-
-  // just recurse into the pieces
-  const vector<PseudoJet> & pieces = jet.pieces();
-  for (vector<PseudoJet>::const_iterator it = pieces.begin(); 
-       it!=pieces.end(); it++)
-    _set_filtered_elements_cafilt(*it, filtered_elements, Rfilt);
-}
-
-
-// set the filtered elements in the generic re-clustering case (wo
-// subtraction)
-void Filter::_set_filtered_elements_generic(const PseudoJet & jet, 
-                                            vector<PseudoJet> & filtered_elements,
-                                            const JetDefinition & subjet_def,
-					    bool do_areas) const{
-  // create a new, internal, ClusterSequence from the jet constituents
-  // get the subjets directly from there
-  //
-  // If the jet has area support then we separate the ghosts from the
-  // "regular" particles so the subjets will also have area
-  // support. Note that we do this regardless of whether rho is zero
-  // or not.
-  //
-  // Note that to be able to separate the ghosts, one needs explicit
-  // ghosts!!
-  // ---------------------------------------------------------------
-  if (do_areas){
-    vector<PseudoJet> all_constituents = jet.constituents();
-    vector<PseudoJet> regular_constituents, ghosts;  
-
-    for (vector<PseudoJet>::iterator it = all_constituents.begin(); 
-         it != all_constituents.end(); it++){
-      if (it->is_pure_ghost())
-        ghosts.push_back(*it);
-      else
-        regular_constituents.push_back(*it);
-    }
-
-    // figure out the ghost area from the 1st ghost (if none, any value
-    // would probably do as the area will be 0 and subtraction will have
-    // no effect!)
-    double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
-    ClusterSequenceActiveAreaExplicitGhosts * csa
-      = new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents, 
-                                                    subjet_def, 
-                                                    ghosts, ghost_area);
-
-    // get the subjets: we use the subtracted or unsubtracted ones
-    // depending on rho or _subtractor being non-zero
-    if (_uses_subtraction()) {
-      if (_subtractor) {
-        filtered_elements = (*_subtractor)(csa->inclusive_jets());
-      } else {
-        filtered_elements = csa->subtracted_jets(_rho);
-      }
-    } else {
-      filtered_elements = csa->inclusive_jets();
-    }
-
-    // allow the cs to be deleted when it's no longer used
-    csa->delete_self_when_unused();
-  } else {
-    ClusterSequence * cs = new ClusterSequence(jet.constituents(), subjet_def);
-    filtered_elements = cs->inclusive_jets();
-    // allow the cs to be deleted when it's no longer used
-    cs->delete_self_when_unused();
-  }
-}
-
 
 // gather the information about what is kept and rejected under the
@@ -295,118 +141,36 @@
                             vector<PseudoJet> & kept, 
                             vector<PseudoJet> & rejected,
-                            const JetDefinition & subjet_def,
-                            const bool discard_area) const {
-  // figure out which recombiner to use
-  const JetDefinition::Recombiner &rec = *(subjet_def.recombiner());
+			    bool ca_optimisation_used) const {
+  PseudoJet filtered_jet;
 
-  // create an appropriate structure and transfer the info to it
-  PseudoJet filtered_jet = join<StructureType>(kept, rec);
+  if (kept.size()+rejected.size()>0){
+    // figure out which recombiner to use
+    const JetDefinition::Recombiner &rec = (kept.size()>0)
+      ? *(kept[0].associated_cs()->jet_def().recombiner())
+      : *(rejected[0].associated_cs()->jet_def().recombiner());
+
+    // create an appropriate structure and transfer the info to it
+    filtered_jet = join<StructureType>(kept, rec);
+  } else {
+    filtered_jet = join<StructureType>(kept);
+  }
   StructureType *fs = (StructureType*) filtered_jet.structure_non_const_ptr();
-//  fs->_original_jet = jet;
   fs->_rejected = rejected;
+  
+  // if we've used C/A optimisation, we need to get rid of the area
+  // information if it comes from a non-explicit-ghost clustering.
+  // (because in that case it can be erroneous due the lack of
+  // information about empty areas)
+  if ((ca_optimisation_used) && (kept.size()+rejected.size()>0)){
+    bool has_non_explicit_ghost_area = (kept.size()>0)
+      ? (kept[0].has_area()     && kept[0].validated_csab()->has_explicit_ghosts())
+      : (rejected[0].has_area() && rejected[0].validated_csab()->has_explicit_ghosts());
+    if (has_non_explicit_ghost_area)
+      fs->discard_area();
+  }
 
-  if (discard_area){
-    // safety check: make sure there is an area to discard!!!
-    assert(fs->_area_4vector_ptr);
-    delete fs->_area_4vector_ptr;
-    fs->_area_4vector_ptr=0;
-  }
-  
   return filtered_jet;
 }
 
-// various checks
-//----------------------------------------------------------------------
-
-// get the pieces down to the fundamental pieces
-// 
-// Note that this just checks that there is an associated CS to the
-// fundamental pieces, not that it is still valid
-bool Filter::_get_all_pieces(const PseudoJet &jet, vector<PseudoJet> &all_pieces) const{
-  if (jet.has_associated_cluster_sequence()){
-    all_pieces.push_back(jet);
-    return true;
-  }
-
-  if (jet.has_pieces()){
-    const vector<PseudoJet> pieces = jet.pieces();
-    for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
-      if (!_get_all_pieces(*it, all_pieces)) return false;
-    return true;
-  }
-
-  return false;
-}
-
-// get the common recombiner to all pieces (NULL if none)
-//
-// Note that if the jet has an associated cluster sequence that is no
-// longer valid, an error will be thrown (needed since it could be the
-// 1st check called after the enumeration of the pieces)
-const JetDefinition::Recombiner* Filter::_get_common_recombiner(const vector<PseudoJet> &all_pieces) const{
-  const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
-  for (unsigned int i=1; i<all_pieces.size(); i++)
-    if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref)) return NULL;
-
-  return jd_ref.recombiner();
-}
-
-// check if the jet (or all its pieces) have explicit ghosts
-// (assuming the jet has area support
-//
-// Note that if the jet has an associated cluster sequence that is no
-// longer valid, an error will be thrown (needed since it could be the
-// 1st check called after the enumeration of the pieces)
-bool Filter::_check_explicit_ghosts(const vector<PseudoJet> &all_pieces) const{
-  for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
-    if (! it->validated_csab()->has_explicit_ghosts()) return false;
-  return true;
-}
-
-// check if one can apply the simplification for C/A subjets
-//
-// This includes:
-//  - the subjet definition asks for C/A subjets
-//  - all the pieces share the same CS
-//  - that CS is C/A with the same recombiner as the subjet def
-//  - the filtering radius is not larger than any of the pairwise
-//    distance between the pieces
-//
-// Note that if the jet has an associated cluster sequence that is no
-// longer valid, an error will be thrown (needed since it could be the
-// 1st check called after the enumeration of the pieces)
-bool Filter::_check_ca(const vector<PseudoJet> &all_pieces) const{
-  if (_subjet_def.jet_algorithm() != cambridge_algorithm) return false;
-
-  // check that the 1st of all the pieces (we're sure there is at
-  // least one) is coming from a C/A clustering. Then check that all
-  // the following pieces share the same ClusterSequence
-  const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
-  if (cs_ref->jet_def().jet_algorithm() != cambridge_algorithm) return false;
-  for (unsigned int i=1; i<all_pieces.size(); i++)
-    if (all_pieces[i].validated_cs() != cs_ref) return false;
-
-  // check that the 1st peice has the same recombiner as the one used
-  // for the subjet clustering
-  // Note that since they share the same CS, checking the 2st one is enough
-  if (!cs_ref->jet_def().has_same_recombiner(_subjet_def)) return false;
-
-  // we also have to make sure that the filtering radius is not larger
-  // than any of the inter-pieces distance
-  double Rfilt2 = _subjet_def.R();
-  Rfilt2 *= Rfilt2;
-  for (unsigned int i=0; i<all_pieces.size()-1; i++){
-    for (unsigned int j=i+1; j<all_pieces.size(); j++){
-      if (all_pieces[i].squared_distance(all_pieces[j]) <  Rfilt2) return false;
-    }
-  }
-
-  return true;
-}
-
-//----------------------------------------------------------------------
-// FilterInterface implementation 
-//----------------------------------------------------------------------
-
 
 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
Index: external/fastjet/tools/Filter.hh
===================================================================
--- external/fastjet/tools/Filter.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/Filter.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -2,8 +2,8 @@
 #define __FASTJET_TOOLS_FILTER_HH__
 
-//STARTHEADER
-// $Id: Filter.hh 2694 2011-11-14 22:27:51Z salam $
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: Filter.hh 3494 2014-07-30 20:38:48Z soyez $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -16,7 +16,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -28,5 +30,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/ClusterSequence.hh>
@@ -98,5 +100,5 @@
   /// Note: this is just for derived classes
   ///       a Filter initialised through this constructor will not work!
-  Filter() : _Rfiltfunc(0){};
+  Filter() : _Rfiltfunc(0), _initialised(false){};
 
   /// define a filter that decomposes a jet into subjets using a
@@ -113,5 +115,5 @@
   /// ghosts
   Filter(JetDefinition subjet_def, Selector selector, double rho = 0.0) : 
-    _subjet_def(subjet_def), _Rfiltfunc(0), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0) {}
+    _subjet_def(subjet_def), _Rfiltfunc(0), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0), _initialised(true) {}
 
   /// Same as the full constructor (see above) but just specifying the radius
@@ -121,5 +123,5 @@
   ///  \param Rfilt   the filtering radius
   Filter(double Rfilt, Selector selector, double rho = 0.0) : 
-    _Rfiltfunc(0), _Rfilt(Rfilt), _selector(selector), _rho(rho), _subtractor(0) { 
+    _Rfiltfunc(0), _Rfilt(Rfilt), _selector(selector), _rho(rho), _subtractor(0), _initialised(true) { 
     if (_Rfilt<0)
       throw Error("Attempt to create a Filter with a negative filtering radius");
@@ -133,5 +135,5 @@
   ///  \param Rfilt_func   the filtering radius function of a PseudoJet
   Filter(FunctionOfPseudoJet<double> *Rfilt_func, Selector selector, double rho = 0.0) : 
-    _Rfiltfunc(Rfilt_func), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0) {}
+    _Rfiltfunc(Rfilt_func), _Rfilt(-1), _selector(selector), _rho(rho), _subtractor(0), _initialised(true) {}
 
   /// default dtor
@@ -140,5 +142,5 @@
   /// Set a subtractor that is applied to all individual subjets before
   /// deciding which ones to keep. It takes precedence over a non-zero rho.
-  void set_subtractor(const Transformer * subtractor) {_subtractor = subtractor;}
+  void set_subtractor(const FunctionOfPseudoJet<PseudoJet> * subtractor) {_subtractor = subtractor;}
 
   /// runs the filtering and sets kept and rejected to be the jets of interest
@@ -159,45 +161,20 @@
   /// It also sets the subjet_def to be used in joining things (the bit of
   /// subjet def that is of interest for later is the recombiner).
-  void _set_filtered_elements(const PseudoJet & jet,
-                              std::vector<PseudoJet> & filtered_elements,
-                              JetDefinition & subjet_def,
-                              bool & discard_area) const;
+  ///
+  /// this returns true if teh optimisation trick for C/A reclustering has been used
+  bool _set_filtered_elements(const PseudoJet & jet,
+                              std::vector<PseudoJet> & filtered_elements) const;
   
-  /// set the filtered elements in the simple case of C/A+C/A
-  void _set_filtered_elements_cafilt(const PseudoJet & jet,
-                                     std::vector<PseudoJet> & filtered_elements,
-                                     double Rfilt) const;
-
-  /// set the filtered elements in the generic re-clustering case
-  void _set_filtered_elements_generic(const PseudoJet & jet, 
-                                      std::vector<PseudoJet> & filtered_elements,
-                                      const JetDefinition & subjet_def,
-				      bool do_areas) const;
-
   /// gather the information about what is kept and rejected under the
   /// form of a PseudoJet with a special ClusterSequenceInfo
+  ///
+  /// The last argument (ca_optimisation_used) should be true if the
+  /// optimisation trick for C/A reclustering has been used (in which
+  /// case some extra tests have to be run for non-explicit-ghost
+  /// areas)
   PseudoJet _finalise(const PseudoJet & jet, 
                       std::vector<PseudoJet> & kept, 
                       std::vector<PseudoJet> & rejected,
-                      const JetDefinition & subjet_def,
-                      const bool discard_area) const;
-
-  // a series of checks
-  //--------------------------------------------------------------------
-  /// get the pieces down to the fundamental pieces
-  bool _get_all_pieces(const PseudoJet &jet, std::vector<PseudoJet> &all_pieces) const;
-
-  /// get the common recombiner to all pieces (NULL if none)
-  const JetDefinition::Recombiner* _get_common_recombiner(const std::vector<PseudoJet> &all_pieces) const;
-
-  /// check if one can apply the simplified trick for C/A subjets
-  bool _check_ca(const std::vector<PseudoJet> &all_pieces) const;
-
-  /// check if the jet (or all its pieces) have explicit ghosts
-  /// (assuming the jet has area support
-  ///
-  /// Note that if the jet has an associated cluster sequence that is no
-  /// longer valid, an error will be thrown
-  bool _check_explicit_ghosts(const std::vector<PseudoJet> &all_pieces) const;
+		      bool ca_optimisation_used) const;
 
   bool _uses_subtraction() const {return (_subtractor || _rho != 0);}
@@ -209,5 +186,7 @@
   Selector _selector;  ///< the subjet selection criterium
   double _rho;                 ///< the background density (used for subtraction when possible)
-  const Transformer * _subtractor; ///< for subtracting bkgd density from subjets
+  const FunctionOfPseudoJet<PseudoJet> * _subtractor; ///< for subtracting bkgd density from subjets
+
+  bool _initialised;    ///< true when the Filter has been properly intialised
 };
 
Index: external/fastjet/tools/GridMedianBackgroundEstimator.cc
===================================================================
--- external/fastjet/tools/GridMedianBackgroundEstimator.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/GridMedianBackgroundEstimator.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: GridMedianBackgroundEstimator.cc 3555 2014-08-11 09:56:35Z salam $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 
@@ -32,4 +34,5 @@
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
 
 //----------------------------------------------------------------------
@@ -39,15 +42,87 @@
 // of the specified particles.
 void GridMedianBackgroundEstimator::set_particles(const vector<PseudoJet> & particles) {
-  fill(_scalar_pt.begin(), _scalar_pt.end(), 0.0);
-  for (unsigned i = 0; i < particles.size(); i++) {
-    int j = igrid(particles[i]);
-    if (j >= 0){
-      if (_rescaling_class == 0)
-        _scalar_pt[j] += particles[i].perp();
-      else
-        _scalar_pt[j] += particles[i].perp()/(*_rescaling_class)(particles[i]);
+  vector<double> scalar_pt(n_tiles(), 0.0);
+
+#ifdef FASTJET_GMBGE_USEFJGRID
+  assert(all_tiles_equal_area());
+  //assert(n_good_tiles() == n_tiles()); // not needed now that we have an implementation
+#endif
+
+  // check if we need to compute only rho or both rho and rho_m
+  if (_enable_rho_m){
+    // both rho and rho_m
+    //
+    // this requires a few other variables
+    vector<double> scalar_dt(n_tiles(), 0.0);
+    double pt, dt;
+    for (unsigned i = 0; i < particles.size(); i++) {
+      int j = tile_index(particles[i]);
+      if (j >= 0){
+	pt = particles[i].pt();
+	dt = particles[i].mt() - pt;
+	if (_rescaling_class == 0){
+	  scalar_pt[j] += pt;
+	  scalar_dt[j] += dt;
+	} else {
+	  double r = (*_rescaling_class)(particles[i]);
+	  scalar_pt[j] += pt/r;
+	  scalar_dt[j] += dt/r;
+	}
+      }
+    }
+    // sort things for _percentile
+    sort(scalar_dt.begin(), scalar_dt.end());
+
+    // compute rho_m and sigma_m (see comment below for the
+    // normaliosation of sigma)
+    double p50 = _percentile(scalar_dt, 0.5);
+    _rho_m   = p50 / mean_tile_area();
+    _sigma_m = (p50-_percentile(scalar_dt, (1.0-0.6827)/2.0))/sqrt(mean_tile_area());
+  } else {
+    // only rho
+    //fill(_scalar_pt.begin(), _scalar_pt.end(), 0.0);
+    for (unsigned i = 0; i < particles.size(); i++) {
+      int j = tile_index(particles[i]);
+      if (j >= 0){
+	if (_rescaling_class == 0){
+	  scalar_pt[j] += particles[i].pt();
+	} else {
+	  scalar_pt[j] += particles[i].pt()/(*_rescaling_class)(particles[i]);
+	}
+      }
     }
   }
-  sort(_scalar_pt.begin(), _scalar_pt.end());
+
+  // if there are some "bad" tiles, then we need to exclude them from
+  // the calculation of the median. We'll do this by condensing the
+  // scalar_pt vector down to just the values for the tiles that are
+  // good.
+  //
+  // tested answers look right in "issue" 2014-08-08-testing-rect-grid
+  if (n_good_tiles() != n_tiles()) {
+    int newn = 0;
+    for (unsigned i = 0; i < scalar_pt.size(); i++) {
+      if (tile_is_good(i)) {
+        // clang gets confused with the SharedPtr swap if we don't
+        // have std:: here
+        std::swap(scalar_pt[i],scalar_pt[newn]);
+        newn++;
+      }
+    }
+    scalar_pt.resize(newn);
+  }
+
+  // in all cases, carry on with the computation of rho
+  // 
+  // first sort
+  sort(scalar_pt.begin(), scalar_pt.end());
+
+  // then compute rho
+  //
+  // watch out: by definition, our sigma is the standard deviation of
+  // the pt density multiplied by the square root of the cell area
+  double p50 = _percentile(scalar_pt, 0.5);
+  _rho   = p50 / mean_tile_area();
+  _sigma = (p50-_percentile(scalar_pt, (1.0-0.6827)/2.0))/sqrt(mean_tile_area());
 
   _has_particles = true;
@@ -61,5 +136,5 @@
 double GridMedianBackgroundEstimator::rho() const {
   verify_particles_set();
-  return _percentile(_scalar_pt, 0.5) / _cell_area;
+  return _rho;
 }
 
@@ -71,9 +146,5 @@
 double GridMedianBackgroundEstimator::sigma() const{
   verify_particles_set();
-  // watch out: by definition, our sigma is the standard deviation of
-  // the pt density multiplied by the square root of the cell area
-  return (_percentile(_scalar_pt, 0.5) -
-          _percentile(_scalar_pt, (1.0-0.6827)/2.0)
-          )/sqrt(_cell_area);
+  return _sigma; 
 }
 
@@ -85,5 +156,5 @@
 // determination.
 double GridMedianBackgroundEstimator::rho(const PseudoJet & jet)  {
-  verify_particles_set();
+  //verify_particles_set();
   double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
   return rescaling*rho();
@@ -95,7 +166,49 @@
 // the position of a given jet. As for rho(jet), it is non-const.
 double GridMedianBackgroundEstimator::sigma(const PseudoJet & jet){
-  verify_particles_set();
+  //verify_particles_set();
   double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
   return rescaling*sigma();
+}
+
+//----------------------------------------------------------------------
+// returns rho_m (particle-masses contribution to the 4-vector density)
+double GridMedianBackgroundEstimator::rho_m() const {
+  if (! _enable_rho_m){
+    throw Error("GridMediamBackgroundEstimator: rho_m requested but rho_m calculation has been disabled.");
+  }
+  verify_particles_set();
+  return _rho_m;
+}
+
+
+//----------------------------------------------------------------------
+// returns sigma_m (particle-masses contribution to the 4-vector
+// density); must be multipled by sqrt(area) to get fluctuations
+// for a region of a given area.
+double GridMedianBackgroundEstimator::sigma_m() const{
+  if (! _enable_rho_m){
+    throw Error("GridMediamBackgroundEstimator: sigma_m requested but rho_m/sigma_m calculation has been disabled.");
+  }
+  verify_particles_set();
+  return _sigma_m; 
+}
+
+//----------------------------------------------------------------------
+// returns rho_m locally at the position of a given jet. As for
+// rho(jet), it is non-const.
+double GridMedianBackgroundEstimator::rho_m(const PseudoJet & jet)  {
+  //verify_particles_set();
+  double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
+  return rescaling*rho_m();
+}
+
+
+//----------------------------------------------------------------------
+// returns sigma_m locally at the position of a given jet. As for
+// rho(jet), it is non-const.
+double GridMedianBackgroundEstimator::sigma_m(const PseudoJet & jet){
+  //verify_particles_set();
+  double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
+  return rescaling*sigma_m();
 }
 
@@ -112,6 +225,11 @@
 string GridMedianBackgroundEstimator::description() const { 
   ostringstream desc;
+#ifdef FASTJET_GMBGE_USEFJGRID
+  desc << "GridMedianBackgroundEstimator, with " << RectangularGrid::description();
+#else
   desc << "GridMedianBackgroundEstimator, with grid extension |y| < " << _ymax 
-       << " and requested grid spacing = " << _requested_grid_spacing;
+       << ", and grid cells of size dy x dphi = " << _dy << " x " << _dphi
+       << " (requested size = " << _requested_grid_spacing << ")";
+#endif
   return desc.str();
 }       
@@ -144,4 +262,5 @@
 
 
+#ifndef FASTJET_GMBGE_USEFJGRID
 //----------------------------------------------------------------------
 // protected material
@@ -168,12 +287,12 @@
 
   _ntotal = _nphi * _ny;
-  _scalar_pt.resize(_ntotal);
-  _cell_area = _dy * _dphi;
-}
-
-
-//----------------------------------------------------------------------
-// retrieve the grid cell index for a given PseudoJet
-int GridMedianBackgroundEstimator::igrid(const PseudoJet & p) const {
+  //_scalar_pt.resize(_ntotal);
+  _tile_area = _dy * _dphi;
+}
+
+
+//----------------------------------------------------------------------
+// retrieve the grid tile index for a given PseudoJet
+int GridMedianBackgroundEstimator::tile_index(const PseudoJet & p) const {
   // directly taking int does not work for values between -1 and 0
   // so use floor instead
@@ -193,8 +312,10 @@
   if (iphi == _nphi) iphi = 0; // just in case of rounding errors
 
-  int igrid_res = iy*_nphi + iphi;
-  assert (igrid_res >= 0 && igrid_res < _ny*_nphi);
-  return igrid_res;
-}
+  int index_res = iy*_nphi + iphi;
+  assert (index_res >= 0 && index_res < _ny*_nphi);
+  return index_res;
+}
+#endif // FASTJET_GMBGE_USEFJGRID
+
 
 
Index: external/fastjet/tools/GridMedianBackgroundEstimator.hh
===================================================================
--- external/fastjet/tools/GridMedianBackgroundEstimator.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/GridMedianBackgroundEstimator.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -2,8 +2,8 @@
 #define __GRID_MEDIAN_BACKGROUND_ESTIMATOR_HH__
 
-//STARTHEADER
-// $Id: GridMedianBackgroundEstimator.hh 2580 2011-09-13 17:25:43Z salam $
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: GridMedianBackgroundEstimator.hh 3610 2014-08-13 09:49:28Z salam $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -16,7 +16,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -28,8 +30,20 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 
 #include "fastjet/tools/BackgroundEstimatorBase.hh"
+
+// if defined then we'll use the RectangularGrid class
+//
+// (For FastJet 3.2, maybe remove the symbol and simply clean up the
+// code below to use exclusively the RectangularGrid)
+#define FASTJET_GMBGE_USEFJGRID
+
+#ifdef FASTJET_GMBGE_USEFJGRID
+#include "fastjet/RectangularGrid.hh"
+#endif
+
+
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
@@ -61,8 +75,14 @@
 ///   rho() [Without rescaling, they are identical]
 ///
-class GridMedianBackgroundEstimator : public BackgroundEstimatorBase {
+class GridMedianBackgroundEstimator : public BackgroundEstimatorBase
+#ifdef FASTJET_GMBGE_USEFJGRID
+                                                                    , RectangularGrid
+#endif 
+{
+
 public:
   /// @name  constructors and destructors
   //\{
+#ifdef FASTJET_GMBGE_USEFJGRID
   //----------------------------------------------------------------
   ///   \param ymax   maximal absolute rapidity extent of the grid
@@ -71,7 +91,26 @@
   ///             periodicity in azimuthal angle (size, not area)
   GridMedianBackgroundEstimator(double ymax, double requested_grid_spacing) :
+    RectangularGrid(ymax, requested_grid_spacing),
+    _has_particles(false), _enable_rho_m(true) {} 
+
+  //----------------------------------------------------------------
+  /// Constructor based on a user's fully specified RectangularGrid
+  GridMedianBackgroundEstimator(const RectangularGrid & grid) :
+    RectangularGrid(grid),
+    _has_particles(false), _enable_rho_m(true) {
+    if (!RectangularGrid::is_initialised()) 
+      throw Error("attempt to construct GridMedianBackgroundEstimator with uninitialised RectangularGrid");
+  }    
+
+#else  // alternative in old framework where we didn't have the rectangular grid
+  GridMedianBackgroundEstimator(double ymax, double requested_grid_spacing) :
     _ymin(-ymax), _ymax(ymax), 
     _requested_grid_spacing(requested_grid_spacing),
-    _has_particles(false){setup_grid();}
+    _has_particles(false), _enable_rho_m(true)
+  {
+     setup_grid();
+  }
+#endif // FASTJET_GMBGE_USEFJGRID
+
   //\}
 
@@ -84,4 +123,8 @@
   /// of the specified particles.
   void set_particles(const std::vector<PseudoJet> & particles);
+
+  /// determine whether the automatic calculation of rho_m and sigma_m
+  /// is enabled (by default true)
+  void set_compute_rho_m(bool enable){ _enable_rho_m = enable;}
 
   //\}
@@ -114,7 +157,32 @@
   bool has_sigma() {return true;}
 
+  //-----------------------------------------------------------------
+  /// Returns rho_m, the purely longitudinal, particle-mass-induced
+  /// component of the background density per unit area
+  double rho_m() const;
+
+  /// returns sigma_m, a measure of the fluctuations in the purely
+  /// longitudinal, particle-mass-induced component of the background
+  /// density per unit area; must be multipled by sqrt(area) to get
+  /// fluctuations for a region of a given area.
+  double sigma_m() const;
+
+  /// Returns rho_m locally at the jet position. As for rho(jet), it is non-const.
+  double rho_m(const PseudoJet & jet);
+
+  /// Returns sigma_m locally at the jet position. As for rho(jet), it is non-const.
+  double sigma_m(const PseudoJet & jet);
+
+  /// Returns true if this background estimator has support for
+  /// determination of rho_m.
+  ///
+  /// Note that support for sigma_m is automatic is one has sigma and
+  /// rho_m support.
+  bool has_rho_m() const {return _enable_rho_m;}
+
+
   /// returns the area of the grid cells (all identical, but
   /// referred to as "mean" area for uniformity with JetMedianBGE).
-  double mean_area() const {return _cell_area;}
+  double mean_area() const {return mean_tile_area();}
   //\}
 
@@ -134,4 +202,8 @@
   /// Note that this has to be called BEFORE any attempt to do an
   /// actual computation
+  ///
+  /// The same profile will be used for both pt and mt (this is
+  /// probabaly a good approximation since the particle density
+  /// changes is what dominates the rapidity profile)
   virtual void set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class);
 
@@ -149,22 +221,35 @@
 
 private:
+
+#ifndef FASTJET_GMBGE_USEFJGRID
+
   /// configure the grid
   void setup_grid();
 
   /// retrieve the grid cell index for a given PseudoJet
-  int igrid(const PseudoJet & p) const;
+  int tile_index(const PseudoJet & p) const;
+
+  // information about the grid
+  double _ymin, _ymax, _dy, _dphi, _requested_grid_spacing, _tile_area;
+  int _ny, _nphi, _ntotal;
+
+  int n_tiles() const {return _ntotal;}
+  int n_good_tiles() const {return n_tiles();}
+  int tile_is_good(int /* itile */) const {return true;}
+
+  double mean_tile_area() const {return _tile_area;}
+#endif // FASTJET_GMBGE_USEFJGRID
+
 
   /// verify that particles have been set and throw an error if not
   void verify_particles_set() const;
 
-  // information about the grid
-  double _ymin, _ymax, _dy, _dphi, _requested_grid_spacing, _cell_area;
-  int _ny, _nphi, _ntotal;
-
   // information abotu the event
-  std::vector<double> _scalar_pt;
+  //std::vector<double> _scalar_pt;
+  double _rho, _sigma, _rho_m, _sigma_m;
   bool _has_particles;
-
-  // various warnings to let people aware of potential dangers
+  bool _enable_rho_m;
+
+  // various warnings to inform people of potential dangers
   LimitedWarning _warning_rho_of_jet;
   LimitedWarning _warning_rescaling;
Index: external/fastjet/tools/JHTopTagger.cc
===================================================================
--- external/fastjet/tools/JHTopTagger.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/JHTopTagger.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: JHTopTagger.cc 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/tools/JHTopTagger.hh>
Index: external/fastjet/tools/JHTopTagger.hh
===================================================================
--- external/fastjet/tools/JHTopTagger.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/JHTopTagger.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -2,8 +2,8 @@
 #define __FASTJET_JH_TOP_TAGGER_HH__
 
-//STARTHEADER
-// $Id: JHTopTagger.hh 2689 2011-11-14 14:51:06Z soyez $
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: JHTopTagger.hh 3433 2014-07-23 08:17:03Z salam $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -16,7 +16,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -28,5 +30,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 
Index: external/fastjet/tools/JetMedianBackgroundEstimator.cc
===================================================================
--- external/fastjet/tools/JetMedianBackgroundEstimator.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/JetMedianBackgroundEstimator.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: JetMedianBackgroundEstimator.cc 3517 2014-08-01 14:23:13Z soyez $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
@@ -31,4 +33,5 @@
 #include <fastjet/ClusterSequenceStructure.hh>
 #include <iostream>
+#include <sstream>
 
 FASTJET_BEGIN_NAMESPACE     // defined in fastjet/internal/base.hh
@@ -37,5 +40,7 @@
 
 double BackgroundJetScalarPtDensity::result(const PseudoJet & jet) const {
-  std::vector<PseudoJet> constituents = jet.constituents();
+  // do not include the ghosts in the list of constituents to have a
+  // correct behaviour when _pt_power is <= 0
+  std::vector<PseudoJet> constituents = (!SelectorIsPureGhost())(jet.constituents());
   double scalar_pt = 0;
   for (unsigned i = 0; i < constituents.size(); i++) {
@@ -43,4 +48,12 @@
   }
   return scalar_pt / jet.area();
+}
+
+
+std::string BackgroundJetScalarPtDensity::description() const {
+  ostringstream oss;
+  oss << "BackgroundScalarJetPtDensity";
+  if (_pt_power != 1.0) oss << " with pt_power = " << _pt_power;
+  return oss.str();
 }
 
@@ -72,5 +85,5 @@
                                          const JetDefinition &jet_def,
                                          const AreaDefinition &area_def)
-  : _rho_range(rho_range), _jet_def(jet_def), _area_def(area_def) {
+  : _rho_range(rho_range), _jet_def(jet_def), _area_def(area_def){
 
   // initialise things decently
@@ -261,4 +274,56 @@
 
 //----------------------------------------------------------------------
+// returns rho_m (particle-masses contribution to the 4-vector density)
+double JetMedianBackgroundEstimator::rho_m() const {
+  if (! has_rho_m()){
+    throw Error("JetMediamBackgroundEstimator: rho_m requested but rho_m calculation is disabled (either eplicitly or due to the presence of a jet density class).");
+  }
+  if (_rho_range.takes_reference())
+    throw Error("The background estimation is obtained from a selector that takes a reference jet. rho(PseudoJet) should be used in that case");
+  _recompute_if_needed();
+  return _rho_m;
+}
+
+
+//----------------------------------------------------------------------
+// returns sigma_m (particle-masses contribution to the 4-vector
+// density); must be multipled by sqrt(area) to get fluctuations
+// for a region of a given area.
+double JetMedianBackgroundEstimator::sigma_m() const{
+  if (! has_rho_m()){
+    throw Error("JetMediamBackgroundEstimator: sigma_m requested but rho_m/sigma_m calculation is disabled (either explicitly or due to the presence of a jet density class).");
+  }
+  if (_rho_range.takes_reference())
+    throw Error("The background estimation is obtained from a selector that takes a reference jet. rho(PseudoJet) should be used in that case");
+  _recompute_if_needed();
+  return _sigma_m;
+}
+
+//----------------------------------------------------------------------
+// returns rho_m locally at the position of a given jet. As for
+// rho(jet), it is non-const.
+double JetMedianBackgroundEstimator::rho_m(const PseudoJet & jet)  {
+  _recompute_if_needed(jet);
+  double our_rho = _rho_m;
+  if (_rescaling_class != 0) { 
+    our_rho *= (*_rescaling_class)(jet);
+  }
+  return our_rho;
+}
+
+
+//----------------------------------------------------------------------
+// returns sigma_m locally at the position of a given jet. As for
+// rho(jet), it is non-const.
+double JetMedianBackgroundEstimator::sigma_m(const PseudoJet & jet){
+  _recompute_if_needed(jet);
+  double our_sigma = _sigma_m;
+  if (_rescaling_class != 0) { 
+    our_sigma *= (*_rescaling_class)(jet);
+  }
+  return our_sigma;
+}
+
+//----------------------------------------------------------------------
 // configuring behaviour
 //----------------------------------------------------------------------
@@ -271,6 +336,9 @@
   set_provide_fj2_sigma(false);
 
+  _enable_rho_m = true;
+
   // reset the computed values
   _rho = _sigma = 0.0;
+  _rho_m = _sigma_m = 0.0;
   _n_jets_used = _n_empty_jets = 0;
   _empty_area = _mean_area = 0.0;
@@ -289,5 +357,5 @@
 // is used (as occurs also if this function is not called).
 void JetMedianBackgroundEstimator::set_jet_density_class(const FunctionOfPseudoJet<double> * jet_density_class_in) {
-  _warnings_preliminary.warn("JetMedianBackgroundEstimator::set_jet_density_class: density classes are still preliminary in FastJet 3.0. Their interface may differ in future releases (without guaranteeing backward compatibility).");
+  _warnings_preliminary.warn("JetMedianBackgroundEstimator::set_jet_density_class: density classes are still preliminary in FastJet 3.1. Their interface may differ in future releases (without guaranteeing backward compatibility). Note that since FastJet 3.1, rho_m and sigma_m are accessible direclty in JetMedianBackgroundEstimator and GridMedianBackgroundEstimator(with no need for a density class).");
   _jet_density_class = jet_density_class_in;
   _uptodate = false;
@@ -338,5 +406,6 @@
   // fill the vector of pt/area (or the quantity from the jet density class) 
   //  - in the range
-  vector<double> vector_for_median;
+  vector<double> vector_for_median_pt;
+  vector<double> vector_for_median_dt;
   double total_area  = 0.0;
   _n_jets_used = 0;
@@ -346,20 +415,38 @@
 
   // compute the pt/area for the selected jets
+  double median_input_pt, median_input_dt=0.0;
+  BackgroundJetPtMDensity m_density;
+  bool do_rho_m = has_rho_m();
   for (unsigned i = 0; i < selected_jets.size(); i++) {
     const PseudoJet & current_jet = selected_jets[i];
 
     double this_area = (_use_area_4vector) ? current_jet.area_4vector().perp() : current_jet.area(); 
-
     if (this_area>0){
-      double median_input;
+      // for the pt component, we either use pt or the user-provided
+      // density class
       if (_jet_density_class == 0) {
-        median_input = current_jet.perp()/this_area;
+        median_input_pt = current_jet.perp()/this_area;
       } else {
-        median_input = (*_jet_density_class)(current_jet);
+        median_input_pt = (*_jet_density_class)(current_jet);
       }
+
+      // handle the rho_m part if requested
+      // note that we're using the scalar area as a normalisation inside the
+      // density class!
+      if (do_rho_m) 
+	median_input_dt = m_density(current_jet);
+    
+      // perform rescaling if needed
       if (_rescaling_class != 0) {
-        median_input /= (*_rescaling_class)(current_jet);
+	double resc = (*_rescaling_class)(current_jet);;
+        median_input_pt /= resc;
+        median_input_dt /= resc;
       }
-      vector_for_median.push_back(median_input);
+      
+      // store the result for future computation of the median
+      vector_for_median_pt.push_back(median_input_pt);
+      if (do_rho_m) 
+	vector_for_median_dt.push_back(median_input_dt);
+
       total_area  += this_area;
       _n_jets_used++;
@@ -367,11 +454,12 @@
       _warnings_zero_area.warn("JetMedianBackgroundEstimator::_compute(...): discarded jet with zero area. Zero-area jets may be due to (i) too large a ghost area (ii) a jet being outside the ghost range (iii) the computation not being done using an appropriate algorithm (kt;C/A).");
     }
-      
   }
   
   // there is nothing inside our region, so answer will always be zero
-  if (vector_for_median.size() == 0) {
+  if (vector_for_median_pt.size() == 0) {
     _rho        = 0.0;
     _sigma      = 0.0;
+    _rho_m      = 0.0;
+    _sigma_m    = 0.0;
     _mean_area  = 0.0;
     return;
@@ -392,5 +480,5 @@
 
   double stand_dev;
-  _median_and_stddev(vector_for_median, _n_empty_jets, _rho, stand_dev, 
+  _median_and_stddev(vector_for_median_pt, _n_empty_jets, _rho, stand_dev, 
                      _provide_fj2_sigma);
 
@@ -398,4 +486,11 @@
   _mean_area  = total_area / total_njets;
   _sigma      = stand_dev * sqrt(_mean_area);
+
+  // compute the rho_m part now
+  if (do_rho_m){
+    _median_and_stddev(vector_for_median_dt, _n_empty_jets, _rho_m, stand_dev, 
+		       _provide_fj2_sigma);
+    _sigma_m = stand_dev * sqrt(_mean_area);
+  }
 
   // record that the computation has been performed  
@@ -439,5 +534,4 @@
 
 
-
 FASTJET_END_NAMESPACE
 
Index: external/fastjet/tools/JetMedianBackgroundEstimator.hh
===================================================================
--- external/fastjet/tools/JetMedianBackgroundEstimator.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/JetMedianBackgroundEstimator.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -2,8 +2,8 @@
 #define __FASTJET_BACKGROUND_ESTIMATOR_HH__
 
-//STARTHEADER
-// $Id: JetMedianBackgroundEstimator.hh 2689 2011-11-14 14:51:06Z soyez $
+//FJSTARTHEADER
+// $Id: JetMedianBackgroundEstimator.hh 3517 2014-08-01 14:23:13Z soyez $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -16,7 +16,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -28,5 +30,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/ClusterSequenceAreaBase.hh>
@@ -121,5 +123,6 @@
   ///
   JetMedianBackgroundEstimator(const Selector &rho_range = SelectorIdentity())
-    : _rho_range(rho_range), _jet_def(JetDefinition()) { reset(); }
+    : _rho_range(rho_range), _jet_def(JetDefinition()),
+      _enable_rho_m(true){ reset(); }
   
 
@@ -168,4 +171,8 @@
   }
 
+  /// determine whether the automatic calculation of rho_m and sigma_m
+  /// is enabled (by default true)
+  void set_compute_rho_m(bool enable){ _enable_rho_m = enable;}
+
   //\}
 
@@ -201,4 +208,32 @@
   virtual bool has_sigma() {return true;}
 
+  //----------------------------------------------------------------
+  // now do the same thing for rho_m and sigma_m
+
+  /// returns rho_m, the purely longitudinal, particle-mass-induced
+  /// component of the background density per unit area
+  virtual double rho_m() const;
+
+  /// returns sigma_m, a measure of the fluctuations in the purely
+  /// longitudinal, particle-mass-induced component of the background
+  /// density per unit area; must be multipled by sqrt(area) to get
+  /// fluctuations for a region of a given area.
+  virtual double sigma_m() const;
+
+  /// Returns rho_m locally at the jet position. As for rho(jet), it is non-const.
+  virtual double rho_m(const PseudoJet & /*jet*/);
+
+  /// Returns sigma_m locally at the jet position. As for rho(jet), it is non-const.
+  virtual double sigma_m(const PseudoJet & /*jet*/);
+
+  /// Returns true if this background estimator has support for
+  /// determination of rho_m.
+  ///
+  /// In te presence of a density class, support for rho_m is
+  /// automatically disabled
+  ///
+  /// Note that support for sigma_m is automatic is one has sigma and
+  /// rho_m support.
+  virtual bool has_rho_m() const {return _enable_rho_m && (_jet_density_class == 0);}
   //\}
   
@@ -208,6 +243,9 @@
   /// Returns the mean area of the jets used to actually compute the
   /// background properties in the last call of rho() or sigma()
+  /// If the configuration has changed in the meantime, throw an error.
   double mean_area() const{
-    _recompute_if_needed();
+    if (!_uptodate)
+      throw Error("JetMedianBackgroundEstimator::mean_area(): one may not retrieve information about the last call to rho() or sigma() when the configuration has changed in the meantime.");
+    //_recompute_if_needed();
     return _mean_area;
   }
@@ -215,7 +253,23 @@
   /// returns the number of jets used to actually compute the
   /// background properties in the last call of rho() or sigma()
+  /// If the configuration has changed in the meantime, throw an error.
   unsigned int n_jets_used() const{
-    _recompute_if_needed();
+    if (!_uptodate)
+      throw Error("JetMedianBackgroundEstimator::n_jets_used(): one may not retrieve information about the last call to rho() or sigma() when the configuration has changed in the meantime.");
+    //_recompute_if_needed();
     return _n_jets_used;
+  }
+
+  /// returns the jets used to actually compute the background
+  /// properties
+  std::vector<PseudoJet> jets_used() const{
+    if (!_uptodate) throw Error("JetMedianBackgroundEstimator::n_jets_used(): one may not retrieve information about the last call to rho() or sigma() when the configuration has changed in the meantime.");
+    _check_csa_alive();
+    std::vector<PseudoJet> tmp_jets = _rho_range(_included_jets);
+    std::vector<PseudoJet> used_jets;
+    for (unsigned int i=0; i<tmp_jets.size(); i++){
+      if (tmp_jets[i].area()>0) used_jets.push_back(tmp_jets[i]);
+    }
+    return used_jets;
   }
 
@@ -223,4 +277,5 @@
   /// the selector) that is not occupied by jets. The value is that
   /// for the last call of rho() or sigma()
+  /// If the configuration has changed in the meantime, throw an error.
   ///
   /// The answer is defined to be zero if the area calculation
@@ -234,5 +289,7 @@
   /// call to the ClusterSequenceAreaBase function.
   double empty_area() const{
-    _recompute_if_needed();
+    if (!_uptodate)
+      throw Error("JetMedianBackgroundEstimator::empty_area(): one may not retrieve information about the last call to rho() or sigma() when the configuration has changed in the meantime.");
+    //_recompute_if_needed();
     return _empty_area;
   }
@@ -241,4 +298,5 @@
   /// background properties. The value is that for the last call of
   /// rho() or sigma().
+  /// If the configuration has changed in the meantime, throw an error.
   ///
   /// If the area has explicit ghosts the result is zero; for active
@@ -250,5 +308,7 @@
   /// call to the ClusterSequenceAreaBase function.
   double n_empty_jets() const{
-    _recompute_if_needed();
+    if (!_uptodate)
+      throw Error("JetMedianBackgroundEstimator::n_empty_jets(): one may not retrieve information about the last call to rho() or sigma() when the configuration has changed in the meantime.");
+    //_recompute_if_needed();
     return _n_empty_jets;
   }
@@ -361,5 +421,5 @@
   /// Issue a warning otherwise
   void _check_jet_alg_good_for_median() const;
-  
+
   // the basic parameters of this class (passed through the variou ctors)
   Selector _rho_range;                   ///< range to compute the background in
@@ -368,13 +428,16 @@
   std::vector<PseudoJet> _included_jets; ///< jets to be used
   
-  // the tunable aprameters of the class
+  // the tunable parameters of the class
   bool _use_area_4vector;
   bool _provide_fj2_sigma;
   const FunctionOfPseudoJet<double> * _jet_density_class;
   //SharedPtr<BackgroundRescalingBase> _rescaling_class_sharedptr;
+  bool _enable_rho_m;
   
   // the actual results of the computation
   mutable double _rho;               ///< background estimated density per unit area
   mutable double _sigma;             ///< background estimated fluctuations
+  mutable double _rho_m;             ///< "mass" background estimated density per unit area
+  mutable double _sigma_m;           ///< "mass" background estimated fluctuations
   mutable double _mean_area;         ///< mean area of the jets used to estimate the background
   mutable unsigned int _n_jets_used; ///< number of jets used to estimate the background
@@ -429,5 +492,5 @@
   virtual double result(const PseudoJet & jet) const;
 
-  virtual std::string description() const {return "BackgroundScalarJetPtDensity";}
+  virtual std::string description() const;
 
 private:
Index: external/fastjet/tools/MassDropTagger.cc
===================================================================
--- external/fastjet/tools/MassDropTagger.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/MassDropTagger.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: MassDropTagger.cc 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/tools/MassDropTagger.hh>
Index: external/fastjet/tools/MassDropTagger.hh
===================================================================
--- external/fastjet/tools/MassDropTagger.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/MassDropTagger.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id: MassDropTagger.hh 2731 2011-11-21 12:15:21Z soyez $
+//FJSTARTHEADER
+// $Id: MassDropTagger.hh 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #ifndef __FASTJET_MASS_DROP_TAGGER_HH__
Index: external/fastjet/tools/Pruner.cc
===================================================================
--- external/fastjet/tools/Pruner.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/Pruner.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: Pruner.cc 3481 2014-07-29 17:24:12Z soyez $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include "fastjet/tools/Pruner.hh"
@@ -47,10 +49,10 @@
 //----------------------------------------------------------------------
 // alternative (dynamic) ctor
-//  \param jet_def the jet definition for the internal clustering
+//  \param jet_def     the jet definition for the internal clustering
 //  \param zcut_dyn    dynamic pt-fraction cut in the pruning
 //  \param Rcut_dyn    dynamic angular distance cut in the pruning
 Pruner::Pruner(const JetDefinition &jet_def, 
-         FunctionOfPseudoJet<double> *zcut_dyn,
-         FunctionOfPseudoJet<double> *Rcut_dyn)
+         const FunctionOfPseudoJet<double> *zcut_dyn,
+         const FunctionOfPseudoJet<double> *Rcut_dyn)
   : _jet_def(jet_def), _zcut(0), _Rcut_factor(0),
     _zcut_dyn(zcut_dyn), _Rcut_dyn(Rcut_dyn), _get_recombiner_from_jet(false)  {
@@ -74,25 +76,17 @@
   double zcut = (_zcut_dyn) ? (*_zcut_dyn)(jet) : _zcut;
   PruningPlugin * pruning_plugin;
+
   // for some constructors, we get the recombiner from the 
-  // input jet -- some acrobatics are needed (see plans for FJ3.1
-  // for a hopefully better solution).
+  // input jet -- some acrobatics are needed
   if (_get_recombiner_from_jet) {
-    const JetDefinition::Recombiner * common_recombiner = 
-                                              _get_common_recombiner(jet);
-    if (common_recombiner) {
-      JetDefinition jet_def = _jet_def;
-      if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
-        RecombinationScheme scheme = 
-          static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
-        jet_def.set_recombination_scheme(scheme);
-      } else {
-        jet_def.set_recombiner(common_recombiner);
-      }
-      pruning_plugin = new PruningPlugin(jet_def, zcut, Rcut);
-    } else {
-      // if there wasn't a common recombiner, we just use the default
-      // recombiner that was in _jet_def
-      pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
-    }
+    JetDefinition jet_def = _jet_def;
+
+    // if all the pieces have a shared recombiner, we'll use that
+    // one. Otherwise, use the one from _jet_def as a fallback.
+    JetDefinition jet_def_for_recombiner;
+    if (_check_common_recombiner(jet, jet_def_for_recombiner)){
+      jet_def.set_recombiner(jet_def_for_recombiner);
+    }
+    pruning_plugin = new PruningPlugin(jet_def, zcut, Rcut);
   } else {
     pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
@@ -122,4 +116,6 @@
   PseudoJet result_local = SelectorNHardest(1)(cs->inclusive_jets())[0];
   PrunerStructure * s = new PrunerStructure(result_local);
+  s->_Rcut = Rcut;
+  s->_zcut = zcut;
   result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
   
@@ -155,31 +151,36 @@
 }
 
-// see if there is a common recombiner among the pieces; if there
-// is return a pointer to it; otherwise, return NULL.
-// 
-// NB: this way of doing things is not ideal, because quite some work
-//     is needed to get a correct handling of the final recombiner
-//     (e.g. default v. non-default). In future add
-//     set_recombiner(jet_def) to JetDefinition, maybe also add
-//     an invalid_scheme to the default recombiner and then 
-//     do all the work below directly with a JetDefinition directly
-//     together with JD::has_same_recombiner(...)
-const JetDefinition::Recombiner * Pruner::_get_common_recombiner(const PseudoJet &jet) const{
-  if (jet.has_associated_cluster_sequence())
-    return jet.validated_cs()->jet_def().recombiner();
+// see if there is a common recombiner among the pieces; if there is
+// return true and set jet_def_for_recombiner so that the recombiner
+// can be taken from that JetDefinition. Otherwise, return
+// false. 'assigned' is initially false; when true, each time we meet
+// a new jet definition, we'll check it shares the same recombiner as
+// jet_def_for_recombiner.
+bool Pruner::_check_common_recombiner(const PseudoJet &jet, 
+				      JetDefinition &jet_def_for_recombiner,
+				      bool assigned) const{
+  if (jet.has_associated_cluster_sequence()){
+    // if the jet def for recombination has already been assigned, check if we have the same
+    if (assigned)
+      return jet.validated_cs()->jet_def().has_same_recombiner(jet_def_for_recombiner);
+
+    // otherwise, assign it.
+    jet_def_for_recombiner = jet.validated_cs()->jet_def();
+    assigned = true;
+    return true;
+  }
 
   // if the jet has pieces, recurse in the pieces
   if (jet.has_pieces()){
     vector<PseudoJet> pieces = jet.pieces();
-    if (pieces.size() == 0) return 0;
-    const JetDefinition::Recombiner * reco = _get_common_recombiner(pieces[0]);
-    for (unsigned int i=1;i<pieces.size(); i++)
-      if (_get_common_recombiner(pieces[i]) != reco) return 0;
+    if (pieces.size() == 0) return false;
+    for (unsigned int i=0;i<pieces.size(); i++)
+      if (!_check_common_recombiner(pieces[i], jet_def_for_recombiner, assigned)) return false;
     // never returned false, so we're OK.
-    return reco;
+    return true;
   }
 
   // return false for any other (unknown) structure
-  return 0;
+  return false;
 }
 
Index: external/fastjet/tools/Pruner.hh
===================================================================
--- external/fastjet/tools/Pruner.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/Pruner.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -2,8 +2,8 @@
 #define __FASTJET_TOOLS_PRUNER_HH__
 
-//STARTHEADER
-// $Id: Pruner.hh 2616 2011-09-30 18:03:40Z salam $
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: Pruner.hh 3481 2014-07-29 17:24:12Z soyez $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -16,7 +16,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -28,5 +30,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include "fastjet/ClusterSequence.hh"
@@ -43,4 +45,9 @@
 class PruningRecombiner;
 class PruningPlugin;
+
+// This tells third-party code that the pruner structure 
+// stores Rcut info; the alternative is for the user to 
+// get the information from the version number
+#define FASTJET_PRUNER_STRUCTURE_STORES_RCUT
 
 //----------------------------------------------------------------------
@@ -135,6 +142,6 @@
   ///  \param Rcut_dyn    dynamic angular distance cut in the pruning
   Pruner(const JetDefinition &jet_def, 
-         FunctionOfPseudoJet<double> *zcut_dyn,
-         FunctionOfPseudoJet<double> *Rcut_dyn);
+         const FunctionOfPseudoJet<double> *zcut_dyn,
+         const FunctionOfPseudoJet<double> *Rcut_dyn);
 
   /// action on a single jet
@@ -152,13 +159,19 @@
   bool _check_explicit_ghosts(const PseudoJet &jet) const;
 
-  /// return a pointer to a "common" recombiner if there is one,
-  /// alternatively a null pointer.
-  const JetDefinition::Recombiner * _get_common_recombiner(const PseudoJet &jet) const;
+  /// see if there is a common recombiner among the pieces; if there
+  /// is return true and set jet_def_for_recombiner so that the
+  /// recombiner can be taken from that JetDefinition. Otherwise,
+  /// return false. 'assigned' is initially false; when true, each
+  /// time we meet a new jet definition, we'll check it shares the
+  /// same recombiner as jet_def_for_recombiner.
+  bool _check_common_recombiner(const PseudoJet &jet, 
+				JetDefinition &jet_def_for_recombiner,
+				bool assigned=false) const;
 
   JetDefinition _jet_def; ///< the internal jet definition
   double _zcut;        	  ///< the pt-fraction cut
   double _Rcut_factor;    ///< the angular separation cut factor
-  FunctionOfPseudoJet<double> *_zcut_dyn; ///< dynamic zcut
-  FunctionOfPseudoJet<double> *_Rcut_dyn; ///< dynamic Rcut
+  const FunctionOfPseudoJet<double> *_zcut_dyn; ///< dynamic zcut
+  const FunctionOfPseudoJet<double> *_Rcut_dyn; ///< dynamic Rcut
   bool   _get_recombiner_from_jet; ///< true for minimal constructor,
                                    ///< causes recombiner to be set equal 
@@ -194,6 +207,15 @@
   std::vector<PseudoJet> extra_jets() const;
 
+  /// return the value of Rcut that was used for this specific pruning.
+  double Rcut() const {return _Rcut;}
+
+  /// return the value of Rcut that was used for this specific pruning.
+  double zcut() const {return _zcut;}
+
 protected:
   friend class Pruner; ///< to allow setting the internal information
+
+private:
+  double _Rcut, _zcut;
 };
 
Index: external/fastjet/tools/Recluster.cc
===================================================================
--- external/fastjet/tools/Recluster.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/tools/Recluster.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,417 @@
+// $Id: Recluster.cc 3629 2014-08-14 17:21:15Z salam $
+//
+// Copyright (c) 2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+
+#include "fastjet/tools/Recluster.hh"
+#include "fastjet/CompositeJetStructure.hh"
+#include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
+#include <sstream>
+#include <typeinfo>
+
+using namespace std;
+
+// Comments:
+//
+//  - If the jet comes from a C/A clustering (or is a composite jet
+//    made of C/A clusterings) and we recluster it with a C/A
+//    algorithm, we just use exclusive jets instead of doing the
+//    clustering explicitly. In this specific case, we need to make
+//    sure that all the pieces share the same cluster sequence.
+//
+//  - If the recombiner has to be taken from the original jet and that
+//    jet is composite, we need to check that all the pieces were
+//    obtained with the same recombiner.
+//
+//  - Note that a preliminary version of this code has been
+//    distributed in the RecursiveTools fastjet-contrib. The present
+//    version has a few minor fixed
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+LimitedWarning Recluster::_explicit_ghost_warning;
+
+// ctor
+Recluster::Recluster(JetAlgorithm new_jet_alg, double new_jet_radius, Keep keep_in)
+  : _new_jet_def(JetDefinition(new_jet_alg, new_jet_radius)), 
+    _acquire_recombiner(true), _keep(keep_in), _cambridge_optimisation_enabled(true){}
+
+Recluster::Recluster(JetAlgorithm new_jet_alg, Keep keep_in)
+    : _acquire_recombiner(true), _keep(keep_in), _cambridge_optimisation_enabled(true){
+  switch (JetDefinition::n_parameters_for_algorithm(new_jet_alg)){
+  case 0: _new_jet_def = JetDefinition(new_jet_alg); break;
+  case 1: _new_jet_def = JetDefinition(new_jet_alg, JetDefinition::max_allowable_R); break;
+  default:
+    throw Error("Recluster(): tried to construct specifying only a jet algorithm ("+JetDefinition::algorithm_description(new_jet_alg)+") which takes more than 1 parameter");
+  };
+}
+
+
+// class description
+string Recluster::description() const {
+  ostringstream ostr;
+  ostr << "Recluster with new_jet_def = ";
+  if (_acquire_recombiner){
+    ostr << _new_jet_def.description_no_recombiner();
+    ostr << ", using a recombiner obtained from the jet being reclustered";
+  } else {
+    ostr << _new_jet_def.description();
+  }
+
+  if (_keep == keep_only_hardest)
+    ostr << " and keeping the hardest inclusive jet";
+  else
+    ostr << " and joining all inclusive jets into a composite jet";
+
+  return ostr.str();
+}
+
+
+// the main piece of code that performs the reclustering
+PseudoJet Recluster::result(const PseudoJet &jet) const {
+  // get the incljets and the exact jet definition that has been used
+  // to get them
+  vector<PseudoJet> incljets;
+  bool ca_optimised = get_new_jets_and_def(jet, incljets);
+
+  return generate_output_jet(incljets, ca_optimised);
+}
+
+
+// a lower-level method that does the actual work of reclustering the
+// input jet. The resulting incljets are stored in output_jets and the
+// jet definition that has been used can be deduced from their
+// associated ClusterSequence
+//
+//  - input_jet       the (input) jet that one wants to recluster
+//  - output_jets     incljets resulting from the new clustering
+//
+// returns true if the C/A optimisation has been used (this means
+// that generate_output_jet will watch out for non-explicit-ghost
+// areas that might be leftover)
+bool Recluster::get_new_jets_and_def(const PseudoJet & input_jet, 
+                                     vector<PseudoJet> & output_jets) const{
+  // generic sanity checks
+  //-------------------------------------------------------------------
+  // make sure that the jet has constituents
+  if (! input_jet.has_constituents())
+    throw Error("Recluster can only be applied on jets having constituents");
+
+  // tests particular to certain configurations
+  //-------------------------------------------------------------------
+
+  // for a whole variety of tests, we shall need the "recursive"
+  // pieces of the jet (the jet itself or recursing down to its most
+  // fundamental pieces). So we do compute these once and for all.
+  //
+  // Note that the pieces are always needed (either for C/A or for the
+  // area checks)
+  vector<PseudoJet> all_pieces;
+  if ((!_get_all_pieces(input_jet, all_pieces)) || (all_pieces.size()==0)){
+    throw Error("Recluster: failed to retrieve all the pieces composing the jet.");
+  }
+
+  // decide which jet definition to use
+  //-------------------------------------------------------------------
+  JetDefinition new_jet_def = _new_jet_def;
+  if (_acquire_recombiner){
+    _acquire_recombiner_from_pieces(all_pieces, new_jet_def);
+  }
+
+  // the vector that will ultimately hold the incljets
+  output_jets.clear();
+
+  // check if we can apply the simplification for C/A jets reclustered
+  // with C/A
+  //
+  // we apply C/A clustering iff
+  //  - the requested new_jet_def is C/A
+  //  - the jet is either directly coming from C/A or if it is a
+  //    superposition of C/A jets from the same cluster sequence
+  //  - the pieces agree with the recombination scheme of new_jet_def
+  //
+  // In this case area support will be automatically inherited so we
+  // can only worry about this later
+  // -------------------------------------------------------------------
+  if (_check_ca(all_pieces, new_jet_def)){
+    _recluster_ca(all_pieces, output_jets, new_jet_def.R());
+    output_jets = sorted_by_pt(output_jets);
+    return true;
+  }
+
+  // decide if area support has to be kept
+  //-------------------------------------------------------------------
+  bool include_area_support = input_jet.has_area();
+  if ((include_area_support) &&  (!_check_explicit_ghosts(all_pieces))){
+    _explicit_ghost_warning.warn("Recluster: the original cluster sequence is lacking explicit ghosts; area support will no longer be available after re-clustering");
+    include_area_support = false;
+  }
+
+  // extract the incljets
+  //-------------------------------------------------------------------
+  _recluster_generic(input_jet, output_jets, new_jet_def, include_area_support);
+  output_jets = sorted_by_pt(output_jets);
+
+  return false;
+}
+
+// given a set of incljets and a jet definition used, create the
+// resulting PseudoJet
+PseudoJet Recluster::generate_output_jet(std::vector<PseudoJet> & incljets, 
+                                         bool ca_optimisation_used) const{
+  // first handle the case where we only need to keep the hardest incljet
+  if (_keep == keep_only_hardest) {
+    if (incljets.size() > 0) {
+      return incljets[0];
+    } else {
+      return PseudoJet();
+    }
+  }
+
+  // now the case where all incljets have to be joined
+
+  // safekeeper
+  if (incljets.size()==0) return join(incljets);
+
+  PseudoJet reclustered = join(incljets, 
+			       *(incljets[0].associated_cluster_sequence()->jet_def().recombiner()));
+
+  // if we've used C/A optimisation, we need to get rid of the area
+  // information if it comes from a non-explicit-ghost clustering.
+  // (because in that case it can be erroneous due to the lack of
+  // information about empty areas)
+  if (ca_optimisation_used){
+    if (reclustered.has_area() &&
+        (incljets.size() > 0) &&
+        (! incljets[0].validated_csab()->has_explicit_ghosts())){
+      CompositeJetStructure *css = dynamic_cast<CompositeJetStructure *>(reclustered.structure_non_const_ptr());
+      assert(css);
+      css->discard_area();
+    }
+  }
+  
+  return reclustered;
+}
+
+//----------------------------------------------------------------------
+// the parts that really do the reclustering
+//----------------------------------------------------------------------
+
+// get the subjets in the simple case of C/A+C/A
+void Recluster::_recluster_ca(const vector<PseudoJet> & all_pieces, 
+                              vector<PseudoJet> & subjets, 
+                              double Rfilt) const{
+  subjets.clear();
+
+  // each individual piece should have a C/A cluster sequence
+  // associated to it. So a simple loop would do the job
+  for (vector<PseudoJet>::const_iterator piece_it = all_pieces.begin(); 
+       piece_it!=all_pieces.end(); piece_it++){
+    // just extract the exclusive subjets of 'jet'
+    const ClusterSequence *cs = piece_it->associated_cluster_sequence(); 
+    vector<PseudoJet> local_subjets;
+
+    double dcut = Rfilt / cs->jet_def().R();
+    if (dcut>=1.0){
+      // remember that in this case all the pairwise interpiece
+      // distances are supposed to be larger than Rfilt (this was
+      // tested in _check_ca), which means that they can never
+      // recombine with each other.
+      local_subjets.push_back(*piece_it);
+    } else {
+      local_subjets = piece_it->exclusive_subjets(dcut*dcut);
+    }
+
+    copy(local_subjets.begin(), local_subjets.end(), back_inserter(subjets));
+  }
+}
+
+
+// perform the reclustering itself for all cases where the "C/A trick"
+// does not apply
+void Recluster::_recluster_generic(const PseudoJet & jet, 
+                                   vector<PseudoJet> & incljets,
+                                   const JetDefinition & new_jet_def,
+                                   bool do_areas) const{
+  // create a new, internal, ClusterSequence from the jet constituents
+  // get the incljets directly from there
+  //
+  // If the jet has area support then we separate the ghosts from the
+  // "regular" particles so the incljets will also have area
+  // support. Note that we do this regardless of whether rho is zero
+  // or not.
+  //
+  // Note that to be able to separate the ghosts, one needs explicit
+  // ghosts!!
+  // ---------------------------------------------------------------
+  if (do_areas){
+    //vector<PseudoJet> all_constituents = jet.constituents();
+    vector<PseudoJet> regular_constituents, ghosts;  
+    SelectorIsPureGhost().sift(jet.constituents(), ghosts, regular_constituents);
+
+    // figure out the ghost area from the 1st ghost (if none, any value
+    // would probably do as the area will be 0 and subtraction will have
+    // no effect!)
+    double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
+    ClusterSequenceActiveAreaExplicitGhosts * csa
+      = new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents, 
+                                                    new_jet_def, 
+                                                    ghosts, ghost_area);
+
+    incljets = csa->inclusive_jets();
+
+    // allow the cs to be deleted when it's no longer used
+    // 
+    // Note that there is at least one constituent in the jet so there
+    // is in principle at least one incljet. But one may have used a
+    // nasty recombiner or jet def that left an empty set of incljets,
+    // so we'd rather play it safe (e.g. GridJetPlugin, with the
+    // constituents outside the range of the grid)
+    if (incljets.size())
+      csa->delete_self_when_unused();
+    else
+      delete csa;
+  } else {
+    ClusterSequence * cs = new ClusterSequence(jet.constituents(), new_jet_def);
+    incljets = cs->inclusive_jets();
+    // allow the cs to be deleted when it's no longer used (again, we
+    // add an extra safety check)
+    if (incljets.size())
+      cs->delete_self_when_unused();
+    else 
+      delete cs;
+  }
+}
+
+//----------------------------------------------------------------------
+// various checks and internal constructs
+//----------------------------------------------------------------------
+
+// fundamental info for CompositeJets
+//----------------------------------------------------------------------
+ 
+// get the pieces down to the fundamental pieces
+// 
+// Note that this just checks that there is an associated CS to the
+// fundamental pieces, not that it is still valid
+bool Recluster::_get_all_pieces(const PseudoJet &jet, vector<PseudoJet> &all_pieces) const{
+  if (jet.has_associated_cluster_sequence()){
+    all_pieces.push_back(jet);
+    return true;
+  }
+
+  if (jet.has_pieces()){
+    const vector<PseudoJet> pieces = jet.pieces();
+    for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
+      if (!_get_all_pieces(*it, all_pieces)) return false;
+    return true;
+  }
+
+  return false;
+}
+
+// construct the re-clustering jet definition using the recombiner
+// from whatever definition has been used to obtain the original jet
+//----------------------------------------------------------------------
+void Recluster::_acquire_recombiner_from_pieces(const vector<PseudoJet> &all_pieces, 
+                                                JetDefinition &new_jet_def) const{
+  // a quick safety check
+  assert(_acquire_recombiner);
+
+  // check that all the pieces have the same recombiner
+  //
+  // Note that if the jet has an associated cluster sequence that is no
+  // longer valid, an error will be thrown (needed since it could be the
+  // 1st check called after the enumeration of the pieces)
+  const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
+  for (unsigned int i=1; i<all_pieces.size(); i++){
+    if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref)){
+      throw Error("Recluster instance is configured to determine the recombination scheme (or recombiner) from the original jet, but different pieces of the jet were found to have non-equivalent recombiners.");
+    }
+  }
+
+  // get the recombiner from the original jet_def
+  new_jet_def.set_recombiner(jd_ref);
+}
+
+// area support
+//----------------------------------------------------------------------
+
+// check if the jet (or all its pieces) have explicit ghosts
+// (assuming the jet has area support).
+//
+// Note that if the jet has an associated cluster sequence that is no
+// longer valid, an error will be thrown (needed since it could be the
+// 1st check called after the enumeration of the pieces)
+bool Recluster::_check_explicit_ghosts(const vector<PseudoJet> &all_pieces) const{
+  for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
+    if (! it->validated_csab()->has_explicit_ghosts()) return false;
+  return true;
+}
+
+// C/A specific tests
+//----------------------------------------------------------------------
+
+// check if one can apply the simplification for C/A incljets
+//
+// This includes:
+//  - the incljet definition asks for C/A incljets
+//  - all the pieces share the same CS
+//  - that CS is C/A with the same recombiner as the incljet def
+//  - the re-clustering radius is not larger than any of the pairwise
+//    distance between the pieces
+//
+// Note that if the jet has an associated cluster sequence that is no
+// longer valid, an error will be thrown (needed since it could be the
+// 1st check called after the enumeration of the pieces)
+bool Recluster::_check_ca(const vector<PseudoJet> &all_pieces, 
+                          const JetDefinition &new_jet_def) const{
+  // check that optimisation is enabled
+  if (!_cambridge_optimisation_enabled) return false;
+
+  // check that we're reclustering with C/A
+  if (new_jet_def.jet_algorithm() != cambridge_algorithm) return false;
+
+  // check that the 1st of all the pieces (we're sure there is at
+  // least one) is coming from a C/A clustering. Then check that all
+  // the following pieces share the same ClusterSequence
+  const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
+  if (cs_ref->jet_def().jet_algorithm() != cambridge_algorithm) return false;
+  for (unsigned int i=1; i<all_pieces.size(); i++)
+    if (all_pieces[i].validated_cs() != cs_ref) return false;
+
+  // check that the 1st piece has the same recombiner as the one used
+  // for the incljet clustering
+  // Note that since they share the same CS, checking the 1st one is enough
+  if (!cs_ref->jet_def().has_same_recombiner(new_jet_def)) return false;
+
+  // we also have to make sure that the reclustering radius is not larger
+  // than any of the inter-piece distances
+  double Rnew2 = new_jet_def.R();
+  Rnew2 *= Rnew2;
+  for (unsigned int i=0; i<all_pieces.size()-1; i++){
+    for (unsigned int j=i+1; j<all_pieces.size(); j++){
+      if (all_pieces[i].squared_distance(all_pieces[j]) <  Rnew2) return false;
+    }
+  }
+
+  return true;
+}
+
+FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
Index: external/fastjet/tools/Recluster.hh
===================================================================
--- external/fastjet/tools/Recluster.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/tools/Recluster.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,262 @@
+#ifndef __FASTJET_TOOLS_RECLUSTER_HH__
+#define __FASTJET_TOOLS_RECLUSTER_HH__
+
+// $Id: Recluster.hh 3632 2014-08-15 10:42:53Z salam $
+//
+// Copyright (c) 2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+
+#include <fastjet/JetDefinition.hh>
+#include <fastjet/FunctionOfPseudoJet.hh>   // to derive Recluster from FOfPJ<PJ>
+#include <iostream>
+#include <string>
+
+// TODO:
+//
+//  - maintain Voronoi areas? Requires CSAB:has_voronoi_area()    {->UNASSIGNED}
+//  - make sure the description of the class is OK
+
+
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+//----------------------------------------------------------------------
+/// \class Recluster
+/// Recluster a jet's constituents with a new jet definition.
+///
+/// When Recluster is constructed from a JetDefinition, it is that
+/// definition that will be used to obtain the new jets. The user may
+/// then decide if the recombiner should be the one from that jet
+/// definition or if it should be acquired from the jet being
+/// processed (the default).
+///
+/// Alternatively, Recluster can be constructed from a jet algorithm
+/// and an optional radius. In that case the recombiner is
+/// systematically obtained fromn the jet being processed (unless you
+/// call set_acquire_recombiner(false)). If only the jet algorithm is
+/// specified, a default radius of max_allowable_R will be assumed if
+/// needed.
+///
+/// Recluster has two possible behaviours:
+///
+///  - if it is constructed with keep=keep_only_hardest the hardest
+///    inclusive jet is returned as a "standard" jet with an
+///    associated cluster sequence (unless there were no inclusive
+///    jets, in which case a zero jet is returned, with no associated
+///    cluster sequence)
+///
+///  - if it is constructed with keep=keep_all
+///    all the inclusive jets are joined into a composite jet
+///
+/// [Note that since the structure of the resulting PseudoJet depends
+/// on its usage, this class inherits from
+/// FunctionOfPseudoJet<PseudoJet> (including a description) rather
+/// than being a full-fledged Transformer]
+///
+class Recluster : public FunctionOfPseudoJet<PseudoJet> {
+public:
+  /// the various options for the output of Recluster
+  enum Keep{
+    /// keep only the hardest inclusive jet and return a "standard" jet with
+    /// an associated ClusterSequence [this will be the default]
+    keep_only_hardest,
+    /// keep all the inclusive jets. result() will join them into a composite
+    /// jet
+    keep_all
+  };
+
+  /// default constructor (uses an undefined JetDefinition, and so cannot
+  /// be used directly).
+  Recluster() : _new_jet_def(), _acquire_recombiner(true),
+                _keep(keep_only_hardest), _cambridge_optimisation_enabled(true){}
+
+  /// Constructs a Recluster object that reclusters a jet into a new jet
+  /// using a generic JetDefinition
+  ///
+  ///  \param new_jet_def   the jet definition applied to do the reclustering
+  ///  \param acquire_recombiner
+  ///                       when true, the reclustering will guess the
+  ///                       recombiner from the input jet instead of
+  ///                       the one in new_jet_def. An error is then
+  ///                       thrown if no consistent recombiner is found
+  ///  \param keep_in       Recluster::keep_only_hardest: the result is
+  ///                       the hardest inclusive jet after reclustering,
+  ///                       returned as a "standard" jet.
+  ///                       Recluster::keep_all: the result is a
+  ///                       composite jet with the inclusive jets as pieces.
+  Recluster(const JetDefinition & new_jet_def, 
+            bool acquire_recombiner_in = false, 
+            Keep keep_in = keep_only_hardest)
+    : _new_jet_def(new_jet_def), _acquire_recombiner(acquire_recombiner_in), 
+      _keep(keep_in), _cambridge_optimisation_enabled(true) {}
+
+  /// Constructs a Recluster object that reclusters a jet into a new jet
+  /// using a JetAlgorithm and its parameters
+  ///
+  ///  \param new_jet_alg    the jet algorithm applied to obtain the new clustering
+  ///  \param new_jet_radius the jet radius
+  ///  \param keep_in        Recluster::keep_only_hardest: the result is
+  ///                        the hardest inclusive jet after reclustering,
+  ///                        returned as a "standard" jet.
+  ///                        Recluster::keep_all: the result is a
+  ///                        composite jet with the inclusive jets as pieces.
+  /// 
+  /// This ctor will always acquire the recombiner from the jet being
+  /// reclustered (it will throw if none can be found).  If you wish
+  /// to use Recluster with an algorithm that requires an extra
+  /// parameter (like the genkt algorithm), please specify the jet
+  /// definition fully using the constructor above.
+  Recluster(JetAlgorithm snew_jet_alg, double new_jet_radius, Keep keep_in = keep_only_hardest);
+
+  /// constructor with just a jet algorithm, but no jet radius. If the
+  /// algorithm requires a jet radius, JetDefinition::max_allowable_R will be used. 
+  ///
+  Recluster(JetAlgorithm new_jet_alg, Keep keep_in = keep_only_hardest);
+
+  /// default dtor
+  virtual ~Recluster(){}
+
+  //----------------------------------------------------------------------
+  // tweaking the behaviour and corresponding enquiry functions
+
+  /// set whether the reclustering should attempt to acquire a
+  /// recombiner from the input jet
+  void set_acquire_recombiner(bool acquire) {_acquire_recombiner = acquire;}
+
+  /// returns true if this reclusterer is set to acquire the
+  /// recombiner from the input jet
+  bool acquire_recombiner() const{ return _acquire_recombiner;}
+
+
+  /// sets whether to try to optimise reclustering with
+  /// Cambridge/Aachen algorithms (by not reclustering if the
+  /// requested C/A reclustering can be obtained by using subjets of
+  /// an input C/A jet or one composed of multiple C/A pieces from the
+  /// same clustering sequence). By default this is enabled, and
+  /// _should_ always be correct; disable it to test this statement!
+  void set_cambridge_optimisation(bool enabled){ _cambridge_optimisation_enabled = enabled;}
+  /// sets whether to try to optimise reclustering with Cambridge/Aachen algorithms (US spelling!)
+  void set_cambridge_optimization(bool enabled){ _cambridge_optimisation_enabled = enabled;}
+
+  /// returns true if the reclusterer tries to optimise reclustering
+  /// with Cambridge/Aachen algorithms
+  bool cambridge_optimization(){return _cambridge_optimisation_enabled;}
+  bool cambridge_optimisation(){return _cambridge_optimisation_enabled;}
+
+  /// set the behaviour with regards to keeping all resulting jets or
+  /// just the hardest.
+  void set_keep(Keep keep_in) {_keep = keep_in;}
+
+  /// returns the current "keep" mode i.e. whether only the hardest
+  /// inclusive jet is returned or all of them (see the Keep enum above)
+  Keep keep() const{ return _keep;}
+
+
+  //----------------------------------------------------------------------
+  // retrieving info about the behaviour
+
+  /// class description
+  virtual std::string description() const;
+
+
+  //----------------------------------------------------------------------
+  // core action of ths class
+
+  /// runs the reclustering and sets kept and rejected to be the jets
+  /// of interest (with non-zero rho, they will have been
+  /// subtracted). Normally this will be accessed through the base
+  /// class's operator().
+  ///
+  /// \param jet    the jet that gets reclustered
+  /// \return the reclustered jet
+  virtual PseudoJet result(const PseudoJet & jet) const;
+
+  /// A lower-level method that does the actual work of reclustering
+  /// the input jet. The resulting jets are stored in output_jets.
+  /// The jet definition that has been used can be accessed from the
+  /// output_jets' ClusterSequence.
+  ///
+  /// \param input_jet       the (input) jet that one wants to recluster
+  /// \param output_jets     inclusive jets resulting from the new clustering
+  ///
+  /// Returns true if the C/A optimisation has been used (this means
+  /// that generate_output_jet then has to watch out for non-explicit-ghost
+  /// areas that might be leftover)
+  bool get_new_jets_and_def(const PseudoJet & input_jet, 
+                            std::vector<PseudoJet> & output_jets) const;
+
+  /// given a set of inclusive jets and a jet definition used, create the
+  /// resulting PseudoJet;
+  /// 
+  /// If ca_optimisation_used then special care will be taken in
+  /// deciding whether the final jet can legitimately have an area.
+  PseudoJet generate_output_jet(std::vector<PseudoJet> & incljets,
+                                bool ca_optimisation_used) const;
+
+
+private:
+  /// set the reclustered elements in the simple case of C/A+C/A
+  void _recluster_ca(const std::vector<PseudoJet> & all_pieces,
+                     std::vector<PseudoJet> & incljets,
+                     double Rfilt) const;
+
+  /// set the reclustered elements in the generic re-clustering case
+  void _recluster_generic(const PseudoJet & jet, 
+                          std::vector<PseudoJet> & incljets,
+                          const JetDefinition & new_jet_def,
+                          bool do_areas) const;
+  
+  // a series of checks
+  //--------------------------------------------------------------------
+  /// get the pieces down to the fundamental pieces
+  bool _get_all_pieces(const PseudoJet &jet, std::vector<PseudoJet> &all_pieces) const;
+
+  /// associate the proper recombiner taken from the underlying pieces
+  /// (an error is thrown if the pieces do no share a common
+  /// recombiner)
+  void _acquire_recombiner_from_pieces(const std::vector<PseudoJet> &all_pieces, 
+                                       JetDefinition &new_jet_def) const;
+
+  /// check if one can apply the simplified trick for C/A subjets
+  bool _check_ca(const std::vector<PseudoJet> &all_pieces, 
+                 const JetDefinition &new_jet_def) const;
+
+  /// check if the jet (or all its pieces) have explicit ghosts
+  /// (assuming the jet has area support
+  ///
+  /// Note that if the jet has an associated cluster sequence that is no
+  /// longer valid, an error will be thrown
+  bool _check_explicit_ghosts(const std::vector<PseudoJet> &all_pieces) const;
+
+  JetDefinition _new_jet_def;  ///< the jet definition to use to extract the jets
+  bool _acquire_recombiner;    ///< get the recombiner from the input
+                               ///< jet rather than from _new_jet_def
+  Keep _keep;                  ///< dictates which inclusive jets are kept and
+                               ///< which are returned (see Keep above)
+
+  bool _cambridge_optimisation_enabled; ///<enable the checks to
+                                        ///< perform optimisation when
+                                        ///< C/A reclustering is asked
+
+  static LimitedWarning   _explicit_ghost_warning;
+};
+
+FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
+
+#endif   // __FASTJET_TOOLS_RECLUSTER_HH__
Index: external/fastjet/tools/RestFrameNSubjettinessTagger.cc
===================================================================
--- external/fastjet/tools/RestFrameNSubjettinessTagger.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/RestFrameNSubjettinessTagger.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: RestFrameNSubjettinessTagger.cc 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/tools/RestFrameNSubjettinessTagger.hh>
Index: external/fastjet/tools/RestFrameNSubjettinessTagger.hh
===================================================================
--- external/fastjet/tools/RestFrameNSubjettinessTagger.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/RestFrameNSubjettinessTagger.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -2,8 +2,8 @@
 #define __FASTJET_RESTFRAMENSUBJETTINESS_TAGGER_HH__
 
-//STARTHEADER
-// $Id: RestFrameNSubjettinessTagger.hh 2689 2011-11-14 14:51:06Z soyez $
+//FJSTARTHEADER
+// $Id: RestFrameNSubjettinessTagger.hh 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -16,7 +16,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -28,5 +30,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/PseudoJet.hh>
Index: external/fastjet/tools/Subtractor.cc
===================================================================
--- external/fastjet/tools/Subtractor.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/Subtractor.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: Subtractor.cc 3594 2014-08-12 15:25:11Z soyez $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include "fastjet/tools/Subtractor.hh"
@@ -38,13 +40,127 @@
 
 
+//----------------------------------------------------------------------
+// ctor
 Subtractor::Subtractor(double rho) : _bge(0), _rho(rho) {
   assert(_rho>0.0);
+  set_defaults();
 }
 
+//----------------------------------------------------------------------
+void Subtractor::set_defaults(){
+  _use_rho_m = false; // likely to change in future releases!!
+  _safe_mass = false; // likely to change in future releases!!
+
+  _sel_known_vertex = Selector();
+  _sel_leading_vertex = Selector();
+}
+
+//----------------------------------------------------------------------
+// perform the subtraction of a given jet
 PseudoJet Subtractor::result(const PseudoJet & jet) const {
   if (!jet.has_area()){
     throw Error("Trying to subtract a jet without area support");
   }
-  
+
+  PseudoJet known_lv, known_pu;
+  PseudoJet unknown = jet;
+  if (_sel_known_vertex.worker()){
+    // separate the jet constituents in 3 groups:
+    //   unknown vertex
+    //   known vertex, leading vertex
+    //   known vertex, non-leading vertex (PU)
+    vector<PseudoJet> constits_unknown, constits_known;
+    _sel_known_vertex.sift(jet.constituents(), 
+			   constits_known,
+			   constits_unknown);
+    vector<PseudoJet> constits_known_lv, constits_known_pu;
+    _sel_leading_vertex.sift(constits_known,
+			     constits_known_lv,
+			     constits_known_pu);
+
+    // For the parts related to the known vertices (LV or PU), we just
+    // sum the 4-momenta. For the unknown part, we assign it the full
+    // jet area.
+    known_lv = (constits_known_lv.size()!=0) 
+      ? SelectorIdentity().sum(constits_known_lv) : 0.0*jet;
+    known_pu = (constits_known_pu.size()!=0) 
+      ? SelectorIdentity().sum(constits_known_pu) : 0.0*jet;
+    if (constits_unknown.size()==0){
+      // no need for any form of subtraction!
+      PseudoJet subtracted_jet = jet;
+      subtracted_jet.reset_momentum(known_lv);
+      return subtracted_jet;
+    }
+    unknown = jet; // that keeps all info including area
+    unknown.reset_momentum(SelectorIdentity().sum(constits_unknown));
+  } else {
+    known_lv = jet; // ensures correct rap-phi!
+    known_lv *= 0.0;
+    known_pu = known_lv;
+  }
+
+  // prepare for the subtraction and compute the 4-vector to be
+  // subtracted
+  PseudoJet subtracted_jet = jet;
+  PseudoJet to_subtract = known_pu + _amount_to_subtract(unknown);
+
+  // sanity check for the transverse momentum
+  if (to_subtract.pt2() < jet.pt2() ) { 
+    // this subtraction should retain the jet's structural
+    // information
+    subtracted_jet -= to_subtract;
+  } else { 
+    // this sets the jet's momentum while maintaining all of the jet's
+    // structural information
+    subtracted_jet.reset_momentum(known_lv);
+    return subtracted_jet;
+  }
+
+  // make sure that in the end the pt is at least the one known to
+  // come from the leading vertex
+  if (subtracted_jet.pt2() < known_lv.pt2()){
+    subtracted_jet.reset_momentum(known_lv);
+    return subtracted_jet;
+  }
+
+  // sanity check for the mass (if needed)
+  if ((_safe_mass) && (subtracted_jet.m2() < known_lv.m2())){
+    // in this case, we keep pt and phi as obtained from the
+    // subtraction above and take rap and m from the part that comes
+    // from the leading vertex (or the original jet if nothing comes
+    // from the leading vertex)
+    subtracted_jet.reset_momentum(PtYPhiM(subtracted_jet.pt(),
+					  known_lv.rap(),
+					  subtracted_jet.phi(),
+					  known_lv.m()));
+  }
+
+  return subtracted_jet;
+}
+
+//----------------------------------------------------------------------
+std::string Subtractor::description() const{
+  if (_bge != 0) {
+    string desc = "Subtractor that uses the following background estimator to determine rho: "+_bge->description();
+    if (use_rho_m()) desc += "; including the rho_m correction";
+    if (safe_mass()) desc += "; including mass safety tests";
+    if (_sel_known_vertex.worker()){
+      desc += "; using known vertex selection: "+_sel_known_vertex.description()+" and leading vertex selection: "+_sel_leading_vertex.description();
+    }
+    return desc;
+  } else if (_rho != _invalid_rho) {
+    ostringstream ostr;
+    ostr << "Subtractor that uses a fixed value of rho = " << _rho;
+    return ostr.str();
+  } else {
+    return "Uninitialised subtractor";
+  }
+}
+
+//----------------------------------------------------------------------
+// compute the 4-vector that should be subtracted from the given
+// jet
+PseudoJet Subtractor::_amount_to_subtract(const PseudoJet &jet) const{
+  // the "transverse momentum" part
   double rho;
   if (_bge != 0) {
@@ -56,31 +172,22 @@
   }
 
-  PseudoJet subtracted_jet = jet;
-  PseudoJet area4vect = jet.area_4vector();
-  // sanity check
-  if (rho*area4vect.perp() < jet.perp() ) { 
-    // this subtraction should retain the jet's structural
-    // information
-    subtracted_jet -= rho*area4vect;
-  } else { 
-    // this sets the jet's momentum to zero while
-    // maintaining all of the jet's structural information
-    subtracted_jet *= 0;
+  PseudoJet area = jet.area_4vector();
+  PseudoJet to_subtract = rho*area;
+
+  double const rho_m_warning_threshold = 1e-5;
+
+  // add an optional contribution from the unknown particles masses
+  if (_use_rho_m){
+    assert(_bge != 0); // test done in "set_use_rho_m()"
+    to_subtract += _bge->rho_m(jet) * PseudoJet(0.0, 0.0, area.pz(), area.E());
+  } else if (_bge && 
+             _bge->has_rho_m() && 
+             _bge->rho_m(jet) > rho_m_warning_threshold * rho) {
+    _unused_rho_m_warning.warn("Background estimator indicates non-zero rho_m, but use_rho_m()==false in subtractor; consider calling set_use_rho_m(true) to include the rho_m information");
   }
-  return subtracted_jet;
+
+  return to_subtract;
 }
 
-//----------------------------------------------------------------------
-std::string Subtractor::description() const{
-  if (_bge != 0) {
-    return "Subtractor that uses the following background estimator to determine rho: "+_bge->description();
-  } else if (_rho != _invalid_rho) {
-    ostringstream ostr;
-    ostr << "Subtractor that uses a fixed value of rho = " << _rho;
-    return ostr.str();
-  } else {
-    return "Uninitialised subtractor";
-  }
-}
 
 FASTJET_END_NAMESPACE
Index: external/fastjet/tools/Subtractor.hh
===================================================================
--- external/fastjet/tools/Subtractor.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/Subtractor.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id: Subtractor.hh 2577 2011-09-13 15:11:38Z salam $
-//
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//FJSTARTHEADER
+// $Id: Subtractor.hh 3584 2014-08-12 12:40:10Z soyez $
+//
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #ifndef __FASTJET_TOOLS_SUBTRACTOR_HH__
@@ -61,5 +63,5 @@
   /// define a subtractor based on a BackgroundEstimator
   Subtractor(BackgroundEstimatorBase * bge) : 
-    _bge(bge), _rho(-1.0) {}
+    _bge(bge), _rho(-1.0) { set_defaults(); }
 
   /// define a subtractor that uses a fixed value of rho, the background
@@ -68,8 +70,91 @@
 
   /// default constructor
-  Subtractor() : _bge(0), _rho(_invalid_rho) {}
+  Subtractor() : _bge(0), _rho(_invalid_rho) { set_defaults(); }
 
   /// default dtor
   virtual ~Subtractor(){};
+
+  /// @name configuring the behaviour
+  //\{
+  //----------------------------------------------------------------
+
+  /// reset all parameters to default values
+  ///
+  /// Note: by default, the rho_m term is not included and the safety
+  /// test for the mass is not done. This is mostly for backwards
+  /// compatibility with FastJet 3.0 and is highly likely to change in
+  /// a future release of FastJet
+  void set_defaults();
+
+  /// when 'use_rho_m' is true, include in the subtraction the
+  /// correction from rho_m, the purely longitudinal,
+  /// particle-mass-induced component of the background density per
+  /// unit area
+  ///
+  /// Note: this will be switched off by default (for backwards
+  /// compatibility with FastJet 3.0) but is highly likely to change
+  /// in a future release of FastJet
+  void set_use_rho_m(bool use_rho_m_in = true){
+    if (_bge == 0) {
+      throw Error("Subtractor: rho_m support works only for Subtractors constructed with background estimator");
+    }
+    _use_rho_m=use_rho_m_in;
+  }
+  
+  /// returns whether or not the rho_m component is used
+  bool use_rho_m() const{ return _use_rho_m;}
+
+  /// when 'safe_mass' is true, ensure that the mass of the subtracted
+  /// 4-vector remain positive
+  ///
+  /// when true, if the subtracted mass is negative, we return a
+  /// 4-vector with 0 mass, pt and phi from the subtracted 4-vector
+  /// and the rapidity of the original, unsubtracted jet.
+  ///
+  /// Note: this will be switched off by default (for backwards
+  /// compatibility with FastJet 3.0) but is highly likely to change
+  /// in a future release of FastJet
+  void set_safe_mass(bool safe_mass_in=true){ _safe_mass=safe_mass_in;}
+
+  /// returns whether or not safety tests on the mass are included
+  bool safe_mass() const{ return _safe_mass;}
+
+  /// This is mostly intended for cherge-hadron-subtracted type of
+  /// events where we wich to use vertex information to improve the
+  /// subtraction.
+  ///
+  /// Given the following parameters:
+  ///   \param sel_known_vertex    selects the particles with a
+  ///                              known vertex origin
+  ///   \param sel_leading_vertex  amongst the particles with a
+  ///                              known vertex origin, select those
+  ///                              coming from the leading vertex
+  /// Momentum identified as coming from the leading vertex will be
+  /// kept, momentum identified as coming from a non-leading vertex
+  /// will be eliminated and a regular area-median subtraction will be
+  /// applied on the 4-vector sum of the particles with unknown vertex
+  /// origin.
+  ///
+  /// When this is set, we shall ensure that the pt of the subtracted
+  /// 4-vector is at least the pt of the particles that are known to
+  /// come from the leading vertex (if it fails, subtraction returns
+  /// the component that is known to come from the leading vertex ---
+  /// or, the original unsubtracted jet if it contains no particles
+  /// from the leading vertex).  Furthermore, when safe_mass() is on, we
+  /// also impose a similar constraint on the mass of the subtracted
+  /// 4-vector (if the test fails, the longitudinal part of the
+  /// subtracted 4-vector is taken from the component that is known to
+  /// come from the leading vertex).
+  void set_known_selectors(const Selector &sel_known_vertex,
+			   const Selector &sel_leading_vertex){
+    _sel_known_vertex   = sel_known_vertex;
+    _sel_leading_vertex = sel_leading_vertex;
+  }
+
+  //\}
+
+  /// @name description and action
+  //\{
+  //----------------------------------------------------------------
 
   /// returns a jet that's subtracted
@@ -82,5 +167,9 @@
   virtual std::string description() const;
 
+  //\}
 protected:
+  /// compute the 4-vector that should be subtracted from the given
+  /// jet
+  PseudoJet _amount_to_subtract(const PseudoJet &jet) const;
 
   /// the tool used to estimate the background
@@ -89,4 +178,14 @@
   /// the fixed value of rho to use if the user has selected that option
   double _rho;
+
+  // configuration parameters/flags
+  bool _use_rho_m;   ///< include the rho_m correction
+  bool _safe_mass;   ///< ensures that the subtracted mass is +ve
+
+  Selector _sel_known_vertex;   ///< selects the particles with a
+				///< known vertex origin
+  Selector _sel_leading_vertex; ///< amongst the particles with a
+				///< known vertex origin, select those
+				///< coming from the leading vertex
 
   /// a value of rho that is used as a default to label that the stored
@@ -98,4 +197,6 @@
   // that's not allowed in an include file.
   static const double _invalid_rho;
+
+  mutable LimitedWarning _unused_rho_m_warning;
 };
 
Index: external/fastjet/tools/TopTaggerBase.cc
===================================================================
--- external/fastjet/tools/TopTaggerBase.cc	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/TopTaggerBase.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id$
+//FJSTARTHEADER
+// $Id: TopTaggerBase.cc 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/tools/TopTaggerBase.hh>
Index: external/fastjet/tools/TopTaggerBase.hh
===================================================================
--- external/fastjet/tools/TopTaggerBase.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/TopTaggerBase.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -2,8 +2,8 @@
 #define __FASTJET_TOP_TAGGER_BASE_HH__
 
-//STARTHEADER
-// $Id: TopTaggerBase.hh 2689 2011-11-14 14:51:06Z soyez $
+//FJSTARTHEADER
+// $Id: TopTaggerBase.hh 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -16,7 +16,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -28,5 +30,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #include <fastjet/internal/base.hh>
Index: external/fastjet/tools/Transformer.hh
===================================================================
--- external/fastjet/tools/Transformer.hh	(revision d69dfe4aebfe1009f10ecbbabd92ce813ba93f3e)
+++ external/fastjet/tools/Transformer.hh	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -1,6 +1,6 @@
-//STARTHEADER
-// $Id: Transformer.hh 2577 2011-09-13 15:11:38Z salam $
+//FJSTARTHEADER
+// $Id: Transformer.hh 3433 2014-07-23 08:17:03Z salam $
 //
-// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
@@ -13,7 +13,9 @@
 //
 //  The algorithms that underlie FastJet have required considerable
-//  development and are described in hep-ph/0512210. If you use
+//  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
-//  include a citation to the FastJet paper.
+//  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,
@@ -25,5 +27,5 @@
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
-//ENDHEADER
+//FJENDHEADER
 
 #ifndef __FASTJET_TRANSFORMER_HH__
