Index: external/fastjet/contribs/RecursiveTools/AUTHORS
===================================================================
--- external/fastjet/contribs/RecursiveTools/AUTHORS	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/AUTHORS	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -4,4 +4,6 @@
    Gregory Soyez <soyez@fastjet.fr>
    Jesse Thaler <jthaler@jthaler.net>
+   Kevin Zhou <knzhou@mit.edu>
+   Frederic Dreyer <fdreyer@mit.edu>
 
 The physics is based on:
@@ -16,2 +18,14 @@
    Andrew J. Larkoski, Simone Marzani, Gregory Soyez, and Jesse Thaler.
    JHEP 1405:146 (2014), arXiv:1402.2657
+
+   [IteratedSoftDrop]
+   Casimir Meets Poisson: Improved Quark/Gluon Discrimination with Counting Observables.
+   Christopher Frye, Andrew J. Larkoski, Jesse Thaler, Kevin Zhou.
+   JHEP 1709:083 (2017), arXiv:1704.06266
+
+   [RecursiveSoftDrop]
+   [BottomUpSoftDrop]
+   Recursive Soft Drop.
+   Frederic A. Dreyer, Lina Necib, Gregory Soyez, and Jesse Thaler
+   arXiv:1804.03657
+   
Index: external/fastjet/contribs/RecursiveTools/BottomUpSoftDrop.cc
===================================================================
--- external/fastjet/contribs/RecursiveTools/BottomUpSoftDrop.cc	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
+++ external/fastjet/contribs/RecursiveTools/BottomUpSoftDrop.cc	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -0,0 +1,324 @@
+// $Id: BottomUpSoftDrop.cc 1064 2017-09-08 09:19:57Z gsoyez $
+//
+// Copyright (c) 2017-, Gavin P. Salam, Gregory Soyez, Jesse Thaler,
+// Kevin Zhou, Frederic Dreyer
+//
+//----------------------------------------------------------------------
+// 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 "BottomUpSoftDrop.hh"
+#include <cassert>
+#include <algorithm>
+#include <sstream>
+#include <typeinfo>
+#include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
+#include "fastjet/Selector.hh"
+#include "fastjet/config.h"
+
+
+using namespace std;
+
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+//----------------------------------------------------------------------
+// BottomUpSoftDrop class
+//----------------------------------------------------------------------
+  
+// action on a single jet
+PseudoJet BottomUpSoftDrop::result(const PseudoJet &jet) const{
+  // soft drop can only be applied to jets that have constituents
+  if (!jet.has_constituents()){
+    throw Error("BottomUpSoftDrop: trying to apply the Soft Drop transformer to a jet that has no constituents");
+  }
+
+  // if the jet has area support and there are explicit ghosts, we can
+  // transfer that support to the internal re-clustering
+  //
+  // Note that this is just meant to maintain the information since
+  // all the jes will have a 0 area
+  bool do_areas = jet.has_area() && _check_explicit_ghosts(jet);
+
+  // build the soft drop plugin
+  BottomUpSoftDropPlugin * softdrop_plugin;
+
+  // for some constructors, we get the recombiner from the 
+  // input jet -- some acrobatics are needed
+  if (_get_recombiner_from_jet) {
+    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)){
+#if FASTJET_VERSION_NUMBER >= 30100
+      // Note that this is better than the option directly passing the
+      // recombiner (for cases where th ejet def own its recombiner)
+      // but it's only available in FJ>=3.1
+      jet_def.set_recombiner(jet_def_for_recombiner);
+#else
+      jet_def.set_recombiner(jet_def_for_recombiner.recombiner());
+#endif
+    }
+    softdrop_plugin = new BottomUpSoftDropPlugin(jet_def, _beta, _symmetry_cut, _R0);
+  } else {
+    softdrop_plugin = new BottomUpSoftDropPlugin(_jet_def, _beta, _symmetry_cut, _R0);
+  }
+
+  // now recluster the constituents of the jet with that plugin
+  JetDefinition internal_jet_def(softdrop_plugin);
+  // flag the plugin for automatic deletion _before_ we make
+  // copies (so that as long as the copies are also present
+  // it doesn't get deleted).
+  internal_jet_def.delete_plugin_when_unused();
+
+  ClusterSequence * cs;
+  if (do_areas){
+    vector<PseudoJet> particles, ghosts;
+    SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles);
+    // determine the ghost area from the 1st ghost (if none, any value
+    // will do, as the area will be 0 and subtraction will have
+    // no effect!)
+    double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
+    cs = new ClusterSequenceActiveAreaExplicitGhosts(particles, internal_jet_def, 
+  						     ghosts, ghost_area);
+  } else {
+    cs = new ClusterSequence(jet.constituents(), internal_jet_def);
+  }
+  
+  PseudoJet result_local = SelectorNHardest(1)(cs->inclusive_jets())[0];
+  BottomUpSoftDropStructure * s = new BottomUpSoftDropStructure(result_local);
+  s->_beta = _beta;
+  s->_symmetry_cut = _symmetry_cut;
+  s->_R0   = _R0;
+  result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
+  
+  // make sure things remain persistent -- i.e. tell the jet definition
+  // and the cluster sequence that it is their responsibility to clean 
+  // up memory once the "result" reaches the end of its life in the user's
+  // code. (The CS deletes itself when the result goes out of scope and
+  // that also triggers deletion of the plugin)
+  cs->delete_self_when_unused();
+
+  return result_local;  
+}
+
+// global grooming on a full event
+// note: does not support jet areas
+vector<PseudoJet> BottomUpSoftDrop::global_grooming(const vector<PseudoJet> & event) const {
+  // start by reclustering the event into one very large jet
+  ClusterSequence cs(event, _jet_def);
+  std::vector<PseudoJet> global_jet = SelectorNHardest(1)(cs.inclusive_jets());
+  // if the event is empty, do nothing
+  if (global_jet.size() == 0) return vector<PseudoJet>();
+  
+  // apply the groomer to the large jet
+  PseudoJet result = this->result(global_jet[0]);
+  return result.constituents();
+}
+  
+// check if the jet has explicit_ghosts (knowing that there is an
+// area support)
+bool BottomUpSoftDrop::_check_explicit_ghosts(const PseudoJet &jet) const {
+  // if the jet comes from a Clustering check explicit ghosts in that
+  // clustering
+  if (jet.has_associated_cluster_sequence())
+    return jet.validated_csab()->has_explicit_ghosts();
+
+  // if the jet has pieces, recurse in the pieces
+  if (jet.has_pieces()){
+    vector<PseudoJet> pieces = jet.pieces();
+    for (unsigned int i=0;i<pieces.size(); i++)
+      if (!_check_explicit_ghosts(pieces[i])) return false;
+    // never returned false, so we're OK.
+    return true;
+  }
+
+  // return false for any other (unknown) structure
+  return false;
+}
+
+// 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 BottomUpSoftDrop::_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 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 true;
+  }
+  // return false for any other (unknown) structure
+  return false;
+}
+
+
+// description
+string BottomUpSoftDrop::description() const{
+  ostringstream oss;
+  oss << "BottomUpSoftDrop with jet_definition = (" << _jet_def.description() << ")"
+      << ", symmetry_cut = " << _symmetry_cut	<< ", beta = " << _beta << ", R0 = " << _R0;
+  return oss.str();
+}
+  
+//----------------------------------------------------------------------
+// BottomUpSoftDropRecombiner class
+//----------------------------------------------------------------------
+
+// perform a recombination taking into account the Soft Drop
+// conditions
+void BottomUpSoftDropRecombiner::recombine(const PseudoJet &pa, 
+                                           const PseudoJet &pb,
+                                           PseudoJet &pab) const {
+  PseudoJet p;
+  _recombiner->recombine(pa, pb, p);
+
+  // calculate symmetry parameter
+  double symmetry_cut_fn =_symmetry_cut * pow(pa.squared_distance(pb)/_R0sqr, 0.5*_beta);
+  double pta = pa.pt();
+  double ptb = pb.pt();
+  // make sure denominator is non-zero
+  double sym = pta + ptb;
+  if (sym == 0){
+    pab = p;
+    return;
+  }
+  sym = min(pta, ptb) / sym;
+    
+  // if the 2 particles pass the soft drop condition, do the
+  // recombination
+  if (sym > symmetry_cut_fn){
+    pab=p;
+    return;
+  }
+
+  // if the soft drop condition is not passed, return the particle
+  // with largest pt
+  if (pta < ptb) {
+    pab = pb;
+    _rejected.push_back(pa.cluster_hist_index());
+  } else {
+    pab = pa;
+    _rejected.push_back(pb.cluster_hist_index());
+  }
+}
+
+//----------------------------------------------------------------------
+// BottomUpSoftDropPlugin class
+//----------------------------------------------------------------------
+
+// the actual clustering work for the plugin
+void BottomUpSoftDropPlugin::run_clustering(ClusterSequence &input_cs) const {
+  // declare a soft drop recombiner
+  BottomUpSoftDropRecombiner softdrop_recombiner(_beta, _symmetry_cut, _R0, _jet_def.recombiner());
+  JetDefinition jet_def = _jet_def;
+  jet_def.set_recombiner(&softdrop_recombiner);
+
+  // cluster the particles using that recombiner
+  ClusterSequence internal_cs(input_cs.jets(), jet_def);
+  const vector<ClusterSequence::history_element> & internal_hist = internal_cs.history();
+
+  // transfer the list of "childless" elements into a bool vector
+  vector<bool> kept(internal_hist.size(), true);
+  const vector<unsigned int> &sd_rej = softdrop_recombiner.rejected();
+  for (unsigned int i=0;i<sd_rej.size(); i++) kept[sd_rej[i]]=false;
+
+  // browse the history, keeping only the elements that have not been
+  // vetoed.
+  //
+  // In the process we build a map for the history indices
+  vector<unsigned int> internal2input(internal_hist.size());
+  for (unsigned int i=0; i<input_cs.jets().size(); i++)
+    internal2input[i] = i;
+
+  for (unsigned int i=input_cs.jets().size(); i<internal_hist.size(); i++){
+    const ClusterSequence::history_element &he = internal_hist[i];
+
+    // deal with recombinations with the beam
+    if (he.parent2 == ClusterSequence::BeamJet){
+      int internal_jetp_index = internal_hist[he.parent1].jetp_index;
+      int internal_hist_index = internal_cs.jets()[internal_jetp_index].cluster_hist_index();
+
+      int input_jetp_index = input_cs.history()[internal2input[internal_hist_index]].jetp_index;
+
+      // cout << "Beam recomb for internal " << internal_hist_index
+      //            << " (input jet index=" << input_jetp_index << endl;
+
+      input_cs.plugin_record_iB_recombination(input_jetp_index, he.dij);
+      continue;
+    }
+
+    // now, deal with two-body recombinations
+    if (!kept[he.parent1]){ // 1 is rejected, we keep only 2
+      internal2input[i]=internal2input[he.parent2];
+      // cout << "rejecting internal " << he.parent1
+      //            << ", mapping internal " << i 
+      //            << " to internal " << he.parent2
+      //            << " i.e. " << internal2input[i] << endl;
+    } else if (!kept[he.parent2]){ // 2 is rejected, we keep only 1
+      internal2input[i]=internal2input[he.parent1];
+      // cout << "rejecting internal " << he.parent2 
+      //            << ", mapping internal " << i 
+      //            << " to internal " << he.parent1
+      //            << " i.e. " << internal2input[i] << endl;
+    } else { // do the recombination
+      int new_index;
+      input_cs.plugin_record_ij_recombination(input_cs.history()[internal2input[he.parent1]].jetp_index,
+					      input_cs.history()[internal2input[he.parent2]].jetp_index,
+					      he.dij, internal_cs.jets()[he.jetp_index], new_index);
+      internal2input[i]=input_cs.jets()[new_index].cluster_hist_index();
+      // cout << "merging " << internal2input[he.parent1] << " (int: " << he.parent1 << ")"
+      //      << " and "    << internal2input[he.parent2] << " (int: " << he.parent2 << ")"
+      //      << " into "   << internal2input[i] << " (int: " << i << ")" << endl;
+    }
+  }
+}
+
+// description of the plugin
+string BottomUpSoftDropPlugin::description() const {
+  ostringstream oss;
+  oss << "BottomUpSoftDropPlugin with jet_definition = (" << _jet_def.description()
+      <<"), symmetry_cut = " << _symmetry_cut << ", beta = "
+      << _beta << ", R0 = " << _R0;
+  return oss.str();
+}
+
+
+}
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/RecursiveTools/BottomUpSoftDrop.hh
===================================================================
--- external/fastjet/contribs/RecursiveTools/BottomUpSoftDrop.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
+++ external/fastjet/contribs/RecursiveTools/BottomUpSoftDrop.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -0,0 +1,303 @@
+// $Id: BottomUpSoftDrop.hh 1085 2017-10-11 02:16:59Z jthaler $
+//
+// Copyright (c) 2017-, Gavin P. Salam, Gregory Soyez, Jesse Thaler,
+// Kevin Zhou, Frederic Dreyer
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+
+#ifndef __BOTTOMUPSOFTDROP_HH__
+#define __BOTTOMUPSOFTDROP_HH__
+
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/WrappedStructure.hh"
+#include "fastjet/tools/Transformer.hh"
+
+#include <iostream>
+#include <string>
+
+// TODO
+//
+//  - missing class description
+//
+//  - check what to do when pta=ptb=0
+//    for the moment, we recombine both for multiple reasons
+//     . this avois breakingteh symemtry between pa and pb
+//     . it would be groomed in later steps anyway
+//    Note that this is slightly inconsistent with our use of
+//    > (instead of >=) in the cdt 
+
+FASTJET_BEGIN_NAMESPACE
+
+namespace contrib{
+  
+// fwd declarations
+class BottomUpSoftDrop;
+class BottomUpSoftDropStructure;
+class BottomUpSoftDropRecombiner;
+class BottomUpSoftDropPlugin;
+
+//----------------------------------------------------------------------
+/// \class BottomUpSoftDrop
+/// Implementation of the BottomUpSoftDrop transformer
+///
+/// Bottom-Up Soft drop grooms a jet by applying a modified
+/// recombination scheme, where particles are recombined only if they
+/// pass the Soft Drop condition
+///
+/// \f[
+///   z < z_{\rm cut} (\theta/R0)^\beta
+/// \f]
+///
+/// the groomed jet contains the particles remaining after this
+/// pair-wise recombination
+///
+/// Note:
+///  - one can use BottomUpSoftDrop on a full event with the
+///    global_grooming(event) method.
+///  - if two recombined particles a and b have momentum pta=ptb=0,
+///    we recombine both.
+/// 
+
+  
+class BottomUpSoftDrop : public Transformer {
+public:
+  /// minimal constructor, which the jet algorithm to CA, sets the radius
+  /// to JetDefinition::max_allowable_R (practically equivalent to
+  /// infinity) and also tries to use a recombiner based on the one in
+  /// the jet definition of the particular jet being Soft Dropped.
+  ///
+  ///  \param beta           the value for beta
+  ///  \param symmetry_cut   the value for symmetry_cut
+  ///  \param R0             the value for R0
+  BottomUpSoftDrop(double beta, double symmetry_cut, double R0 = 1.0) 
+    : _jet_def(cambridge_algorithm, JetDefinition::max_allowable_R),
+      _beta(beta),_symmetry_cut(symmetry_cut), _R0(R0),
+      _get_recombiner_from_jet(true) {}
+
+  /// alternative constructor which takes a specified jet algorithm
+  ///
+  ///  \param jet_alg        the jet algorithm for the internal clustering (uses R=infty)
+  ///  \param symmetry_cut   the value of symmetry_cut
+  ///  \param beta           the value for beta
+  ///  \param R0             the value for R0
+  BottomUpSoftDrop(const JetAlgorithm jet_alg, double beta, double symmetry_cut,
+		   double R0 = 1.0) 
+    : _jet_def(jet_alg, JetDefinition::max_allowable_R),
+      _beta(beta), _symmetry_cut(symmetry_cut), _R0(R0),
+      _get_recombiner_from_jet(true) {}
+
+
+  /// alternative ctor in which the full reclustering jet definition can
+  /// be specified.
+  ///
+  ///  \param jet_def        the jet definition for the internal clustering
+  ///  \param symmetry_cut   the value of symmetry_cut
+  ///  \param beta           the value for beta
+  ///  \param R0             the value for R0
+  BottomUpSoftDrop(const JetDefinition &jet_def, double beta, double symmetry_cut,
+		   double R0 = 1.0)
+    : _jet_def(jet_def), _beta(beta), _symmetry_cut(symmetry_cut), _R0(R0),
+      _get_recombiner_from_jet(false) {}
+
+  /// action on a single jet
+  virtual PseudoJet result(const PseudoJet &jet) const;
+
+  /// global grooming on a full event
+  /// note: does not support jet areas
+  virtual std::vector<PseudoJet> global_grooming(const std::vector<PseudoJet> & event) const;
+  
+  /// description
+  virtual std::string description() const;
+  
+  // the type of the associated structure
+  typedef BottomUpSoftDropStructure StructureType;
+
+private:
+  /// check if the jet has explicit_ghosts (knowing that there is an
+  /// area support)
+  bool _check_explicit_ghosts(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 _beta;           ///< the value of beta
+  double _symmetry_cut;   ///< the value of symmetry_cut
+  double _R0;             ///< the value of R0
+  bool   _get_recombiner_from_jet; ///< true for minimal constructor,
+                                   ///< causes recombiner to be set equal 
+                                   ///< to that already used in the jet 
+                                   ///< (if it can be deduced)
+};
+
+//----------------------------------------------------------------------
+/// The structure associated with a PseudoJet thas has gone through a
+/// bottom/up SoftDrop transformer
+class BottomUpSoftDropStructure : public WrappedStructure{
+public:
+  /// default ctor
+  ///  \param result_jet  the jet for which we have to keep the structure
+  BottomUpSoftDropStructure(const PseudoJet & result_jet)
+    : WrappedStructure(result_jet.structure_shared_ptr()){}
+
+  /// description
+  virtual std::string description() const{
+    return "Bottom/Up Soft Dropped PseudoJet";
+  }
+
+  /// return the constituents that have been rejected
+  std::vector<PseudoJet> rejected() const{ 
+    return validated_cs()->childless_pseudojets();
+  }
+
+  /// return the other jets that may have been found along with the
+  /// result of the bottom/up Soft Drop
+  /// The resulting vector is sorted in pt
+  std::vector<PseudoJet> extra_jets() const {
+    return sorted_by_pt((!SelectorNHardest(1))(validated_cs()->inclusive_jets()));
+  }
+
+  /// return the value of beta that was used for this specific Soft Drop.
+  double beta() const {return _beta;}
+
+  /// return the value of symmetry_cut that was used for this specific Soft Drop.
+  double symmetry_cut() const {return _symmetry_cut;}
+
+  /// return the value of R0 that was used for this specific Soft Drop.
+  double R0() const {return _R0;}
+
+protected:
+  friend class BottomUpSoftDrop; ///< to allow setting the internal information
+
+private:
+  double _beta, _symmetry_cut, _R0;
+};
+
+//----------------------------------------------------------------------
+/// Class for Soft Drop recombination
+/// recombines the objects that are not vetoed by Bottom-Up SoftDrop
+///
+/// This recombiner only recombines, using the provided 'recombiner',
+/// objects (i and j) that pass the following SoftDrop criterion:
+///
+///   min(pti, ptj) > zcut (pti+ptj) (theta_ij/R0)^beta
+///
+/// If the criterion fail, the hardest of i and j is kept and the
+/// softest is rejected.
+///
+/// Note that this in not meant for standalone use [in particular
+/// because it could lead to memory issues due to the rejected indices
+/// stored internally].
+///
+/// This class is a direct adaptation of PruningRecombiner in Fastjet tools
+class BottomUpSoftDropRecombiner : public JetDefinition::Recombiner {
+public:
+  /// ctor
+  ///  \param symmetry_cut   value of cut on symmetry measure
+  ///  \param beta   avalue of beta parameter
+  ///  \param recomb pointer to a recombiner to use to cluster pairs
+  BottomUpSoftDropRecombiner(double beta, double symmetry_cut, double R0,
+                             const JetDefinition::Recombiner *recombiner)
+    : _beta(beta), _symmetry_cut(symmetry_cut), _R0sqr(R0*R0),
+      _recombiner(recombiner) {}
+
+  /// perform a recombination taking into account the Soft Drop
+  /// conditions
+  virtual void recombine(const PseudoJet &pa, 
+			 const PseudoJet &pb,
+			 PseudoJet &pab) const;
+
+  /// returns the description of the recombiner
+  virtual std::string description() const {
+    std::ostringstream oss;
+    oss << "SoftDrop recombiner with symmetry_cut = " << _symmetry_cut
+	<< ", beta = " << _beta
+	<< ", and underlying recombiner = " << _recombiner->description();
+    return oss.str();    
+  }
+
+  /// return the history indices that have been soft dropped away
+  const std::vector<unsigned int> & rejected() const{ return _rejected;}
+
+  /// clears the list of rejected indices
+  ///
+  /// If one decides to use this recombiner standalone, one has to
+  /// call this after each clustering in order for the rejected() vector
+  /// to remain sensible and not grow to infinite size.
+  void clear_rejected(){ _rejected.clear();}
+
+private:
+  double _beta;          ///< beta parameter
+  double _symmetry_cut;  ///< value of symmetry_cut
+  double _R0sqr;         ///< normalisation of the angular distance
+  const JetDefinition::Recombiner *_recombiner; ///< the underlying recombiner to use
+  mutable std::vector<unsigned int> _rejected;  ///< list of rejected history indices
+};
+
+//----------------------------------------------------------------------
+/// \class BottomUpSoftDropPlugin
+/// Class for a bottom/up Soft Drop algorithm, based on the Pruner plugin
+///
+/// This is an internal plugin that clusters the particles using the
+/// BottomUpRecombiner. 
+///
+/// See BottomUpRecombiner for a description of what bottom-up
+/// SoftDrop does.
+///
+/// Note that this is an internal class used by the BottomUpSoftDrop
+/// transformer and it is not meant to be used as a standalone
+/// clustering tool.
+class BottomUpSoftDropPlugin : public JetDefinition::Plugin {
+public:
+  /// ctor
+  ///  \param jet_def the jet definition to be used for the 
+  ///                 internal clustering
+  ///  \param symmetry_cut    value of cut on symmetry measure
+  ///  \param beta    value of beta parameter
+  BottomUpSoftDropPlugin(const JetDefinition &jet_def, double beta, double symmetry_cut,
+                         double R0 = 1.0)
+    : _jet_def(jet_def), _beta(beta), _symmetry_cut(symmetry_cut), _R0(R0) {}
+
+  /// the actual clustering work for the plugin
+  virtual void run_clustering(ClusterSequence &input_cs) const;
+
+  /// description of the plugin
+  virtual std::string description() const;
+
+  /// returns the radius
+  virtual double R() const {return _jet_def.R();}
+
+private:
+  JetDefinition _jet_def; ///< the internal jet definition
+  double _beta;           ///< beta parameter
+  double _symmetry_cut;   ///< value of symmetry_cut 
+  double _R0;             ///< normalisation of the angular distance
+};
+
+}
+
+FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
+#endif // __BOTTOMUPSOFTDROP_HH__
Index: external/fastjet/contribs/RecursiveTools/COPYING
===================================================================
--- external/fastjet/contribs/RecursiveTools/COPYING	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/COPYING	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -14,4 +14,5 @@
   arXiv:1307.0007
   arXiv:1402.2657
+  arXiv:1704.06266
 
 ======================================================================
Index: external/fastjet/contribs/RecursiveTools/ChangeLog
===================================================================
--- external/fastjet/contribs/RecursiveTools/ChangeLog	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/ChangeLog	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,2 +1,396 @@
+2018-05-29  Jesse Thaler  <jthaler@jthaler.net>
+
+	* VERSION
+	* NEWS
+	Changed to 2.0.0-beta2, noted in news
+
+2018-04-21  Jesse Thaler  <jthaler@jthaler.net>
+
+	* AUTHORS: updated arxiv number for RecursiveSoftDrop
+
+2018-04-21  Gavin Salam  <gavin.salam@cern.ch>
+
+	* README: updated arxiv number for RecursiveSoftDrop & Bottom-up
+	Soft Drop.
+
+2018-04-04  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSoftDrop.cc (contrib):
+	fixed syntax of calls to structure_of<...>
+	(thanks to Attila Krasznahorkay)
+
+2018-01-24  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSoftDrop.cc:
+	for the (dafault) dynamical R0 implementation, the dynamical R0 is
+	evolved independently in each branch.
+
+2017-10-10  Jesse Thaler  <jthaler@jthaler.net>
+	* AUTHORS
+	Added mention of BUSD
+
+	* README
+	Some tweaks to the wording
+
+2017-10-11  Gregory Soyez  <soyez@fastjet.fr>
+
+	* IteratedSoftDrop.hh:
+	IteratedSoftDropInfo::size() and multiplicity() now return an
+	unsigned int instead of a double
+
+2017-10-10  Jesse Thaler  <jthaler@jthaler.net>
+	* AUTHORS
+	Updated journal reference for ISD
+
+	* example_isd.{cc,ref}:
+	Included soft drop multiplicity in example
+
+	* README
+	Added warning at top that documentation has not been updated
+
+	* VERSION
+	Changed to 2.0.0-beta1
+
+2017-10-10  Gregory Soyez  <soyez@fastjet.fr>
+
+	* example_isd.ref:
+	updated to reflect the bugfix below
+
+	* IteratedSoftDrop.cc:
+	fixed issue in description (was taking sqrt of -ve number when
+	there were no angular cut)
+
+	* NEWS:
+	drafted for RecursiveTools-2.0.0-beta1
+
+	* TODO:
+	updated list in view of a beta release
+
+	* example_isd.cc:
+	pointed to the right header for IteratedSoftDrop
+
+	* RecursiveSoftDrop.{hh,cc}:
+	* IteratedSoftDrop.{hh,cc}:
+	moved IteratedSoftDrop to its own file (replacing the old
+	implementation)
+
+	* example_advanced_usage.ref:
+	updated reference file following the fix below.
+
+2017-09-28  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSymmetryCutBase.cc:
+	when no substructure is found, keep the _symmetry, _delta_R and
+	_mu structure variables at -1.  This for example allows one to
+	trigger on substructure by checking if delta_R>=0.
+
+	Note: we could set it to 0 (as it was before) and trigger using
+	_delta_R>0 but there might be some genuine substructure exactly
+	collinear.
+
+2017-09-19  Gregory Soyez  <soyez@fastjet.fr>
+
+	* example_isd.ref:
+	updated to the latest version of the example
+
+2017-09-18  Gregory Soyez  <soyez@fastjet.fr>
+
+	* Makefile:
+	updating make check to use all the examples (failing on ISD as
+	expected, see below)
+
+	* example_bottomup_softdrop.cc:
+	fixed typo
+
+	* example_bottomup_softdrop.ref:
+	* example_recursive_softdrop.ref:
+	added reference output for this example
+
+	* RecursiveSoftDrop.{hh,cc}:
+	* RecursiveSymmetryCutBase.{cc,hh}:
+	moved the "all_prongs" method from the base structure t oa
+	standalone function in RecursiveSoftDrop.hh
+
+	In practice, this is irrelevant for mMDT and SD (since pieces()
+	gets the job done, and the substructure class does not have (as
+	is) reliable info to get the full structure)
+
+	* RecursiveSymmetryCutBase.cc:
+	revamped a series of requests for substructure info to better
+	handle possible recursion into deeper jet substructure.
+
+	* RecursiveSoftDrop.{hh,cc}:
+	updated "description" to reuse the info from the base class
+
+	* example_isd.cc:
+	updated to use the newer implementation of ISD. Checked that it
+	gives the same results as the former implementation.
+
+	Note: make check still needs fixing because the example now
+	      computes a different set of angularities
+
+	* RecursiveSoftDrop.hh:
+	added a few helpers to IteratedSoftDropInfo ([] operator and size,
+	meaning it can be used as a vector<pair<double,double> >)
+
+	* RecursiveSymmetryCutBase.cc:
+	. fixed bugs in the calculation of the geometric distances for ee
+	  coordinates
+	. fixed bug in the computation of the (zg,thetag) pairs [it was
+	  returning the groomed ones instead of the ones passing the
+	  condition]
+
+	* example_recursive_softdrop.cc:
+	set the R0 parameter to the original jet radius
+
+2017-09-13  Gregory Soyez  <soyez@fastjet.fr>
+
+	* example_recursive_softdrop.cc:
+	tied up a few comments and the code output
+
+	* RecursiveSymmetryCutBase.{hh,cc}:
+	removed the unneeded _is_composite
+
+	* RecursiveSoftDrop.cc:
+	fixed issue with "verbose" dropped info on branches with no
+	further substructure
+
+2017-09-11  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSoftDrop.{hh,cc}:
+	have IteratedSoftDDrop returning by default an object of type
+	IteratedSoftDropInfo; added several helpers
+
+2017-09-08  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSoftDrop.{hh,cc}:
+	updated IteratedSoftDrop to give it the flexibility of
+	RecursiveSoftDrop
+
+	* RecursiveSymmetryCutBase.hh:
+	fixed typo in comment
+
+	* example_mmdt_ee.cc:     *** ADDED ***
+	added an example to illustrat usage in ee collisions
+
+	* example_isd.cc:
+	* BottomUpSoftDrop.cc:
+	* IteratedSoftDrop.cc:
+	* RecursiveSoftDrop.cc:
+	Fixed compilation issues with FJ 3.0 (mostly the usage of features
+	introduced only in FJ3.1)
+
+	* RecursiveSymmetryCutBase.{hh,cc}:
+	used the internal Recluster class for FJ<3.1.0 and the FJ antive
+	one for FJ>=3.1.0
+
+	* BottomUpSoftDrop.{hh,cc}:
+	moved the implementation of global_grooming to the source file and
+	fixed a few trivial compilation errors
+
+2017-09-08  FrÃ©dÃ©ric Dreyer  <frederic.dreyer@gmail.com>
+
+	* BottomUpSoftDrop.hh:
+	added the global_grooming method to process full event
+
+2017-09-07  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSoftDrop.cc:
+	cleaned (mostly by removing older commented-out code)
+
+	* RecursiveSoftDrop.{hh,cc}:
+	* RecursiveSymmetryCutBase.{hh,cc}:
+	* SoftDrop.cc:
+
+	added support for ee coordinates.  For that, the symmetry measure
+	has to be set to either theta_E (which uses the 3-vector angle
+	theta) or to cos_theta_E which uses sqrt(2*[1-cos(theta)])
+
+	Accordingly, the recursion_choice can be set to larger_E to
+	recurse into the branch with the largest energy. The larger_m
+	mode, recorsing into the larger-mass branch is also possible but
+	not advised (for the same reason as the pp version).
+
+	* RecursiveSymmetryCutBase.{hh,cc}:
+	switched to the Recluster class provided with FastJet. ASlso
+	included the recluster description to RecursiveSymmetryCutBase
+	when it is user-defined.
+
+2017-09-06  Gregory Soyez  <soyez@fastjet.fr>
+
+	* BottomUpSoftDrop.{hh,cc}:
+	. renamed SoftDropStructure  -> BottomUpSoftDropStructure
+	          SoftDropRecombiner -> BottomUpSoftDropRecombiner
+	          SoftDropPlugin     -> BottomUpSoftDropPlugin
+	. moved 'description' to source file (instead of header)
+	. kept the "area" information when available (jets will just
+	appear as having a 0 area)
+	. added most class description (main class still missing)
+
+	* RecursiveSoftDrop.cc:
+	. put internal classes to handle CS history in an internal namespace
+	. replaced the "switch" in the mail loop by a series of if (allows
+	  us a few simplificatins/cleaning)
+	. more uniform treatment of issues in the presence of an angular cut
+	  (as part of the above reorganisation)
+
+	* example_advanced_usage.ref:
+	updated reference output file following the bugfix (missing
+	"grooming mode" initialisation in one of the SoftDrop ctors) on
+	2017-08-01
+
+	* RecursiveSymmetryCutBase.cc:
+	removed redundent code
+
+2017-08-10  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSoftDrop.cc:
+	fixed trivial typo in variable name
+
+>>>>>>> .r1071
+2017-08-04  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSoftDrop.cc:
+	do not impose an angular cut in IterativeSD if it is -ve.
+
+2017-08-01  Gregory Soyez  <soyez@fastjet.fr>
+
+	* example_recursive_softdrop.cc:
+	added a series of optional flags
+
+	* RecursiveSoftDrop.cc:
+	fixed a few issues with the fixed depth version
+
+	* RecursiveSymmetryCutBase.hh:
+	a jet is now considered as havig substructure if deltaR>0
+	(coherent with released version)
+
+	* SoftDrop.hh:
+	bugfix: set the "grooming mode" by default in all ctors
+	EDIT: caused issue with make check, fixed on 2017-09-069 (see
+	      above)
+
+	* RecursiveSoftDrop.{hh,cc}:
+	added support for the "same depth" variant
+
+	* RecursiveSymmetryCutBase.cc:
+	also associate a RecursiveSymmetryCutBase::StructureType structure
+	to the result jet in case it is just a single particle (in
+	grooming mode)
+
+2017-07-31  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSymmetryCutBase.{hh,cc}:
+	added the option to pass an extra parameter to the symmetry cut
+	function
+
+	* RecursiveSymmetryCutBase.{hh,cc}:
+	* ModifiedMassDropTagger.hh
+	* SoftDrop.hh:
+	minor adaptions due to the above change + added a few methods to
+	query the class information (symmetry cut, beta, ...)
+
+	* RecursiveSoftDrop.{hh,cc}:
+	Added support for
+	- a dynamical R0
+	- recursing only in the hardest branch
+	- imposing a min deltaR cut
+
+	Added a tentative IterativeSoftDrop class
+
+2017-07-28  Gregory Soyez  <soyez@fastjet.fr>
+
+	* RecursiveSoftDrop.cc:
+	adapted to the latest changes in RecursiveSymmetryCutBase
+
+	* RecursiveSymmetryCutBase.{hh,cc}:
+	reorganised the output of the recursion step (recurse_one_step)
+	using an enum to make iot more readable (and to fix issues where
+	the dropped prong is actually 0, e.g. after subtraction)
+
+2017-07-26  Gregory Soyez  <soyez@fastjet.fr>
+
+	* example_recursive_softdrop.cc:   *** ADDED ***
+	added a draft example for the use of RecursiveSoftDrop
+
+	* RecursiveSoftDrop.{hh,cc}:       *** ADDED ***
+	added a first pass at an implementation of RecursiveSoftDrop.
+
+	This is based on Frederic's initial version but keeps the
+	branching structure of the jet. Some of the features, like
+	dynamical R0, direct access to the substructure or the same depth
+	variant, are still unsupported.
+
+	* SoftDrop.hh:
+	declared _beta, _symmetry_cut and _R0sqr as protected (was
+	private) so they ca n be used by RecursiveSoftDrop
+
+	* RecursiveSymmetryCutBase.{hh,cc}:
+	extracted from result() the part that performs one step of the
+	resursion (implemented as recurse_one_step()). This is repeatedly
+	called by result(). It has specific conventions to indicate
+	whether or not some substructure has been found or if one ran into
+	some issue.
+
+2017-04-25  Kevin Zhou  <knzhou@mit.edu>
+
+	* IteratedSoftDrop.hh
+	  . Added Doxygen documentation
+
+	* RecursiveSymmetryCutBase.hh
+	  . Added references to ISD
+
+2017-04-25  Jesse Thaler  <jthaler@jthaler.net>
+	* AUTHORS, COPYING:
+	  . Added ISD arXiv number
+
+	* example_isd.{cc,ref}
+	  . Added ISD arXiv number
+	  . Changing z_cut to be optimal value (with Lambda = 1 GeV)
+	  . Tweaked formatting
+
+	* IteratedSoftDrop.{hh,cc}
+	  . Added ISD arXiv number
+	  . Assert added if recluster does not return one jet.
+
+	* README
+	  . Added ISD arXiv number and tweaked wording
+
+
+2017-04-20  Kevin Zhou  <knzhou@mit.edu>
+
+	* example_isd.{cc,ref}      ** ADDED **
+	* IteratedSoftDrop.{cc,hh}  ** ADDED **
+
+	* Makefile
+	  . Added IteratedSoftDrop (ISD) as appropriate
+
+	* AUTHORS
+	  . Added myself to author list
+          . Added placeholder for ISD paper
+
+	* COPYING
+	  . Added placeholder for ISD paper
+
+	* README
+	  . Added description of ISD
+
+	* TODO
+	  . Added tasks to integrate ISD with other classes, e.g. SD, 
+	    Recluster, and the future RecursiveSoftDrop (RSD)
+
+	* NEWS
+	  . Added dummy for release of 1.1.0
+
+	* VERSION:
+	  . Switched version number to 1.1.0-dev
+
+	* example_advanced_usage.cc:
+	  . Added $Id$ tag, didn't change anything else
+
+2014-07-30  Gregory Soyez  <soyez@fastjet.fr>
+
+	* Recluster.hh: fixed the name of the #define for the header
+
 2014-07-09  Gregory Soyez  <soyez@fastjet.fr> + Jesse
 
Index: external/fastjet/contribs/RecursiveTools/IteratedSoftDrop.cc
===================================================================
--- external/fastjet/contribs/RecursiveTools/IteratedSoftDrop.cc	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
+++ external/fastjet/contribs/RecursiveTools/IteratedSoftDrop.cc	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -0,0 +1,128 @@
+// $Id: IteratedSoftDrop.cc 1084 2017-10-10 20:36:50Z gsoyez $
+//
+// Copyright (c) 2017-, Jesse Thaler, Kevin Zhou, Gavin P. Salam
+// andGregory Soyez
+//
+// based on arXiv:1704.06266 by Christopher Frye, Andrew J. Larkoski, 
+// Jesse Thaler, Kevin Zhou
+// 
+//----------------------------------------------------------------------
+// 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 "IteratedSoftDrop.hh"
+#include <sstream>
+#include <cmath>
+#include <fastjet/ClusterSequence.hh>
+
+using namespace std;
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+//========================================================================
+// implementation of IteratedSoftDropInfo
+//========================================================================
+
+// returns the angularity with angular exponent alpha and z
+// exponent kappa calculated on the zg's and thetag's found by
+// iterated SoftDrop
+//
+// returns 0 if no substructure was found
+double IteratedSoftDropInfo::angularity(double alpha, double kappa) const{
+  double sum = 0.0;
+  for (unsigned int i=0; i< _all_zg_thetag.size(); ++i)
+    sum += pow(_all_zg_thetag[i].first, kappa) * pow(_all_zg_thetag[i].second, alpha);
+  return sum;    
+}
+  
+//========================================================================
+// implementation of IteratedSoftDrop
+//========================================================================
+
+// Constructor. Takes in the standard Soft Drop parameters, an angular cut \f$\theta_{\rm cut}\f$, 
+// and a choice of angular and symmetry measure. 
+//
+//  - beta           the Soft Drop beta parameter
+//  - symmetry_cut   the Soft Drop symmetry cut
+//  - angular_cut    the angular cutoff to halt Iterated Soft Drop
+//  - R0             the angular distance normalization
+IteratedSoftDrop::IteratedSoftDrop(double beta, double symmetry_cut, double angular_cut, double R0,
+                                   const FunctionOfPseudoJet<PseudoJet> * subtractor) :
+  _rsd(beta, symmetry_cut, -1, R0, subtractor){
+  _rsd.set_hardest_branch_only(true);
+  if (angular_cut>0)
+    _rsd.set_min_deltaR_squared(angular_cut*angular_cut);
+}
+
+
+// Full constructor, which takes the following parameters:
+//
+// \param beta               the value of the beta parameter
+// \param symmetry_cut       the value of the cut on the symmetry measure
+// \param symmetry_measure   the choice of measure to use to estimate the symmetry
+// \param angular_cut        the angular cutoff to halt Iterated Soft Drop
+// \param R0                 the angular distance normalisation [1 by default]
+// \param mu_cut             the maximal allowed value of mass drop variable mu = m_heavy/m_parent 
+// \param recursion_choice   the strategy used to decide which subjet to recurse into
+// \param subtractor         an optional pointer to a pileup subtractor (ignored if zero)
+IteratedSoftDrop::IteratedSoftDrop(double  beta,
+                                   double  symmetry_cut, 
+                                   RecursiveSoftDrop::SymmetryMeasure  symmetry_measure,
+                                   double  angular_cut,
+                                   double  R0,
+                                   double  mu_cut, 
+                                   RecursiveSoftDrop::RecursionChoice  recursion_choice,
+                                   const FunctionOfPseudoJet<PseudoJet> * subtractor)
+  : _rsd(beta, symmetry_cut, symmetry_measure, -1, R0, mu_cut, recursion_choice, subtractor){
+  _rsd.set_hardest_branch_only(true);
+  if (angular_cut>0)
+    _rsd.set_min_deltaR_squared(angular_cut*angular_cut);
+}
+
+  
+// returns vector of ISD symmetry factors and splitting angles
+IteratedSoftDropInfo IteratedSoftDrop::result(const PseudoJet& jet) const{
+  PseudoJet rsd_jet = _rsd(jet);
+  if (! rsd_jet.has_structure_of<RecursiveSoftDrop>())
+    return IteratedSoftDropInfo();
+  return IteratedSoftDropInfo(rsd_jet.structure_of<RecursiveSoftDrop>().sorted_zg_and_thetag());
+}
+
+
+std::string IteratedSoftDrop::description() const{
+  std::ostringstream oss;
+  oss << "IteratedSoftDrop with beta =" << _rsd.beta()
+      << ", symmetry_cut=" << _rsd.symmetry_cut()
+      << ", R0=" << _rsd.R0();
+  
+  if (_rsd.min_deltaR_squared() >= 0){
+    oss << " and angular_cut=" << sqrt(_rsd.min_deltaR_squared());
+  } else {
+    oss << " and no angular_cut";
+  }
+  
+  if (_rsd.subtractor()){
+    oss << ", and with internal subtraction using [" << _rsd.subtractor()->description() << "]";
+  }
+  return oss.str();
+}
+
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/RecursiveTools/IteratedSoftDrop.hh
===================================================================
--- external/fastjet/contribs/RecursiveTools/IteratedSoftDrop.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
+++ external/fastjet/contribs/RecursiveTools/IteratedSoftDrop.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -0,0 +1,219 @@
+// $Id: IteratedSoftDrop.hh 1086 2017-10-11 08:07:26Z gsoyez $
+//
+// Copyright (c) 2017-, Jesse Thaler, Kevin Zhou, Gavin P. Salam,
+// Gregory Soyez
+//
+// based on arXiv:1704.06266 by Christopher Frye, Andrew J. Larkoski,
+// Jesse Thaler, Kevin Zhou
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+
+#ifndef __FASTJET_CONTRIB_ITERATEDSOFTDROP_HH__
+#define __FASTJET_CONTRIB_ITERATEDSOFTDROP_HH__
+
+#include "RecursiveSoftDrop.hh"
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+//------------------------------------------------------------------------
+/// \class IteratedSoftDropInfo
+/// helper class that carries all the relevant information one can get
+/// from running IteratedSoftDrop on a given jet (or vector of jets)
+///
+class IteratedSoftDropInfo{
+public:
+  /// ctor without initialisation
+  IteratedSoftDropInfo(){}
+
+  /// ctor with initialisation
+  IteratedSoftDropInfo(std::vector<std::pair<double,double> > zg_thetag_in)
+    : _all_zg_thetag(zg_thetag_in){}
+
+  /// get the raw list of (angular-ordered) zg and thetag
+  const std::vector<std::pair<double,double> > &all_zg_thetag() const{
+    return _all_zg_thetag;
+  }
+
+  /// overloadd the () operator so that it also returns the full (zg,thetag) list
+  const std::vector<std::pair<double,double> > & operator()() const{
+    return _all_zg_thetag;
+  }
+
+  /// overloadd the [] operator to access the ith (zg,thetag) pair
+  const std::pair<double,double> & operator[](unsigned int i) const{
+    return _all_zg_thetag[i];
+  }
+
+  /// returns the angularity with angular exponent alpha and z
+  /// exponent kappa calculated on the zg's and thetag's found by
+  /// iterated SoftDrop
+  ///
+  /// returns 0 if no substructure was found
+  double angularity(double alpha, double kappa=1.0) const;
+
+  /// returns the Iterated SoftDrop multiplicity
+  unsigned int multiplicity() const{ return _all_zg_thetag.size(); }  
+
+  /// returns the Iterated SoftDrop multiplicity (i.e. size)
+  unsigned int size() const{ return _all_zg_thetag.size(); }  
+  
+protected:
+  /// the real information: angular-ordered list of all the zg and
+  /// thetag that passed the (recursive) SD conddition
+  std::vector<std::pair<double,double> > _all_zg_thetag;
+};
+  
+
+  
+//------------------------------------------------------------------------
+/// \class IteratedSoftDrop
+/// implementation of the IteratedSoftDrop procedure
+/// 
+/// This class provides an implementation of the IteratedSoftDrop
+/// procedure.  It is based on the SoftDrop procedure can be used to
+/// define a 'groomed symmetry factor', equal to the symmetry factor
+/// of the two subjets of the resulting groomed jet.  The Iterated
+/// Soft Drop procedure recursively performs Soft Drop on the harder
+/// branch of the groomed jet, halting at a specified angular cut
+/// \f$\theta_{\rm cut}\f$, returning a list of symmetry factors which
+/// can be used to define observables.
+///
+/// Like SoftDrop, the cut applied recursively is 
+///   \f[
+///      z > z_{\rm cut} (\theta/R_0)^\beta
+///   \f]
+/// with z the asymmetry measure and \f$\theta\f$ the geometrical
+/// distance between the two subjets. The procedure halts when
+/// \f$\theta < \theta_{\rm cut}\f$.
+///
+/// By default, this implementation returs the IteratedSoftDropInfo
+/// obtained after running IteratedSoftDrop on a jet
+///
+/// Although all these quantities can be obtained from the returned
+/// IteratedSoftDropInfo, we also provide helpers to directly get the
+/// multiplicity, some (generalised) angularity, or the raw list of
+/// (angular-ordered) (zg, thetag) pairs that passed the (recursive)
+/// SoftDrop condition.
+///
+/// We stress the fact that IteratedSoftDrop is _not_ a Transformer
+/// since it returns an IteratedSoftDropInfo and not a modified
+/// PseudoJet
+///
+class IteratedSoftDrop : public FunctionOfPseudoJet<IteratedSoftDropInfo> {
+public:
+  /// Constructor. Takes in the standard Soft Drop parameters, an angular cut \f$\theta_{\rm cut}\f$, 
+  /// and a choice of angular and symmetry measure. 
+  ///
+  /// \param beta               the Soft Drop beta parameter
+  /// \param symmetry_cut       the Soft Drop symmetry cut
+  /// \param angular_cut        the angular cutoff to halt Iterated Soft Drop
+  /// \param R0                 the angular distance normalization
+  /// \param subtractor         an optional pointer to a pileup subtractor (ignored if zero)
+  IteratedSoftDrop(double beta, double symmetry_cut, double angular_cut, double R0 = 1.0,
+                   const FunctionOfPseudoJet<PseudoJet> * subtractor = 0);
+  
+  /// Full constructor, which takes the following parameters:
+  ///
+  /// \param beta               the value of the beta parameter
+  /// \param symmetry_cut       the value of the cut on the symmetry measure
+  /// \param symmetry_measure   the choice of measure to use to estimate the symmetry
+  /// \param angular_cut        the angular cutoff to halt Iterated Soft Drop
+  /// \param R0                 the angular distance normalisation [1 by default]
+  /// \param mu_cut             the maximal allowed value of mass drop variable mu = m_heavy/m_parent 
+  /// \param recursion_choice   the strategy used to decide which subjet to recurse into
+  /// \param subtractor         an optional pointer to a pileup subtractor (ignored if zero)
+  ///
+  /// Notes: 
+  ///
+  /// - by default, SoftDrop will recluster the jet with the
+  ///   Cambridge/Aachen algorithm if it is not already the case. This
+  ///   behaviour can be changed using the "set_reclustering" method
+  ///   defined below
+  ///
+  IteratedSoftDrop(double  beta,
+                   double  symmetry_cut, 
+                   RecursiveSoftDrop::SymmetryMeasure  symmetry_measure,
+                   double  angular_cut,
+                   double  R0 = 1.0,
+                   double  mu_cut = std::numeric_limits<double>::infinity(), 
+                   RecursiveSoftDrop::RecursionChoice  recursion_choice = RecursiveSoftDrop::larger_pt,
+                   const FunctionOfPseudoJet<PseudoJet> * subtractor = 0);
+
+  /// default destructor
+  virtual ~IteratedSoftDrop(){}
+
+  //----------------------------------------------------------------------
+  // behaviour tweaks (inherited from RecursiveSoftDrop and RecursiveSymmetryCutBase)
+
+  /// switch to using a dynamical R0 (see RecursiveSoftDrop)
+  void set_dynamical_R0(bool value=true) { _rsd.set_dynamical_R0(value); }
+  bool use_dynamical_R0() const { return _rsd.use_dynamical_R0(); }
+
+  /// an alternative way to set the subtractor (see RecursiveSymmetryCutBase)
+  void set_subtractor(const FunctionOfPseudoJet<PseudoJet> * subtractor_) {_rsd.set_subtractor(subtractor_);}
+  const FunctionOfPseudoJet<PseudoJet> * subtractor() const {return _rsd.subtractor();}
+  
+  /// returns the IteratedSoftDropInfo associated with the jet "jet"
+  IteratedSoftDropInfo result(const PseudoJet& jet) const;
+
+  /// Tells the tagger whether to assume that the input jet has
+  /// already been subtracted (relevant only with a non-null
+  /// subtractor, see RecursiveSymmetryCutBase)
+  void set_input_jet_is_subtracted(bool is_subtracted) { _rsd.set_input_jet_is_subtracted(is_subtracted);}
+  bool input_jet_is_subtracted() const {return _rsd.input_jet_is_subtracted();}
+  
+  /// configure the reclustering prior to the recursive de-clustering
+  void set_reclustering(bool do_reclustering=true, const Recluster *recluster=0){
+    _rsd.set_reclustering(do_reclustering, recluster);
+  }
+  
+  //----------------------------------------------------------------------
+  // actions on jets
+  /// returns vector of ISD symmetry factors and splitting angles
+  std::vector<std::pair<double,double> > all_zg_thetag(const PseudoJet& jet) const{
+    return result(jet).all_zg_thetag();
+  }
+
+  /// returns the angularity with angular exponent alpha and z
+  /// exponent kappa calculated on the zg's and thetag's found by
+  /// iterated SoftDrop
+  ///
+  /// returns 0 if no substructure was found
+  double angularity(const PseudoJet& jet, double alpha, double kappa=1.0) const{
+    return result(jet).angularity(alpha, kappa);
+  }
+  
+  /// returns the Iterated SoftDrop multiplicity
+  double multiplicity(const PseudoJet& jet) const{ return result(jet).multiplicity(); }
+
+  /// description of the class
+  std::string description() const;
+
+protected:
+  RecursiveSoftDrop _rsd;
+};
+
+  
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif  // __FASTJET_CONTRIB_ITERATEDSOFTDROP_HH__
Index: external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.hh
===================================================================
--- external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.hh	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/ModifiedMassDropTagger.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,3 +1,3 @@
-// $Id: ModifiedMassDropTagger.hh 688 2014-06-17 14:29:56Z jthaler $
+// $Id: ModifiedMassDropTagger.hh 1032 2017-07-31 14:20:03Z gsoyez $
 //
 // Copyright (c) 2014-, Gavin P. Salam
@@ -105,4 +105,8 @@
   virtual ~ModifiedMassDropTagger(){}
 
+  //----------------------------------------------------------------------
+  // access to class info
+  double symmetry_cut() const { return _symmetry_cut; }
+  
 protected:
 
@@ -110,5 +114,6 @@
   /// has no dependence on the subjet kinematics
   virtual double symmetry_cut_fn(const PseudoJet & /* p1 */, 
-                                 const PseudoJet & /* p2 */
+                                 const PseudoJet & /* p2 */,
+                                 void *extra_parameters = 0
                                  ) const {return _symmetry_cut;}
   virtual std::string symmetry_cut_description() const;
Index: external/fastjet/contribs/RecursiveTools/NEWS
===================================================================
--- external/fastjet/contribs/RecursiveTools/NEWS	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/NEWS	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,2 +1,11 @@
+2018/05/31: release of version 2.0.0-beta2 with corrected syntax
+
+2017/10/10: release of version 2.0.0-beta1 including implementations of
+
+* RecursiveSoftDrop (see example_rsd.hh for usage)
+* IteratedSoftDrop (see example_isd.hh for usage)
+* e+e- version of the recursive tools (see example_mmdt_ee.hh for usage)
+* BottomUpSoftDrop (see example_bottomup_softdrop.cc for usage)
+
 2014/07/09: release of version 1.0.0 of RecursiveTools including
             ModifiedMassDropTagger and SoftDrop (as well as Recluster)
Index: external/fastjet/contribs/RecursiveTools/README
===================================================================
--- external/fastjet/contribs/RecursiveTools/README	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/README	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -17,8 +17,19 @@
   Marzani, Gregory Soyez, Jesse Thaler
 
+- RecursiveSoftDrop
+- BottomUpSoftDrop
+  This corresponds to arXiv:1804.03657 by Frederic Dreyer, Lina
+  Necib, Gregory Soyez and Jesse Thaler
+
+- IteratedSoftDrop
+  This corresponds to arXiv:1704.06266 by Christopher Frye, Andrew J.
+  Larkoski, Jesse Thaler, Kevin Zhou
+
 - Recluster
   A generic tool to recluster a given jet into subjets
-  Note: this is largely based on the Filter code in FastJet v3.0 and
-  ultimately, this tool will probably be moved into FastJet
+  Note: a Recluster class is available natively in FastJet since v3.1.
+        Users are therefore encouraged to use the FastJet version
+        rather than this one which is mostly provided for
+        compatibility of this contrib with older versions of FastJet.
 
 The interface for these tools is described in more detail below, with
@@ -74,5 +85,5 @@
 
 The SoftDrop procedure is very similar to mMDT, albeit with a 
-generalized symmetry condition:
+generalised symmetry condition:
  
    z > z_cut * (R / R0)^beta
@@ -108,11 +119,116 @@
 
 ------------------------------------------------------------------------
+RecursiveSoftDrop
+------------------------------------------------------------------------
+
+The RecursiveSoftDrop procedure applies the Soft Drop procedure N times 
+in a jet in order to find up to N+1 prongs.  N=0 makes no modification
+to the jet, and N=1 is equivalent to the original SoftDrop.
+
+Once one has more than one prong, one has to decide which will be
+declustered next.  At each step of the declustering procedure, one
+undoes the clustering which has the largest declustering angle
+(amongst all the branches that are searched for substructure). [see
+"set_fixed_depth" below for an alternative]
+
+Compared to SoftDrop, RecursiveSoftDrop takes an extra argument N
+specifying the number of times the SoftDrop procedure is recursively
+applied. Negative N means that the procedure is applied until no
+further substructure is found (i.e. corresponds to taking N=infinity).
+
+   double z_cut = 0.10;
+   double beta  = 2.0;
+   double R0    = 1.0; // this is the default value
+   int N        = -1;  
+   RecursiveSoftDrop rsd(z_cut, beta, N, R0);
+
+One then acts on a jet as
+
+   PseudoJet groomed_jet = rsd(jet)
+
+and get additional information via
+
+   groomed_jet.structure_of<RecursiveSoftDrop>()
+
+------------------------------------------------------------------------
+IteratedSoftDrop
+------------------------------------------------------------------------
+
+Iterated Soft Drop (ISD) is a repeated variant of SoftDrop.  After
+performing the Soft Drop procedure once, it logs the groomed symmetry
+factor, then recursively performs Soft Drop again on the harder
+branch.  This procedure is repeated down to an (optional) angular cut
+theta_cut, yielding a set of symmetry factors from which observables
+can be built.
+
+An IteratedSoftDrop tool can be created as follows:
+
+  double beta = -1.0;
+  double z_cut = 0.005;
+  double theta_cut = 0.0;
+  double R0 = 0.5;        // characteristic radius of jet algorithm
+  IteratedSoftDrop isd(beta, z_cut, double theta_cut, R0);
+
+By default, ISD applied on a jet gives a result of type
+IteratedSoftDropInfo that can then be probed to obtain physical
+observables
+
+  IteratedSoftDropInfo isd_info = isd(jet);
+
+  unsigned int multiplicity = isd_info.multiplicity();
+  double kappa = 1.0;     // changes angular scale of ISD angularity
+  double isd_width = isd_info.angularity(kappa);
+  vector<pair<double,double> > zg_thetags = isd_info.all_zg_thetag();
+  vector<pair<double,double> > zg_thetags = isd_info();
+  for (unsigned int i=0; i< isd_info.size(); ++i){
+    cout << "(zg, theta_g)_" << i << " = "
+         << isd_info[i].first << " " <<  isd_info[i].second << endl;
+  }
+
+Alternatively, one can directly get the multiplicity, angularity, and
+(zg,thetag) pairs from the IteratedSoftDrop class, at the expense of
+re-running the declustering procedure:
+
+  unsigned int multiplicity = isd.multiplicity(jet);
+  double isd_width = isd.angularity(jet, 1.0);
+  vector<pair<double,double> > zg_thetags = isd.all_zg_thetag(jet);
+
+
+Note: the iterative declustering procedure is the same as what one
+  would obtain with RecursiveSoftDrop with an (optional) angular cut
+  and recursing only in the hardest branch [see the "Changing
+  behaviour" section below for details], except that it returns some
+  information about the jet instead of a modified jet as RSD does.
+
+
+------------------------------------------------------------------------
+BottomUpSoftDrop
+------------------------------------------------------------------------
+
+This is a bottom-up version of the RecursiveSoftDrop procedure, in a
+similar way as Pruning can be seen as a bottom-up version of Trimming.
+
+In practice, the jet is reclustered and at each step of the clustering
+one checks the SoftDrop condition
+
+   z > z_cut * (R / R0)^beta
+
+If the condition is met, the pair is recombined. If the condition is
+not met, only the hardest of the two objects is kept for further
+clustering and the softest is rejected.
+
+------------------------------------------------------------------------
 Recluster
 ------------------------------------------------------------------------
+
+  *** NOTE: this is provided only for backwards compatibility ***
+  *** with FastJet <3.1. For FastJet >=3.1, the native        ***
+  *** fastjet::Recluster is used instead                      ***
 
 The Recluster class allows the constituents of a jet to be reclustered
 with a different recursive clustering algorithm.  This is used
-internally in the mMDT/SoftDrop code in order to recluster the jet using
-the CA algorithm.  This is achieved via
+internally in the mMDT/SoftDrop/RecursiveSoftDrop/IteratedSoftDrop
+code in order to recluster the jet using the CA algorithm.  This is
+achieved via
 
   Recluster ca_reclusterer(cambridge_algorithm,
@@ -123,24 +239,37 @@
 delete_self_when_unused.
 
-
 ------------------------------------------------------------------------
 Changing behaviour
 ------------------------------------------------------------------------
 
-The behaviour of the ModifiedMassDropTagger and SoftDrop classes can
-be tweaked using the following options:
-
-SymmetryMeasure = {scalar_z,vector_z,y}  [constructor argument]
+The behaviour of the all the tools provided here
+(ModifiedMassDropTagger, SoftDrop, RecursiveSoftDrop and
+IteratedSoftDrop) can be tweaked using the following options:
+
+SymmetryMeasure = {scalar_z, vector_z, y, theta_E, cos_theta_E}
+  [constructor argument]
   : The definition of the energy sharing between subjets, with 0
-    corresponding to the most asymmetric
-
-RecursionChoice = {larger_pt,larger_mt,larger_m}  [constructor argument]
+    corresponding to the most asymmetric.
+    . scalar_z = min(pt1,pt2)/(pt1+pt2)   [default]
+    . vector_z = min(pt1,pt2)/pt_{1+2}
+    . y = min(pt1^2,pt2^2)/m_{12}^2  (original y from MDT)
+    . theta_E     = min(E1,E2)/(E1+E2),
+           with angular measure theta_{12}^2
+    . cos_theta_E = min(E1,E2)/(E1+E2),
+           with angular measure 2[1-cos(theta_{12})]
+    The last two variants are meant for use in e+e- collisions,
+    together with the "larger_E" recursion choice (see below)
+
+RecursionChoice = {larger_pt, larger_mt, larger_m, larger_E}
+  [constructor argument]
   : The path to recurse through the tree after the symmetry condition
-    fails
+    fails. Options refer to transverse momentum (pt), transverse mass
+    (mt=sqrt(pt^2+m^2), mass (m) or energy (E). the latter is meant
+    for use in e+e- collisions
 
 mu_cut   [constructor argument]
   : An optional mass drop condition
 
-set_subtractor(subtractor*) [or subtracter as a constructor argument]
+set_subtractor(subtractor*) [or subtractor as a constructor argument]
   : provide a subtractor. When a subtractor is supplied, the
     kinematic constraints are applied on subtracted 4-vectors. In
@@ -165,4 +294,35 @@
     and SoftDrop defaults to grooming mode.
 
+set_verbose_structure(bool)
+  : when set to true, additional information will be stored in the jet
+    structure. This includes in particular values of symmetry,
+    delta_R, and mu of dropped branches
+
+For the specific case of RecursiveSoftDrop, additional tweaking is
+possible via the following methods
+
+set_fixed_depth_mode(bool)
+  : when this is true, RSD will recurse (N times) into all the
+    branches found during the previous iteration [instead of recursing
+    through the largest declustering angle until N prongs have been
+    found]. This yields at most 2^N prong. For infinite N, the two
+    options are equivalent.
+
+set_dynamical_R0(bool)
+  : By default the angles in the SD condition are normalised to the
+    parameter R0. With "dynamical R0", RSD will dynamically adjust R0
+    to be the angle between the two prongs found during the previous
+    iteration.
+
+set_hardest_branch_only(bool)
+  : When substructure is found, only recurse into the hardest of the
+    two branches for further substructure search. This uses the class
+    RecursionChoice.
+
+set_min_deltaR_squared(double):
+  : set a minimal angle (squared) at which we stop the declustering
+    procedure.  This cut is ineffective for negative values of the
+    argument.
+
 ------------------------------------------------------------------------
 Technical Details
Index: external/fastjet/contribs/RecursiveTools/Recluster.hh
===================================================================
--- external/fastjet/contribs/RecursiveTools/Recluster.hh	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/Recluster.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,6 +1,6 @@
-#ifndef __FASTJET_TOOLS_RECLUSTER_HH__
-#define __FASTJET_TOOLS_RECLUSTER_HH__
+#ifndef __FASTJET_CONTRIB_TOOLS_RECLUSTER_HH__
+#define __FASTJET_CONTRIB_TOOLS_RECLUSTER_HH__
 
-// $Id: Recluster.hh 700 2014-07-07 12:50:05Z gsoyez $
+// $Id: Recluster.hh 723 2014-07-30 09:11:01Z gsoyez $
 //
 // Copyright (c) 2014-, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -181,3 +181,3 @@
 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
 
-#endif   // __FASTJET_TOOLS_RECLUSTER_HH__
+#endif   // __FASTJET_CONTRIB_TOOLS_RECLUSTER_HH__
Index: external/fastjet/contribs/RecursiveTools/RecursiveSoftDrop.cc
===================================================================
--- external/fastjet/contribs/RecursiveTools/RecursiveSoftDrop.cc	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
+++ external/fastjet/contribs/RecursiveTools/RecursiveSoftDrop.cc	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -0,0 +1,482 @@
+// $Id: RecursiveSoftDrop.cc 1111 2018-04-04 10:06:11Z gsoyez $
+//
+// Copyright (c) 2017-, Gavin P. Salam, Gregory Soyez, Jesse Thaler,
+// Kevin Zhou, Frederic Dreyer
+//
+//----------------------------------------------------------------------
+// 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 "RecursiveSoftDrop.hh"
+#include "fastjet/ClusterSequence.hh"
+
+using namespace std;
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+namespace contrib{
+
+namespace internal_recursive_softdrop{
+  
+//========================================================================
+/// \class RSDHistoryElement
+/// a helper class to help keeping track od the RSD tree
+///
+/// The element is created at the top of a branch and updated each
+/// time one grooms something away.
+class RSDHistoryElement{
+public:
+  RSDHistoryElement(const PseudoJet &jet, const RecursiveSoftDrop* rsd_ptr, double R0sqr) :
+    R0_squared(R0sqr),
+    child1_in_history(-1), child2_in_history(-1), symmetry(-1.0), mu2(-1.0){
+    reset(jet, rsd_ptr);
+  }
+
+  void reset(const PseudoJet &jet, const RecursiveSoftDrop* rsd_ptr){
+    current_in_ca_tree = jet.cluster_hist_index();
+    PseudoJet piece1, piece2;
+    theta_squared = (jet.has_parents(piece1, piece2))
+      ? rsd_ptr->squared_geometric_distance(piece1,piece2) : 0.0;
+  }
+  
+  int current_in_ca_tree;  ///< (history) index of the current particle in the C/A tree 
+  double theta_squared;    ///< squared angle at which this decays
+  double R0_squared;       ///< squared angle at the top of the branch
+                           ///< (used for RSD with dynamic_R0)
+  int child1_in_history;   ///< hardest of the 2 decay products (-1 if untagged)
+  int child2_in_history;   ///< softest of the 2 decay products (-1 if untagged)
+
+  // info about what has been dropped and the local substructure
+  vector<double> dropped_delta_R;
+  vector<double> dropped_symmetry;
+  vector<double> dropped_mu;
+  double symmetry, mu2;
+};
+
+
+/// \class OrderRSDHistoryElements
+/// angular ordering of (pointers to) the history elements
+///
+/// our priority queue will use pointers to these elements that are
+/// ordered in angle (of the objects they point to)
+class OrderRSDHistoryElements{
+public:
+  bool operator()(const RSDHistoryElement *e1, const RSDHistoryElement *e2) const {
+    return e1->theta_squared < e2->theta_squared;
+  }
+};
+
+} // internal_recursive_softdrop
+  
+//========================================================================
+
+// initialise all the flags and parameters to their default value
+void RecursiveSoftDrop::set_defaults(){
+  set_fixed_depth_mode(false);
+  set_dynamical_R0(false);
+  set_hardest_branch_only(false);
+  set_min_deltaR_squared(-1.0);
+}
+
+// description of the tool
+string RecursiveSoftDrop::description() const{
+  ostringstream res;
+  res << "recursive application of ["
+      << RecursiveSymmetryCutBase::description()
+      << "]";
+  
+  if (_fixed_depth){
+    res << ", recursively applied down to a maximal depth of N=";
+    if (_n==-1) res << "infinity"; else res << _n;
+  } else {
+    res << ", applied N=";
+    if (_n==-1) res << "infinity"; else res << _n;
+    res << " times";
+  }
+  
+  if (_dynamical_R0)
+    res << ", with R0 dynamically scaled";
+  else
+    res << ", with R0 kept fixed";
+    
+  if (_hardest_branch_only)
+    res << ", following only the hardest branch";
+
+  if (_min_dR2>0)
+    res << ", with minimal angle (squared) = " << _min_dR2;
+
+  return res.str();
+}
+  
+  
+// action on a single jet with RecursiveSoftDrop.
+//
+// uses "result_fixed_tags" by default (i.e. recurse from R0 to
+// smaller angles until n SD conditions have been met), or
+// "result_fixed_depth" where each of the previous SD branches are
+// recirsed into down to a depth of n.
+PseudoJet RecursiveSoftDrop::result(const PseudoJet &jet) const{
+  return _fixed_depth ? result_fixed_depth(jet) : result_fixed_tags(jet);
+}
+
+// this routine applies the Soft Drop criterion recursively on the
+// CA tree until we find n subjets (or until it converges), and
+// adds them together into a groomed PseudoJet
+PseudoJet RecursiveSoftDrop::result_fixed_tags(const PseudoJet &jet) const {
+  // start by reclustering jet with C/A algorithm
+  PseudoJet ca_jet = _recluster_if_needed(jet);
+
+  if (! ca_jet.has_valid_cluster_sequence()){
+    throw Error("RecursiveSoftDrop can only be applied to jets associated to a (valid) cluster sequence");
+  }
+  
+  const ClusterSequence *cs = ca_jet.validated_cluster_sequence();
+  const vector<ClusterSequence::history_element> &cs_history = cs->history();
+  const vector<PseudoJet> &cs_jets = cs->jets();
+
+  // initialize counter to 1 subjet (i.e. the full ca_jet)
+  int n_tagged = 0;
+  int max_njet = ca_jet.constituents().size();
+
+  // create the list of branches
+  unsigned int max_history_size = 2*max_njet;
+  if ((_n>0) && (_n<max_njet-1)){ max_history_size = 2*(_n+1); }
+
+  // we need to pre-allocate the space for the vector so that the
+  // pointers are not invalidated
+  vector<internal_recursive_softdrop::RSDHistoryElement> history;
+  history.reserve(max_history_size);  // could be one shorter
+  history.push_back(internal_recursive_softdrop::RSDHistoryElement(ca_jet, this, _R0sqr));
+  
+  // create a priority queue containing the subjets and a comparison definition
+  // initialise to the full ca_jet
+  priority_queue<internal_recursive_softdrop::RSDHistoryElement*, vector<internal_recursive_softdrop::RSDHistoryElement*>, internal_recursive_softdrop::OrderRSDHistoryElements> active_branches;
+  active_branches.push(& (history[0]));
+
+  PseudoJet parent, piece1, piece2;
+  double sym, mu2;
+
+  // which R0 to use
+  //double R0sqr = _R0sqr;
+  
+  // loop over C/A tree until we reach the appropriate number of subjets
+  while ((continue_grooming(n_tagged)) && (active_branches.size())) {
+    // get the element corresponding to the max dR
+    // and the associated PJ
+    internal_recursive_softdrop::RSDHistoryElement * elm = active_branches.top();
+    PseudoJet parent = cs_jets[cs_history[elm->current_in_ca_tree].jetp_index];
+    
+    // do one step of SD
+    RecursionStatus status = recurse_one_step(parent, piece1, piece2, sym, mu2, &elm->R0_squared);
+
+    // check if we passed the SD condition
+    if (status==recursion_success){
+      // check for the optional angular cut
+      if ((_min_dR2 > 0) && (squared_geometric_distance(piece1,piece2) < _min_dR2))
+        break;
+
+      // both subjets are kept in the list for potential further de-clustering
+      elm->child1_in_history = history.size();
+      elm->child2_in_history = history.size()+1;
+      elm->symmetry = sym;
+      elm->mu2      = mu2;
+      active_branches.pop();
+
+      // update the history
+      double next_R0_squared = (_dynamical_R0)
+        ? piece1.squared_distance(piece2) : elm->R0_squared;
+
+      internal_recursive_softdrop::RSDHistoryElement elm1(piece1, this, next_R0_squared);
+      history.push_back(elm1);
+      active_branches.push(&(history.back()));
+      internal_recursive_softdrop::RSDHistoryElement elm2(piece2, this, next_R0_squared);
+      history.push_back(elm2);
+      if (!_hardest_branch_only){
+        active_branches.push(&(history.back()));
+      }
+      
+      ++n_tagged;
+    } else if (status==recursion_dropped){
+       // check for the optional angular cut
+      if ((_min_dR2 > 0) && (squared_geometric_distance(piece1,piece2) < _min_dR2))
+        break;
+
+      active_branches.pop();
+      // tagging failed and the softest branch should be dropped
+      // keep track of what has een groomed away
+      max_njet -= piece2.constituents().size();
+      elm->dropped_delta_R .push_back((elm->theta_squared >= 0) ? sqrt(elm->theta_squared) : -sqrt(elm->theta_squared));
+      elm->dropped_symmetry.push_back(sym);
+      elm->dropped_mu      .push_back((mu2>=0) ? sqrt(mu2) : -sqrt(mu2));
+      
+      // keep the hardest bhanch in the recursion
+      elm->reset(piece1, this);
+      active_branches.push(elm);
+    } else if (status==recursion_no_parents){
+      if (_min_dR2 > 0) break;
+      active_branches.pop();
+      // nothing specific to do: we just keep the curent jet as a "leaf"
+    } else { // recursion_issue
+      active_branches.pop();
+      // we've met an issue
+      // if the piece2 is null as well, it means we've had a critical problem.
+      // In that case, return an empty PseudoJet
+      if (piece2 == 0) return PseudoJet();
+
+      // otherwise, we should consider "piece2" as a final particle
+      // not to be recursed into
+      if (_min_dR2 > 0) break;
+      max_njet -= (piece2.constituents().size()-1);
+      break;
+    }
+    
+    // If the missing number of tags is exactly the number of objects
+    // we have left in the recursion, stop
+    if (n_tagged == max_njet) break;
+  }
+
+  // now we have a bunch of history elements that we can use to build the final jet
+  vector<PseudoJet> mapped_to_history(history.size());
+  unsigned int history_index = history.size();
+  do {
+    --history_index;
+    const internal_recursive_softdrop::RSDHistoryElement & elm = history[history_index];
+
+    // two kinds of events: either just a final leave, poteitially  with grooming
+    // or a brandhing (also with potential grooming at the end)
+    if (elm.child1_in_history<0){
+      // this is a leaf, i.e. with no further sustructure
+      PseudoJet & subjet = mapped_to_history[history_index]
+        = cs_jets[cs_history[elm.current_in_ca_tree].jetp_index];
+
+      StructureType * structure = new StructureType(subjet);
+      if (has_verbose_structure()){
+        structure->set_verbose(true);
+        structure->set_dropped_delta_R (elm.dropped_delta_R);
+        structure->set_dropped_symmetry(elm.dropped_symmetry);
+        structure->set_dropped_mu      (elm.dropped_mu);
+      }
+      subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
+    } else {
+      PseudoJet & subjet = mapped_to_history[history_index]
+        = join(mapped_to_history[elm.child1_in_history], mapped_to_history[elm.child2_in_history]);
+      StructureType * structure = new StructureType(subjet, sqrt(elm.theta_squared), elm.symmetry, sqrt(elm.mu2));
+      if (has_verbose_structure()){
+        structure->set_verbose(true);
+        structure->set_dropped_delta_R (elm.dropped_delta_R);
+        structure->set_dropped_symmetry(elm.dropped_symmetry);
+        structure->set_dropped_mu      (elm.dropped_mu);
+      }
+      subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
+    }    
+  } while (history_index>0);
+
+  return mapped_to_history[0];
+}
+
+// this routine applies the Soft Drop criterion recursively on the
+// CA tree, recursing into all the branches found during the previous iteration
+// until n layers have been found (or until it converges)
+PseudoJet RecursiveSoftDrop::result_fixed_depth(const PseudoJet &jet) const {
+  // start by reclustering jet with C/A algorithm
+  PseudoJet ca_jet = _recluster_if_needed(jet);
+
+  if (! ca_jet.has_valid_cluster_sequence()){
+    throw Error("RecursiveSoftDrop can only be applied to jets associated to a (valid) cluster sequence");
+  }
+  
+  const ClusterSequence *cs = ca_jet.validated_cluster_sequence();
+  const vector<ClusterSequence::history_element> &cs_history = cs->history();
+  const vector<PseudoJet> &cs_jets = cs->jets();
+
+  // initialize counter to 1 subjet (i.e. the full ca_jet)
+  int n_depth = 0;
+  int max_njet = ca_jet.constituents().size();
+
+  // create the list of branches
+  unsigned int max_history_size = 2*max_njet;
+  //if ((_n>0) && (_n<max_njet-1)){ max_history_size = 2*(_n+1); }
+
+  // we need to pre-allocate the space for the vector so that the
+  // pointers are not invalidated
+  vector<internal_recursive_softdrop::RSDHistoryElement> history;
+  history.reserve(max_history_size);  // could be one shorter
+  history.push_back(internal_recursive_softdrop::RSDHistoryElement(ca_jet, this, _R0sqr));
+  history.back().theta_squared = _R0sqr;
+  
+  // create a priority queue containing the subjets and a comparison definition
+  // initialize counter to 1 subjet (i.e. the full ca_jet)
+  list<internal_recursive_softdrop::RSDHistoryElement*> active_branches;
+  active_branches.push_back(& (history[0]));
+
+  PseudoJet parent, piece1, piece2;
+  
+  while ((continue_grooming(n_depth)) && (active_branches.size())) {
+    // loop over all the branches and look for substructure
+    list<internal_recursive_softdrop::RSDHistoryElement*>::iterator hist_it=active_branches.begin();
+    while (hist_it!=active_branches.end()){
+      // get the element corresponding to the max dR
+      // and the associated PJ
+      internal_recursive_softdrop::RSDHistoryElement * elm = (*hist_it);
+      PseudoJet parent = cs_jets[cs_history[elm->current_in_ca_tree].jetp_index];
+
+      // we need to iterate this branch until we find some substructure
+      PseudoJet result_sd;
+      if (_dynamical_R0){
+        SoftDrop sd(_beta, _symmetry_cut, symmetry_measure(), sqrt(elm->theta_squared),
+                    mu_cut(), recursion_choice(), subtractor());
+        sd.set_reclustering(false);
+        sd.set_verbose_structure(has_verbose_structure());
+        result_sd = sd(parent);
+      } else {
+        result_sd = SoftDrop::result(parent);
+      }
+
+      // if we had an empty PJ, that means we ran into some problems.
+      // just return an empty PJ ourselves
+      if (result_sd == 0) return PseudoJet();
+      
+      // update the history element to reflect our iteration
+      elm->current_in_ca_tree = result_sd.cluster_hist_index();
+      
+      if (has_verbose_structure()){
+        elm->dropped_delta_R    = result_sd.structure_of<SoftDrop>().dropped_delta_R();
+        elm->dropped_symmetry   = result_sd.structure_of<SoftDrop>().dropped_symmetry();
+        elm->dropped_mu         = result_sd.structure_of<SoftDrop>().dropped_mu();
+      }
+
+      // if some substructure was found:
+      if (result_sd.structure_of<SoftDrop>().has_substructure()){
+        // update the history element to reflect our iteration
+        elm->child1_in_history  = history.size();
+        elm->child2_in_history  = history.size()+1;
+        elm->theta_squared      = result_sd.structure_of<SoftDrop>().delta_R();
+        elm->theta_squared     *= elm->theta_squared;
+        elm->symmetry           = result_sd.structure_of<SoftDrop>().symmetry();
+        elm->mu2                = result_sd.structure_of<SoftDrop>().mu();
+        elm->mu2               *= elm->mu2;
+             
+        // the next iteration will have to handle 2 new history
+        // elements (the R0squared argument here is unused)
+        result_sd.has_parents(piece1, piece2);
+        internal_recursive_softdrop::RSDHistoryElement elm1(piece1, this, _R0sqr);
+        history.push_back(elm1);
+        // insert it in the active branches if needed
+        if (elm1.theta_squared>0)
+          active_branches.insert(hist_it,&(history.back())); // insert just before
+
+        internal_recursive_softdrop::RSDHistoryElement elm2(piece2, this, _R0sqr);
+        history.push_back(elm2);
+        if ((!_hardest_branch_only) && (elm2.theta_squared>0)){
+          active_branches.insert(hist_it,&(history.back())); // insert just before
+        }
+      }
+      // otherwise we've just reached the end of the recursion the
+      // history information has been updated above
+      //
+      // we just need to make sure that we do not recurse into that
+      // element any longer
+
+      list<internal_recursive_softdrop::RSDHistoryElement*>::iterator current = hist_it;
+      ++hist_it;
+      active_branches.erase(current);
+    } // loop over branches at current depth
+    ++n_depth;
+  } // loop over depth
+
+  // now we have a bunch of history elements that we can use to build the final jet
+  vector<PseudoJet> mapped_to_history(history.size());
+  unsigned int history_index = history.size();
+  do {
+    --history_index;
+    const internal_recursive_softdrop::RSDHistoryElement & elm = history[history_index];
+
+    // two kinds of events: either just a final leave, poteitially  with grooming
+    // or a brandhing (also with potential grooming at the end)
+    if (elm.child1_in_history<0){
+      // this is a leaf, i.e. with no further sustructure
+      PseudoJet & subjet = mapped_to_history[history_index]
+        = cs_jets[cs_history[elm.current_in_ca_tree].jetp_index];
+
+      StructureType * structure = new StructureType(subjet);
+      if (has_verbose_structure()){
+        structure->set_verbose(true);
+      }
+      subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
+    } else {
+      PseudoJet & subjet = mapped_to_history[history_index]
+        = join(mapped_to_history[elm.child1_in_history], mapped_to_history[elm.child2_in_history]);
+      StructureType * structure = new StructureType(subjet, sqrt(elm.theta_squared), elm.symmetry, sqrt(elm.mu2));
+      if (has_verbose_structure()){
+        structure->set_verbose(true);
+        structure->set_dropped_delta_R (elm.dropped_delta_R);
+        structure->set_dropped_symmetry(elm.dropped_symmetry);
+        structure->set_dropped_mu      (elm.dropped_mu);
+      }
+      subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
+    }    
+  } while (history_index>0);
+
+  return mapped_to_history[0];
+}
+
+
+//========================================================================
+// implementation of the helpers
+//========================================================================
+
+// helper to get all the prongs in a jet that has been obtained using
+// RecursiveSoftDrop (instead of recursively parsing the 1->2
+// composite jet structure)
+vector<PseudoJet> recursive_soft_drop_prongs(const PseudoJet & rsd_jet){
+  // make sure that the jet has the appropriate RecursiveSoftDrop structure
+  if (!rsd_jet.has_structure_of<RecursiveSoftDrop>())
+    return vector<PseudoJet>();
+
+  // if this jet has no substructure, just return a 1-prong object
+  if (!rsd_jet.structure_of<RecursiveSoftDrop>().has_substructure())
+    return vector<PseudoJet>(1, rsd_jet);
+  
+  // otherwise fill a vector with all the prongs (no specific ordering)
+  vector<PseudoJet> prongs;
+
+  // parse the list of PseudoJet we still need to deal with
+  vector<PseudoJet> to_parse = rsd_jet.pieces();  // valid both for a C/A recombination step or a RSD join
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const PseudoJet &current = to_parse[i_parse];
+
+    if ((current.has_structure_of<RecursiveSymmetryCutBase>()) &&
+        (current.structure_of<RecursiveSymmetryCutBase>().has_substructure())){
+      // if this has some deeper substructure, add it to the list of
+      // things to further process
+      vector<PseudoJet> pieces = current.pieces();
+      assert(pieces.size() == 2);
+      to_parse[i_parse] = pieces[0];
+      to_parse.push_back(pieces[1]);
+    } else {
+      // no further substructure, just add this as a branch
+      prongs.push_back(current);
+      ++i_parse;
+    }
+  }
+   
+  return prongs;
+}
+    
+}
+
+FASTJET_END_NAMESPACE
Index: external/fastjet/contribs/RecursiveTools/RecursiveSoftDrop.hh
===================================================================
--- external/fastjet/contribs/RecursiveTools/RecursiveSoftDrop.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
+++ external/fastjet/contribs/RecursiveTools/RecursiveSoftDrop.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -0,0 +1,207 @@
+// $Id: RecursiveSoftDrop.hh 1082 2017-10-10 12:00:13Z gsoyez $
+//
+// Copyright (c) 2014-, Gavin P. Salam, Gregory Soyez, Jesse Thaler,
+// Kevin Zhou, Frederic Dreyer
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+
+#ifndef __RECURSIVESOFTDROP_HH__
+#define __RECURSIVESOFTDROP_HH__
+
+#include "Recluster.hh"
+#include "SoftDrop.hh"
+#include "fastjet/WrappedStructure.hh"
+
+#include <iostream>
+#include <queue>
+#include <vector>
+
+FASTJET_BEGIN_NAMESPACE
+
+namespace contrib{
+
+//------------------------------------------------------------------------
+/// \class RecursiveSoftDrop
+/// An implementation of the RecursiveSoftDrop.
+///
+/// Recursive Soft Drop will recursively groom a jet, removing
+/// particles that fail the criterion
+/// \f[
+///     z > z_{\rm cut} (\theta/R0)^\beta
+/// \f]
+/// until n subjets have been found.
+///
+/// Several variants are supported:
+///  - set_fixed_depth_mode() switches to fixed depth on all branches
+///    of the clustering tree
+///  - set_dynamical_R0() switches to dynamical R0 implementation of
+///    RSD
+///  - set_hardest_branch_only() switches to following only the
+///    hardest branch (e.g. for Iterated Soft Drop)
+///  - set_min_deltaR_square(val) sets a minimum angle considered for
+///    substructure (e.g. for Iterated Soft Drop)
+///
+/// Notes:
+///
+///  - Even though the calls to "set_tagging_mode()" or
+///    "set_grooming_mode(false)" are allowed, they should not be used
+///    with n=-1, and the default grooming_mode has to remain
+///    untouched (except for beta<0 and finite n).
+///
+//----------------------------------------------------------------------
+class RecursiveSoftDrop : public SoftDrop {
+public:
+  /// Simplified constructor. This takes the value of the "beta"
+  /// parameter and the symmetry cut (applied by default on the
+  /// scalar_z variable, as for the mMDT). It also takes an optional
+  /// subtractor.
+  ///
+  /// n is the number of times we require the SoftDrop condition to be
+  /// satisfied. n=-1 means infinity, i.e. we recurse into the jet
+  /// until individual constituents
+  ///
+  /// If the (optional) pileup subtractor can be supplied, then see
+  /// also the documentation for the set_input_jet_is_subtracted() member
+  /// function.
+  ///
+  /// \param beta               the value of the beta parameter
+  /// \param symmetry_cut       the value of the cut on the symmetry measure
+  /// \param n                  the requested number of iterations
+  /// \param R0                 the angular distance normalisation [1 by default]
+  RecursiveSoftDrop(double beta,
+                    double symmetry_cut,
+                    int n = -1,
+                    double R0 = 1,
+                    const FunctionOfPseudoJet<PseudoJet> * subtractor = 0) : 
+    SoftDrop(beta, symmetry_cut, R0, subtractor), _n(n) { set_defaults(); }
+
+  /// Full constructor, which takes the following parameters:
+  ///
+  /// \param beta               the value of the beta parameter
+  /// \param symmetry_cut       the value of the cut on the symmetry measure
+  /// \param symmetry_measure   the choice of measure to use to estimate the symmetry
+  /// \param n                  the requested number of iterations
+  /// \param R0                 the angular distance normalisation [1 by default]
+  /// \param mu_cut             the maximal allowed value of mass drop variable mu = m_heavy/m_parent 
+  /// \param recursion_choice   the strategy used to decide which subjet to recurse into
+  /// \param subtractor         an optional pointer to a pileup subtractor (ignored if zero)
+  RecursiveSoftDrop(double           beta,
+                    double           symmetry_cut, 
+                    SymmetryMeasure  symmetry_measure,
+                    int              n = -1,
+                    double           R0 = 1.0,
+                    double           mu_cut = std::numeric_limits<double>::infinity(), 
+                    RecursionChoice  recursion_choice = larger_pt,
+                    const FunctionOfPseudoJet<PseudoJet> * subtractor = 0) : 
+    SoftDrop(beta, symmetry_cut, symmetry_measure, R0, mu_cut, recursion_choice, subtractor),
+    _n(n) { set_defaults(); }
+
+  /// default destructor
+  virtual ~RecursiveSoftDrop(){}
+
+  //----------------------------------------------------------------------
+  // access to class info
+  int n() const { return _n; }
+  
+  //----------------------------------------------------------------------
+  // on top of the tweaks that we inherit from SoftDrop (via
+  // RecursiveSymmetryBase):
+  //  - set_verbose_structure()
+  //  - set_subtractor()
+  //  - set_input_jet_is_subtracted()
+  // we provide several other knobs, given below
+
+  /// initialise all the flags below to their default value
+  void set_defaults();
+  
+  /// switch to using the "same depth" variant where instead of
+  /// recursing from large to small angles and requiring n SD
+  /// conditions to be met (our default), we recurse simultaneously in
+  /// all the branches found during the previous iteration, up to a
+  /// maximum depth of n.
+  /// default: false
+  void set_fixed_depth_mode(bool value=true) { _fixed_depth = value; }
+  bool fixed_depth_mode() const { return _fixed_depth; }
+  
+  /// switch to using a dynamical R0 (used for the normalisation of
+  /// the symmetry measure) set by the last deltaR at which some
+  /// substructure was found.
+  /// default: false
+  void set_dynamical_R0(bool value=true) { _dynamical_R0 = value; }
+  bool use_dynamical_R0() const { return _dynamical_R0; }
+
+  /// when finding some substructure, only follow the hardest branch
+  /// for the recursion
+  /// default: false (i.e. recurse in both branches)
+  void set_hardest_branch_only(bool value=true) { _hardest_branch_only = value; }
+  bool use_hardest_branch_only() const { return _hardest_branch_only; }
+
+  /// set the minimum angle (squared) that we should consider for
+  /// substructure
+  /// default: -1.0 (i.e. no minimum)
+  void set_min_deltaR_squared(double value=-1.0) { _min_dR2 = value; }
+  double   min_deltaR_squared() const { return _min_dR2; }
+
+  /// description of the tool
+  virtual std::string description() const;
+
+  //----------------------------------------------------------------------
+  /// action on a single jet with RecursiveSoftDrop.
+  ///
+  /// uses "result_fixed_tags" by default (i.e. recurse from R0 to
+  /// smaller angles until n SD conditions have been met), or
+  /// "result_fixed_depth" where each of the previous SD branches are
+  /// recirsed into down to a depth of n.
+  virtual PseudoJet result(const PseudoJet &jet) const;
+
+  /// this routine applies the Soft Drop criterion recursively on the
+  /// CA tree until we find n subjets (or until it converges), and
+  /// adds them together into a groomed PseudoJet
+  PseudoJet result_fixed_tags(const PseudoJet &jet) const;
+
+  /// this routine applies the Soft Drop criterion recursively on the
+  /// CA tree, recursing into all the branches found during the previous iteration
+  /// until n layers have been found (or until it converges)
+  PseudoJet result_fixed_depth(const PseudoJet &jet) const;
+    
+protected:  
+  /// return false if we reached desired layer of grooming _n
+  bool continue_grooming(int current_n) const {
+    return ((_n < 0) or (current_n < _n));
+  }
+  
+private:
+  int    _n;            ///< the value of n
+
+  // behaviour tweaks
+  bool _fixed_depth;         ///< look in parallel into each all branches until depth n
+  bool _dynamical_R0;        ///< when true, use the last deltaR with substructure as D0
+  bool _hardest_branch_only; ///< recurse only in the hardest branch
+                             ///  when substructure is found
+  double _min_dR2;           ///< the min allowed angle to search for substructure
+};
+
+// helper to get the (linear) list of prongs inside a jet resulting
+// from RecursiveSoftDrop. This would avoid having amnually to go
+// through the successive pairwise compositeness
+std::vector<PseudoJet> recursive_soft_drop_prongs(const PseudoJet & rsd_jet);
+
+}
+
+FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
+#endif // __RECURSIVESOFTDROP_HH__
Index: external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.cc
===================================================================
--- external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.cc	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.cc	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,3 +1,3 @@
-// $Id: RecursiveSymmetryCutBase.cc 700 2014-07-07 12:50:05Z gsoyez $
+// $Id: RecursiveSymmetryCutBase.cc 1080 2017-09-28 07:51:37Z gsoyez $
 //
 // Copyright (c) 2014-, Gavin P. Salam, Gregory Soyez, Jesse Thaler
@@ -23,4 +23,5 @@
 #include "fastjet/JetDefinition.hh"
 #include "fastjet/ClusterSequenceAreaBase.hh"
+#include <cassert>
 #include <algorithm> 
 #include <cstdlib> 
@@ -42,25 +43,17 @@
 PseudoJet RecursiveSymmetryCutBase::result(const PseudoJet & jet) const {
   // construct the input jet (by default, recluster with C/A)
-
   if (! jet.has_constituents()){
     throw Error("RecursiveSymmetryCutBase can only be applied to jets with constituents");
   }
 
-  PseudoJet j = 
-    _do_reclustering 
-      ? _recluster ? (*_recluster)(jet)
-                   : Recluster(cambridge_algorithm, JetDefinition::max_allowable_R)(jet)
-      : jet;
-    
-  // issue a warning if the jet is not obtained through a C/A
-  // clustering
-  // if ((! j.has_associated_cluster_sequence()) ||
-  //     (j.validated_cs()->jet_def().jet_algorithm() != cambridge_algorithm))
-  //  _nonca_warning.warn("RecursiveSymmetryCutBase is designed to be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms at your own risk.");
-
+  PseudoJet j = _recluster_if_needed(jet);
+
+  // sanity check: the jet must have a valid CS
   if (! j.has_valid_cluster_sequence()){
     throw Error("RecursiveSymmetryCutBase can only be applied to jets associated to a (valid) cluster sequence");
   }
 
+  // check that area information is there in case we have a subtractor
+  // GS: do we really need this since subtraction may not require areas?
   if (_subtractor) {
     const ClusterSequenceAreaBase * csab = 
@@ -75,6 +68,4 @@
     subjet = (*_subtractor)(subjet);
   }
-
-  bool use_mu_cut = (_mu_cut != numeric_limits<double>::infinity());
 
   // variables for tracking what will happen
@@ -86,90 +77,41 @@
   std::vector<double> dropped_symmetry;
   std::vector<double> dropped_mu;
+
+  double sym, mu2;
   
   // now recurse into the jet's structure
-  while (subjet.has_parents(piece1, piece2)) {
-    
-    // first sanity check: 
-    // - zero or negative pts are not allowed for the input subjet
-    // - zero or negative masses are not allowed for configurations
-    //   in which the mass will effectively appear in a denominator
-    //   (The masses will be checked later)
-    if (subjet.pt2() <= 0) return PseudoJet();
-
-    if (_subtractor) {
-      piece1 = (*_subtractor)(piece1);
-      piece2 = (*_subtractor)(piece2);
-    }
-    
-    // determine the symmetry parameter
-    double sym;
-
-    if        (_symmetry_measure == y) {
-      // the original d_{ij}/m^2 choice from MDT
-      // first make sure the mass is sensible
-      if (subjet.m2() <= 0) {
-        _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot calculate y, because (sub)jet mass is negative; bailing out");
-        return _result_no_substructure(subjet); //TBC: do we return the hardest parent? A NULL PseudoJet?
-      }
-      sym = piece1.kt_distance(piece2) / subjet.m2();
-
-    } else if (_symmetry_measure == vector_z) {
-      // min(pt1, pt2)/(pt), where the denominator is a vector sum
-      // of the two subjets
-      sym = min(piece1.pt(), piece2.pt()) / subjet.pt();
-
-    } else if (_symmetry_measure == scalar_z) {
-      // min(pt1, pt2)/(pt1+pt2), where the denominator is a scalar sum
-      // of the two subjets
-      double pt1 = piece1.pt();
-      double pt2 = piece2.pt();
-      // make sure denominator is non-zero
-      sym = pt1 + pt2;
-      if (sym == 0) return PseudoJet(); 
-      sym = min(pt1, pt2) / sym;
-
-    } else {
-      throw Error ("Unrecognized choice of symmetry_measure");
-    }
-
-    // determine the symmetry cut
-    // (This function is specified in the derived classes)
-    double this_symmetry_cut = symmetry_cut_fn(piece1, piece2);
-
-    // and make a first tagging decision based on symmetry cut
-    bool tagged = (sym > this_symmetry_cut);
-
-    // if tagged based on symmetry cut, then check the mu cut (if relevant)
-    // and update the tagging decision. Calculate mu^2 regardless, for cases
-    // of users not cutting on mu2, but still interested in its value.
-    double mu2 = max(piece1.m2(), piece2.m2())/subjet.m2();
-    if (tagged && use_mu_cut) {
-      // first a sanity check -- mu2 won't be sensible if the subjet mass 
-      // is negative, so we can't then trust the mu cut - bail out
-      if (subjet.m2() <= 0) {
-        _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot trust mu, because (sub)jet mass is negative; bailing out");
-        return PseudoJet();
-      }
-      if (mu2 > 1) _mu2_gt1_warning.warn("RecursiveSymmetryCutBase encountered mu^2 value > 1");
-      if (mu2 > pow(_mu_cut,2)) tagged = false;
-    }
-
-
-    // if we've tagged the splitting, return the jet with its substructure
-    if (tagged) {
-      // record relevant information
-      StructureType * structure = new StructureType(subjet);
-      structure->_symmetry = sym;
-      structure->_mu       = (mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2);
-      structure->_delta_R  = piece1.delta_R(piece2);
-      if (_verbose_structure) {
-        structure->_has_verbose = true;
-        structure->_dropped_symmetry = dropped_symmetry;
-        structure->_dropped_mu = dropped_mu;
-        structure->_dropped_delta_R = dropped_delta_R;
-      }
-      subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
-      return subjet;
-    }
+  RecursionStatus status;
+  while ((status=recurse_one_step(subjet, piece1, piece2, sym, mu2)) != recursion_success) {
+    // start with sanity checks:
+    if ((status == recursion_issue) || (status == recursion_no_parents)) {
+      // we should return piece1 by our convention for recurse_one_step
+      PseudoJet result;
+      if (status == recursion_issue){
+        result = piece1;
+        if (_verbose) cout << "reached end; returning null jet " << endl;
+      } else {
+        result = _result_no_substructure(piece1);
+        if (_verbose) cout << "no parents found; returning last PJ or empty jet" << endl;
+      }
+      
+      if (result != 0) {
+        // if in grooming mode, add dummy structure information
+        StructureType * structure = new StructureType(result);
+        // structure->_symmetry = 0.0;
+        // structure->_mu       = 0.0;
+        // structure->_delta_R  = 0.0;
+        if (_verbose_structure) { // still want to store verbose information about dropped branches
+          structure->_has_verbose = true;
+          structure->_dropped_symmetry = dropped_symmetry;
+          structure->_dropped_mu = dropped_mu;
+          structure->_dropped_delta_R = dropped_delta_R;
+        }
+        result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
+      }
+  
+      return result;
+    }
+
+    assert(status == recursion_dropped);
     
     // if desired, store information about dropped branches before recursing
@@ -179,48 +121,151 @@
       dropped_mu.push_back((mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2));
     }
-    
-    // otherwise continue unclustering, allowing for the different
-    // ways of choosing which parent to look into
-    int choice;
-    if        (_recursion_choice == larger_mt) {
-      choice = piece1.mt2() > piece2.mt2() ? 1 : 2;
-
-    } else if (_recursion_choice == larger_pt) {
-      choice = piece1.pt2() > piece2.pt2() ? 1 : 2;
-
-    } else if (_recursion_choice == larger_m)  {
-      choice = piece1.m2()  > piece2.m2()  ? 1 : 2;
-
-    } else {
-      throw Error ("Unrecognized value for recursion_choice");
-    }    
-    if (_verbose) cout << "choice is " << choice << endl;;
-    subjet = (choice == 1) ? piece1 : piece2;
-  } // (subjet.has_parents(...))
-
-  if (_verbose) cout << "reached end; returning null jet " << endl;
-  
-  // decide on tagging versus grooming mode here
-  PseudoJet result = _result_no_substructure(subjet);
-  
-  if (result != 0) {
-    // if in grooming mode, add dummy structure information
-    StructureType * structure = new StructureType(result);
-    structure->_symmetry = 0.0;
-    structure->_mu       = 0.0;
-    structure->_delta_R  = 0.0;
-    if (_verbose_structure) { // still want to store verbose information about dropped branches
-      structure->_has_verbose = true;
-      structure->_dropped_symmetry = dropped_symmetry;
-      structure->_dropped_mu = dropped_mu;
-      structure->_dropped_delta_R = dropped_delta_R;
-    }
-    result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
-  }
-  
-  return result;
-}
-
-
+
+    subjet = piece1;
+  }
+  
+
+  // we've tagged the splitting, return the jet with its substructure
+  StructureType * structure = new StructureType(subjet);
+  structure->_symmetry = sym;
+  structure->_mu       = (mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2);
+  structure->_delta_R  = sqrt(squared_geometric_distance(piece1, piece2));
+  if (_verbose_structure) {
+    structure->_has_verbose = true;
+    structure->_dropped_symmetry = dropped_symmetry;
+    structure->_dropped_mu = dropped_mu;
+    structure->_dropped_delta_R = dropped_delta_R;
+  }
+  subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
+  return subjet;
+}
+
+
+  
+//----------------------------------------------------------------------
+// the method below is the one actually performing one step of the
+// recursion.
+//
+// It returns a status code (defined above)
+//
+// In case of success, all the information is filled
+// In case of "no parents", piee1 is the same subjet
+// In case of trouble, piece2 will be a 0 PJ and piece1 is the PJ we
+//   should return (either 0 itself if the issue was critical, or
+//   non-wero in case of a minor issue just causing the recursion to
+//   stop)
+RecursiveSymmetryCutBase::RecursionStatus
+      RecursiveSymmetryCutBase::recurse_one_step(const PseudoJet & subjet,
+                                                 PseudoJet &piece1, PseudoJet &piece2,
+                                                 double &sym, double &mu2,
+                                                 void *extra_parameters) const {
+  if (!subjet.has_parents(piece1, piece2)){
+    piece1 = subjet;    
+    piece2 = PseudoJet();
+    return recursion_no_parents;
+  }
+
+  // first sanity check: 
+  // - zero or negative pts are not allowed for the input subjet
+  // - zero or negative masses are not allowed for configurations
+  //   in which the mass will effectively appear in a denominator
+  //   (The masses will be checked later)
+  if (subjet.pt2() <= 0){ // this is a critical problem, return an empty PJ
+    piece1 = piece2 = PseudoJet();
+    return recursion_issue;
+  }
+
+  if (_subtractor) {
+    piece1 = (*_subtractor)(piece1);
+    piece2 = (*_subtractor)(piece2);
+  }
+    
+  // determine the symmetry parameter
+  if        (_symmetry_measure == y) {
+    // the original d_{ij}/m^2 choice from MDT
+    // first make sure the mass is sensible
+    if (subjet.m2() <= 0) { 
+      _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot calculate y, because (sub)jet mass is negative; bailing out");
+      // since rounding errors can give -ve masses, be a it more
+      // tolerant and consider that no substructure has been found
+      piece1 = _result_no_substructure(subjet);
+      piece2 = PseudoJet();
+      return recursion_issue;
+    }
+    sym = piece1.kt_distance(piece2) / subjet.m2();
+    
+  } else if (_symmetry_measure == vector_z) {
+    // min(pt1, pt2)/(pt), where the denominator is a vector sum
+    // of the two subjets
+    sym = min(piece1.pt(), piece2.pt()) / subjet.pt();    
+  } else if (_symmetry_measure == scalar_z) {
+    // min(pt1, pt2)/(pt1+pt2), where the denominator is a scalar sum
+    // of the two subjets
+    double pt1 = piece1.pt();
+    double pt2 = piece2.pt();
+    // make sure denominator is non-zero
+    sym = pt1 + pt2;
+    if (sym == 0){ // this is a critical problem, return an empty PJ
+      piece1 = piece2 = PseudoJet();
+      return recursion_issue;
+    }
+    sym = min(pt1, pt2) / sym;
+  } else if ((_symmetry_measure == theta_E) || (_symmetry_measure == cos_theta_E)){
+    // min(E1, E2)/(E1+E2)
+    double E1 = piece1.E();
+    double E2 = piece2.E();
+    // make sure denominator is non-zero
+    sym = E1 + E2;
+    if (sym == 0){ // this is a critical problem, return an empty PJ
+      piece1 = piece2 = PseudoJet();
+      return recursion_issue;
+    }
+    sym = min(E1, E2) / sym;
+  } else {
+    throw Error ("Unrecognized choice of symmetry_measure");
+  }
+  
+  // determine the symmetry cut
+  // (This function is specified in the derived classes)
+  double this_symmetry_cut = symmetry_cut_fn(piece1, piece2, extra_parameters);
+  
+  // and make a first tagging decision based on symmetry cut
+  bool tagged = (sym > this_symmetry_cut);
+
+  // if tagged based on symmetry cut, then check the mu cut (if relevant)
+  // and update the tagging decision. Calculate mu^2 regardless, for cases
+  // of users not cutting on mu2, but still interested in its value.
+  bool use_mu_cut = (_mu_cut != numeric_limits<double>::infinity());
+  mu2 = max(piece1.m2(), piece2.m2())/subjet.m2();
+  if (tagged && use_mu_cut) {
+    // first a sanity check -- mu2 won't be sensible if the subjet mass 
+    // is negative, so we can't then trust the mu cut - bail out
+    if (subjet.m2() <= 0) {
+      _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot trust mu, because (sub)jet mass is negative; bailing out");
+      piece1 = piece2 = PseudoJet();
+      return recursion_issue;
+    }
+    if (mu2 > 1) _mu2_gt1_warning.warn("RecursiveSymmetryCutBase encountered mu^2 value > 1");
+    if (mu2 > pow(_mu_cut,2)) tagged = false;
+  }
+
+  // we'll continue unclustering, allowing for the different
+  // ways of choosing which parent to look into
+  if        (_recursion_choice == larger_pt) {
+    if (piece1.pt2() < piece2.pt2()) std::swap(piece1, piece2);    
+  } else if (_recursion_choice == larger_mt) {
+    if (piece1.mt2() < piece2.mt2()) std::swap(piece1, piece2);    
+  } else if (_recursion_choice == larger_m)  {
+    if (piece1.m2() < piece2.m2()) std::swap(piece1, piece2);    
+  } else if (_recursion_choice == larger_E)  {
+    if (piece1.E()  < piece2.E())  std::swap(piece1, piece2);    
+  } else {
+    throw Error ("Unrecognized value for recursion_choice");
+  }    
+
+  return tagged ? recursion_success : recursion_dropped;
+}
+
+  
 //----------------------------------------------------------------------
 string RecursiveSymmetryCutBase::description() const {
@@ -235,4 +280,8 @@
   case vector_z:
     ostr << "vector_z"; break;
+  case theta_E:
+    ostr << "theta_E"; break;
+  case cos_theta_E:
+    ostr << "cos_theta_E"; break;
   default:
     cerr << "failed to interpret symmetry_measure" << endl; exit(-1);
@@ -254,4 +303,6 @@
   case larger_m:
     ostr << "mass"; break;
+  case larger_E:
+    ostr << "energy"; break;
   default:
     cerr << "failed to interpret recursion_choice" << endl; exit(-1);
@@ -259,11 +310,49 @@
 
   if (_subtractor) {
-    ostr << " and subtractor: " << _subtractor->description();
+    ostr << ", subtractor: " << _subtractor->description();
     if (_input_jet_is_subtracted) {ostr << " (input jet is assumed already subtracted)";}
   }
+
+  if (_recluster) {
+    ostr << " and reclustering using " << _recluster->description();
+  }
+  
   return ostr.str();
 }
 
+//----------------------------------------------------------------------
+// helper for handling the reclustering
+PseudoJet RecursiveSymmetryCutBase::_recluster_if_needed(const PseudoJet &jet) const{
+  if (! _do_reclustering) return jet;
+  if (_recluster) return (*_recluster)(jet);
+  if (is_ee()){
+#if FASTJET_VERSION_NUMBER >= 30100
+    return Recluster(JetDefinition(ee_genkt_algorithm, JetDefinition::max_allowable_R, 0.0), true)(jet);
+#else
+    return Recluster(JetDefinition(ee_genkt_algorithm, JetDefinition::max_allowable_R, 0.0))(jet);
+#endif
+  }
+
+  return Recluster(cambridge_algorithm, JetDefinition::max_allowable_R)(jet);
+}  
+  
+//----------------------------------------------------------------------
 // decide what to return when no substructure has been found
+double RecursiveSymmetryCutBase::squared_geometric_distance(const PseudoJet &j1,
+                                                            const PseudoJet &j2) const{
+  if (_symmetry_measure == theta_E){
+    double dot_3d = j1.px()*j2.px() + j1.py()*j2.py() + j1.pz()*j2.pz();
+    double cos_theta = max(-1.0,min(1.0, dot_3d/sqrt(j1.modp2()*j2.modp2())));
+    double theta = acos(cos_theta);
+    return theta*theta;
+  } else if (_symmetry_measure == cos_theta_E){
+    double dot_3d = j1.px()*j2.px() + j1.py()*j2.py() + j1.pz()*j2.pz();
+    return max(0.0, 2*(1-dot_3d/sqrt(j1.modp2()*j2.modp2())));
+  }
+
+  return j1.squared_distance(j2);
+}
+
+//----------------------------------------------------------------------
 PseudoJet RecursiveSymmetryCutBase::_result_no_substructure(const PseudoJet &last_parent) const{
   if (_grooming_mode){
@@ -277,4 +366,276 @@
 
 
+//========================================================================
+// implementation of the details of the structure
+  
+// the number of dropped subjets
+int RecursiveSymmetryCutBase::StructureType::dropped_count(bool global) const {
+  check_verbose("dropped_count()");
+
+  // if this jet has no substructure, just return an empty vector
+  if (!has_substructure()) return _dropped_delta_R.size();
+
+  // deal with the non-global case
+  if (!global) return _dropped_delta_R.size();
+
+  // for the global case, we've unfolded the recursion (likely more
+  // efficient as it requires less copying)
+  unsigned int count = 0;
+  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    count += current->_dropped_delta_R.size();
+
+    // check if we need to recurse deeper in the substructure
+    //
+    // we can have 2 situations here for the underlying structure (the
+    // one we've wrapped around):
+    //  - it's of the clustering type
+    //  - it's a composite jet
+    // only in the 2nd case do we have to recurse deeper
+    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(current->_structure.get());
+    if (css == 0){ ++i_parse; continue; }
+    
+    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
+    assert(prongs.size() == 2);
+    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
+      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
+        RecursiveSymmetryCutBase::StructureType* prong_structure
+          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
+        if (prong_structure->has_substructure())
+          to_parse.push_back(prong_structure);
+      }
+    }
+
+    ++i_parse;
+  }
+  return count;
+}
+
+// the delta_R of all the dropped subjets
+vector<double> RecursiveSymmetryCutBase::StructureType::dropped_delta_R(bool global) const {
+  check_verbose("dropped_delta_R()");
+
+  // if this jet has no substructure, just return an empty vector
+  if (!has_substructure()) return vector<double>();
+
+  // deal with the non-global case
+  if (!global) return _dropped_delta_R;
+
+  // for the global case, we've unfolded the recursion (likely more
+  // efficient as it requires less copying)
+  vector<double> all_dropped;
+  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    all_dropped.insert(all_dropped.end(), current->_dropped_delta_R.begin(), current->_dropped_delta_R.end());
+
+    // check if we need to recurse deeper in the substructure
+    //
+    // we can have 2 situations here for the underlying structure (the
+    // one we've wrapped around):
+    //  - it's of the clustering type
+    //  - it's a composite jet
+    // only in the 2nd case do we have to recurse deeper
+    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(current->_structure.get());
+    if (css == 0){ ++i_parse; continue; }
+    
+    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
+    assert(prongs.size() == 2);
+    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
+      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
+        RecursiveSymmetryCutBase::StructureType* prong_structure
+          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
+        if (prong_structure->has_substructure())
+          to_parse.push_back(prong_structure);
+      }
+    }
+
+    ++i_parse;
+  }
+  return all_dropped;
+}
+
+// the symmetry of all the dropped subjets
+vector<double> RecursiveSymmetryCutBase::StructureType::dropped_symmetry(bool global) const {
+  check_verbose("dropped_symmetry()");
+
+  // if this jet has no substructure, just return an empty vector
+  if (!has_substructure()) return vector<double>();
+
+  // deal with the non-global case
+  if (!global) return _dropped_symmetry;
+
+  // for the global case, we've unfolded the recursion (likely more
+  // efficient as it requires less copying)
+  vector<double> all_dropped;
+  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    all_dropped.insert(all_dropped.end(), current->_dropped_symmetry.begin(), current->_dropped_symmetry.end());
+
+    // check if we need to recurse deeper in the substructure
+    //
+    // we can have 2 situations here for the underlying structure (the
+    // one we've wrapped around):
+    //  - it's of the clustering type
+    //  - it's a composite jet
+    // only in the 2nd case do we have to recurse deeper
+    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(current->_structure.get());
+    if (css == 0){ ++i_parse; continue; }
+    
+    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
+    assert(prongs.size() == 2);
+    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
+      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
+        RecursiveSymmetryCutBase::StructureType* prong_structure
+          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
+        if (prong_structure->has_substructure())
+          to_parse.push_back(prong_structure);
+      }
+    }
+
+    ++i_parse;
+  }
+  return all_dropped;
+}
+
+// the mu of all the dropped subjets
+vector<double> RecursiveSymmetryCutBase::StructureType::dropped_mu(bool global) const {
+  check_verbose("dropped_mu()");
+
+  // if this jet has no substructure, just return an empty vector
+  if (!has_substructure()) return vector<double>();
+
+  // deal with the non-global case
+  if (!global) return _dropped_mu;
+
+  // for the global case, we've unfolded the recursion (likely more
+  // efficient as it requires less copying)
+  vector<double> all_dropped;
+  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    all_dropped.insert(all_dropped.end(), current->_dropped_mu.begin(), current->_dropped_mu.end());
+
+    // check if we need to recurse deeper in the substructure
+    //
+    // we can have 2 situations here for the underlying structure (the
+    // one we've wrapped around):
+    //  - it's of the clustering type
+    //  - it's a composite jet
+    // only in the 2nd case do we have to recurse deeper
+    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(current->_structure.get());
+    if (css == 0){ ++i_parse; continue; }
+    
+    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
+    assert(prongs.size() == 2);
+    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
+      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
+        RecursiveSymmetryCutBase::StructureType* prong_structure
+          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
+        if (prong_structure->has_substructure())
+          to_parse.push_back(prong_structure);
+      }
+    }
+
+    ++i_parse;
+  }
+  return all_dropped;
+}
+
+// the maximum of the symmetry over the dropped subjets
+double RecursiveSymmetryCutBase::StructureType::max_dropped_symmetry(bool global) const {
+  check_verbose("max_dropped_symmetry()");
+
+  // if there is no substructure, just exit
+  if (!has_substructure()){ return 0.0; }
+
+  // local value of the max_dropped_symmetry
+  double local_max = (_dropped_symmetry.size() == 0)
+    ? 0.0 : *max_element(_dropped_symmetry.begin(),_dropped_symmetry.end());
+
+  // recurse down the structure if instructed to do so
+  if (global){
+    // we can have 2 situations here for the underlying structure (the
+    // one we've wrapped around):
+    //  - it's of the clustering type
+    //  - it's a composite jet
+    // only in the 2nd case do we have to recurse deeper
+    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(_structure.get());
+    if (css == 0) return local_max;
+    
+    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
+    assert(prongs.size() == 2);
+    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
+      // check if the prong has further substructure
+      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
+        RecursiveSymmetryCutBase::StructureType* prong_structure
+          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
+        local_max = max(local_max, prong_structure->max_dropped_symmetry(true));
+      }
+    }
+  }
+
+  return local_max;
+}
+
+//------------------------------------------------------------------------
+// helper class to sort by decreasing thetag
+class SortRecursiveSoftDropStructureZgThetagPair{
+public:
+  bool operator()(const pair<double, double> &p1, const pair<double, double> &p2) const{
+    return p1.second > p2.second;
+  }
+};
+//------------------------------------------------------------------------
+
+// the (zg,thetag) pairs of all the splitting that were found and passed the SD condition
+vector<pair<double,double> > RecursiveSymmetryCutBase::StructureType::sorted_zg_and_thetag() const {
+  //check_verbose("sorted_zg_and_thetag()");
+
+  // if this jet has no substructure, just return an empty vector
+  if (!has_substructure()) return vector<pair<double,double> >();
+  
+  // otherwise fill a vector with all the prongs (no specific ordering)
+  vector<pair<double,double> > all;
+  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    all.push_back(pair<double,double>(current->_symmetry, current->_delta_R));
+
+    vector<PseudoJet> prongs = current->pieces(PseudoJet());
+    assert(prongs.size() == 2);
+    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
+      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
+        RecursiveSymmetryCutBase::StructureType* prong_structure
+          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
+        if (prong_structure->has_substructure())
+          to_parse.push_back(prong_structure);
+      }
+    }
+
+    ++i_parse;
+  }
+
+  sort(all.begin(), all.end(), SortRecursiveSoftDropStructureZgThetagPair());
+  return all;
+}
+
 } // namespace contrib
 
Index: external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.hh
===================================================================
--- external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.hh	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,3 +1,3 @@
-// $Id: RecursiveSymmetryCutBase.hh 700 2014-07-07 12:50:05Z gsoyez $
+// $Id: RecursiveSymmetryCutBase.hh 1074 2017-09-18 15:15:20Z gsoyez $
 //
 // Copyright (c) 2014-, Gavin P. Salam, Gregory Soyez, Jesse Thaler
@@ -24,21 +24,29 @@
 
 #include <limits>
+#include <cassert>
 #include <fastjet/internal/base.hh>
 #include "fastjet/tools/Transformer.hh"
 #include "fastjet/WrappedStructure.hh"
-
+#include "fastjet/CompositeJetStructure.hh"
+
+#include "fastjet/config.h"
+
+// we'll use the native FJ class for reculstering if available
+#if FASTJET_VERSION_NUMBER >= 30100
+#include "fastjet/tools/Recluster.hh"
+#else
 #include "Recluster.hh"
+#endif
 
 /** \mainpage RecursiveTools contrib 
 
-    The aims RecursiveTools contrib to provide a set of tools for
-    recursive investigation of the substructure of jets.
-
-    Currently it includes:
+    The RecursiveTools contrib provides a set of tools for
+    recursive investigation jet substructure. Currently it includes:
     - fastjet::contrib::ModifiedMassDropTagger
     - fastjet::contrib::SoftDrop
-    - the two above classes derive from fastjet::contrib::RecursiveSymmetryCutBase
-    - fastjet::contrib::Recluster provides a reclustering transformer
-    - example*.cc provides a usage examples
+    - fastjet::contrib::RecursiveSymmetryCutBase (from which the above two classes derive)
+    - fastjet::contrib::IteratedSoftDropSymmetryFactors (defines ISD procedure)
+    - fastjet::contrib::IteratedSoftDropMultiplicity (defines a useful observable using ISD)  
+    - example*.cc provides usage examples
  */
 
@@ -78,8 +86,10 @@
 
   /// an enum of the different (a)symmetry measures that can be used
-  enum SymmetryMeasure{scalar_z, ///< \f$ \min(p_{ti}, p_{tj})/(p_{ti} + p_{tj}) \f$
-                       vector_z, ///< \f$ \min(p_{ti}, p_{tj})/p_{t(i+j)} \f$
-                       y         ///  \f$ \min(p_{ti}^2,p_{tj}^2) \Delta R_{ij}^2 / m_{ij}^2 \f$
-                       
+  enum SymmetryMeasure{scalar_z,   ///< \f$ \min(p_{ti}, p_{tj})/(p_{ti} + p_{tj}) \f$
+                       vector_z,   ///< \f$ \min(p_{ti}, p_{tj})/p_{t(i+j)} \f$
+                       y,          ///< \f$ \min(p_{ti}^2,p_{tj}^2) \Delta R_{ij}^2 / m_{ij}^2 \f$
+                       theta_E,    ///< \f$ \min(E_i,E_j)/(E_i+E_j) \f$ with 3d angle (ee collisions)
+                       cos_theta_E ///< \f$ \min(E_i,E_j)/(E_i+E_j) \f$ with
+                                   ///  \f$ \sqrt{2[1-cos(theta)]}\f$ for angles (ee collisions)
   };
 
@@ -87,5 +97,6 @@
   enum RecursionChoice{larger_pt, ///< choose the subjet with larger \f$ p_t \f$
                        larger_mt, ///< choose the subjet with larger \f$ m_t \equiv (m^2+p_t^2)^{\frac12}] \f$
-                       larger_m   ///  choose the subjet with larger mass (deprecated)
+                       larger_m,  ///< choose the subjet with larger mass (deprecated)
+                       larger_E   ///< choose the subjet with larger energy (meant for ee collisions)
   }; 
 
@@ -119,4 +130,10 @@
   virtual ~RecursiveSymmetryCutBase(){}
 
+  // access to class info
+  //----------------------------------------------------------------------
+  SymmetryMeasure symmetry_measure() const { return _symmetry_measure; }
+  double mu_cut() const { return _mu_cut; }
+  RecursionChoice recursion_choice() const { return _recursion_choice; }
+
   // internal subtraction configuration
   //----------------------------------------------------------------------
@@ -147,5 +164,5 @@
   //----------------------------------------------------------------------
 
-  /// configure the reclustering prior to the SoftDrop de-clustering
+  /// configure the reclustering prior to the recursive de-clustering
   ///  \param do_reclustering   recluster the jet or not?
   ///  \param recluster         how to recluster the jet 
@@ -174,4 +191,5 @@
   /// in particular values of symmetry, delta_R, and mu of dropped branches
   void set_verbose_structure(bool enable=true) { _verbose_structure = enable; }
+  bool has_verbose_structure() const { return _verbose_structure; }
   
   
@@ -185,10 +203,13 @@
   /// subtracted jet.
   virtual PseudoJet result(const PseudoJet & j) const;
-
+  
   /// description of the tool
   virtual std::string description() const;
 
-  /// the type of the associated structure
-  //typedef RecursiveSymmetryCutBaseStructure StructureType;  
+  /// returns the gepometrical distance between the two particles
+  /// depending on the symmetry measure used
+  double squared_geometric_distance(const PseudoJet &j1,
+                                    const PseudoJet &j2) const;
+  
 
   class StructureType;
@@ -203,7 +224,47 @@
   /// apply for a given pair of subjets p1 and p2
   virtual double symmetry_cut_fn(const PseudoJet & /* p1 */, 
-                                 const PseudoJet & /* p2 */) const = 0;
+                                 const PseudoJet & /* p2 */,
+                                 void *extra_parameters = 0) const = 0;
   /// the associated dwescription
   virtual std::string symmetry_cut_description() const = 0;
+
+  //----------------------------------------------------------------------
+  /// this defines status codes when checking for substructure
+  enum RecursionStatus{
+    recursion_success=0,   //< found some substructure
+    recursion_dropped,     //< dropped softest prong; recursion continues
+    recursion_no_parents,  //< down to constituents; bottom of recursion
+    recursion_issue        //< something went wrong; recursion stops
+  };
+  
+  //----------------------------------------------------------------------
+  /// the method below is the one actually performing one step of the
+  /// recursion.
+  ///
+  /// It returns a status code (defined above)
+  ///
+  /// In case of success, all the information is filled
+  /// In case of "no parents", piee1 is the same subjet
+  /// In case of trouble, piece2 will be a 0 PJ and piece1 is the PJ we
+  ///   should return (either 0 itself if the issue was critical, or
+  ///   non-wero in case of a minor issue just causing the recursion to
+  ///   stop)
+  ///
+  /// The extra_parameter argument allows one to pass extra agruments
+  /// to the symmetry condition
+  RecursionStatus recurse_one_step(const PseudoJet & subjet,
+                                   PseudoJet &piece1, PseudoJet &piece2,
+                                   double &sym, double &mu2,
+                                   void *extra_parameters = 0) const;
+
+  //----------------------------------------------------------------------
+  /// helper for handling the reclustering
+  PseudoJet _recluster_if_needed(const PseudoJet &jet) const;
+  
+  //----------------------------------------------------------------------
+  // helpers for selecting between ee and pp cases
+  bool is_ee() const{
+    return ((_symmetry_measure==theta_E) || (_symmetry_measure==cos_theta_E));
+  }
 
 private:
@@ -238,51 +299,68 @@
 public:
   StructureType(const PseudoJet & j) :
-    WrappedStructure(j.structure_shared_ptr()),
+    WrappedStructure(j.structure_shared_ptr()), _delta_R(-1.0), _symmetry(-1.0), _mu(-1.0),
+    _is_composite(false), _has_verbose(false) // by default, do not store verbose structure
+  {}
+
+  /// construct a structure with
+  ///  - basic info inherited from the reference jet "j"
+  ///  - a given deltaR       for substructure
+  ///  - a given symmetry     for substructure
+  ///  - a given mu parameter for substructure
+  /// If j is a "copmposite jet", it means that it has further
+  /// substructure to potentially recurse  into
+  StructureType(const PseudoJet & j, double delta_R_in, double symmetry_in, double mu_in=-1.0) :
+    WrappedStructure(j.structure_shared_ptr()), _delta_R(delta_R_in), _symmetry(symmetry_in), _mu(mu_in),
+    _is_composite(dynamic_cast<const CompositeJetStructure*>(j.structure_ptr()) != NULL),
     _has_verbose(false) // by default, do not store verbose structure
   {}
   
   // information about kept branch
-  double delta_R()  const {return _delta_R;};
-  double symmetry() const {return _symmetry;};
-  double mu()       const {return _mu;};
+  double delta_R()  const {return _delta_R;}
+  double thetag()   const {return _delta_R;}  // alternative name
+  double symmetry() const {return _symmetry;}
+  double zg()       const {return _symmetry;} // alternative name
+  double mu()       const {return _mu;}
   
   // additional verbose information about dropped branches
   bool has_verbose() const { return _has_verbose;}
-  
-  // number of dropped branches
-  int dropped_count() const {
-    if (!_has_verbose) throw Error("RecursiveSymmetryCutBase::StructureType: Verbose structure must be turned on to get dropped_count() values.");
-    return _dropped_delta_R.size();
-  }
-  
-  // delta_R of dropped branches
-  std::vector<double> dropped_delta_R() const {
-    if (!_has_verbose) throw Error("RecursiveSymmetryCutBase::StructureType: Verbose structure must be turned on to get dropped_delta_R() values.");
-    return _dropped_delta_R;
-  }
-
-  // symmetry values of dropped branches
-  std::vector<double> dropped_symmetry() const {
-    if (!_has_verbose) throw Error("RecursiveSymmetryCutBase::StructureType: Verbose structure must be turned on to get dropped_symmetry() values.");
-    return _dropped_symmetry;
-  }
-
-  // mass drop values of dropped branches
-  std::vector<double> dropped_mu() const {
-    if (!_has_verbose) throw Error("RecursiveSymmetryCutBase::StructureType: Verbose structure must be turned on to get dropped_mu() values.");
-    return _dropped_mu;
-  }
-  
-  // maximum symmetry value dropped
-  double max_dropped_symmetry() const {
-    if (!_has_verbose) throw Error("RecursiveSymmetryCutBase::StructureType: Verbose structure must be turned on to get max_dropped_symmetry().");
-    if (_dropped_symmetry.size() == 0) return 0.0;
-    return *std::max_element(_dropped_symmetry.begin(),_dropped_symmetry.end());
-  }
-  
+  void set_verbose(bool value) { _has_verbose = value;}
+
+  /// returns true if the current jet has some substructure (i.e. has
+  /// been tagged by the resursion) or not
+  ///
+  /// Note that this should include deltaR==0 (e.g. a perfectly
+  /// collinear branching with SoftDrop)
+  bool has_substructure() const { return _delta_R>=0; }
+  
+  /// number of dropped branches
+  int dropped_count(bool global=true) const;
+  
+  /// delta_R of dropped branches
+  /// when "global" is set, recurse into possile further substructure
+  std::vector<double> dropped_delta_R(bool global=true) const;
+  void set_dropped_delta_R(const std::vector<double> &v) { _dropped_delta_R = v; }
+
+  /// symmetry values of dropped branches
+  std::vector<double> dropped_symmetry(bool global=true) const;
+  void set_dropped_symmetry(const std::vector<double> &v) { _dropped_symmetry = v; }
+
+  /// mass drop values of dropped branches
+  std::vector<double> dropped_mu(bool global=true) const;
+  void set_dropped_mu(const std::vector<double> &v) { _dropped_mu = v; }
+  
+  /// maximum symmetry value dropped
+  double max_dropped_symmetry(bool global=true) const;
+
+  /// (global) list of groomed away elements as zg and thetag
+  /// sorted from largest to smallest anlge
+  std::vector<std::pair<double,double> > sorted_zg_and_thetag() const;
+
 private:
   double _delta_R, _symmetry, _mu;
   friend class RecursiveSymmetryCutBase;
 
+  bool _is_composite;
+  
   // additional verbose information
   bool _has_verbose;
@@ -291,4 +369,13 @@
   std::vector<double> _dropped_symmetry;
   std::vector<double> _dropped_mu;
+
+  bool check_verbose(const std::string &what) const{
+    if (!_has_verbose){
+      throw Error("RecursiveSymmetryCutBase::StructureType: Verbose structure must be turned on to get "+what+".");
+      return false;
+    }
+    return true;
+  }
+    
   
 };
Index: external/fastjet/contribs/RecursiveTools/SoftDrop.cc
===================================================================
--- external/fastjet/contribs/RecursiveTools/SoftDrop.cc	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/SoftDrop.cc	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,3 +1,3 @@
-// $Id: SoftDrop.cc 686 2014-06-14 03:25:09Z jthaler $
+// $Id: SoftDrop.cc 1059 2017-09-07 20:48:00Z gsoyez $
 //
 // Copyright (c) 2014-, Gregory Soyez, Jesse. Thaler
@@ -57,6 +57,8 @@
 //----------------------------------------------------------------------
 double SoftDrop::symmetry_cut_fn(const PseudoJet & p1, 
-                                 const PseudoJet & p2) const{
-  return _symmetry_cut * pow(p1.squared_distance(p2)/_R0sqr, 0.5*_beta);
+                                 const PseudoJet & p2,
+                                 void *optional_R0sqr_ptr) const{
+  double R0sqr = (optional_R0sqr_ptr == 0) ? _R0sqr : *((double*) optional_R0sqr_ptr);
+  return _symmetry_cut * pow(squared_geometric_distance(p1,p2)/R0sqr, 0.5*_beta);
 }
 
Index: external/fastjet/contribs/RecursiveTools/SoftDrop.hh
===================================================================
--- external/fastjet/contribs/RecursiveTools/SoftDrop.hh	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/SoftDrop.hh	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,3 +1,3 @@
-// $Id: SoftDrop.hh 686 2014-06-14 03:25:09Z jthaler $
+// $Id: SoftDrop.hh 1034 2017-08-01 10:03:53Z gsoyez $
 //
 // Copyright (c) 2014-, Gregory Soyez, Jesse Thaler
@@ -112,18 +112,28 @@
      RecursiveSymmetryCutBase(symmetry_measure, mu_cut, recursion_choice, subtractor),
      _beta(beta), _symmetry_cut(symmetry_cut), _R0sqr(R0*R0)
-  {}
+  {
+    // change the default: use grooming mode
+    set_grooming_mode();
+  }
 
   /// default destructor
   virtual ~SoftDrop(){}
 
+  //----------------------------------------------------------------------
+  // access to class info
+  double beta()         const { return _beta; }
+  double symmetry_cut() const { return _symmetry_cut; }
+  double R0()           const { return sqrt(_R0sqr); }
+  
 protected:
 
   // Unlike MMDT, the SoftDrop symmetry_cut_fn depends on the subjet kinematics
   // since the symmetry condition depends on the DeltaR between subjets.
-  virtual double symmetry_cut_fn(const PseudoJet & /* p1 */, 
-                                 const PseudoJet & /* p2 */) const;
+  virtual double symmetry_cut_fn(const PseudoJet & p1, 
+                                 const PseudoJet & p2,
+                                 void * optional_R0sqr_ptr = 0) const;
   virtual std::string symmetry_cut_description() const;
 
-private:
+  //private:
   double _beta;         ///< the power of the angular distance to be used
                         ///< in the symmetry condition
Index: external/fastjet/contribs/RecursiveTools/TODO
===================================================================
--- external/fastjet/contribs/RecursiveTools/TODO	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/TODO	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,2 +1,18 @@
+For v2.0
+--------
+
+- CHECK: lines 100-102 of RecursiveSymmetryCutBase.hh: this was
+  initially setting the structure members to 0. I think it's highly
+  preferable to set them to their default of -1, signalling the
+  absence of substructure. [0 dR might happen with a perfectly
+  collinear splitting]
+
+- CHECK: are we happy with how the structure info is stored and the
+  recursive_soft_drop_prongs(...) function in RecursiveSoftDrop.hh?
+
+- QUESTION: do we move set_min_deltaR_squared in the base class?
+
+- Add option to define z wrt to original jet?
+
 For future releases?
 --------------------
@@ -7,2 +23,4 @@
   auto-indenting. [Consider this for the full contrib??]
 - More generic kinematic cuts?
+- KZ: common base class for IteratedSoftDrop and RecursiveSoftDrop?
+- KZ: make IteratedSoftDrop use Recluster
Index: external/fastjet/contribs/RecursiveTools/VERSION
===================================================================
--- external/fastjet/contribs/RecursiveTools/VERSION	(revision 1f1f858be3b07de5df7883fba72ff0a8ab002b56)
+++ external/fastjet/contribs/RecursiveTools/VERSION	(revision 667a02abab10bed835e1186b81ddef9d3e67b782)
@@ -1,1 +1,1 @@
-1.0.0
+2.0.0-beta2
