Index: external/fastjet/tools/Recluster.cc
===================================================================
--- external/fastjet/tools/Recluster.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
+++ external/fastjet/tools/Recluster.cc	(revision 35cdc461e095e63247ab9dded070d76bc2c888e2)
@@ -0,0 +1,417 @@
+// $Id: Recluster.cc 3629 2014-08-14 17:21:15Z salam $
+//
+// Copyright (c) 2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see <http://www.gnu.org/licenses/>.
+//----------------------------------------------------------------------
+
+#include "fastjet/tools/Recluster.hh"
+#include "fastjet/CompositeJetStructure.hh"
+#include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
+#include <sstream>
+#include <typeinfo>
+
+using namespace std;
+
+// Comments:
+//
+//  - If the jet comes from a C/A clustering (or is a composite jet
+//    made of C/A clusterings) and we recluster it with a C/A
+//    algorithm, we just use exclusive jets instead of doing the
+//    clustering explicitly. In this specific case, we need to make
+//    sure that all the pieces share the same cluster sequence.
+//
+//  - If the recombiner has to be taken from the original jet and that
+//    jet is composite, we need to check that all the pieces were
+//    obtained with the same recombiner.
+//
+//  - Note that a preliminary version of this code has been
+//    distributed in the RecursiveTools fastjet-contrib. The present
+//    version has a few minor fixed
+
+FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+LimitedWarning Recluster::_explicit_ghost_warning;
+
+// ctor
+Recluster::Recluster(JetAlgorithm new_jet_alg, double new_jet_radius, Keep keep_in)
+  : _new_jet_def(JetDefinition(new_jet_alg, new_jet_radius)), 
+    _acquire_recombiner(true), _keep(keep_in), _cambridge_optimisation_enabled(true){}
+
+Recluster::Recluster(JetAlgorithm new_jet_alg, Keep keep_in)
+    : _acquire_recombiner(true), _keep(keep_in), _cambridge_optimisation_enabled(true){
+  switch (JetDefinition::n_parameters_for_algorithm(new_jet_alg)){
+  case 0: _new_jet_def = JetDefinition(new_jet_alg); break;
+  case 1: _new_jet_def = JetDefinition(new_jet_alg, JetDefinition::max_allowable_R); break;
+  default:
+    throw Error("Recluster(): tried to construct specifying only a jet algorithm ("+JetDefinition::algorithm_description(new_jet_alg)+") which takes more than 1 parameter");
+  };
+}
+
+
+// class description
+string Recluster::description() const {
+  ostringstream ostr;
+  ostr << "Recluster with new_jet_def = ";
+  if (_acquire_recombiner){
+    ostr << _new_jet_def.description_no_recombiner();
+    ostr << ", using a recombiner obtained from the jet being reclustered";
+  } else {
+    ostr << _new_jet_def.description();
+  }
+
+  if (_keep == keep_only_hardest)
+    ostr << " and keeping the hardest inclusive jet";
+  else
+    ostr << " and joining all inclusive jets into a composite jet";
+
+  return ostr.str();
+}
+
+
+// the main piece of code that performs the reclustering
+PseudoJet Recluster::result(const PseudoJet &jet) const {
+  // get the incljets and the exact jet definition that has been used
+  // to get them
+  vector<PseudoJet> incljets;
+  bool ca_optimised = get_new_jets_and_def(jet, incljets);
+
+  return generate_output_jet(incljets, ca_optimised);
+}
+
+
+// a lower-level method that does the actual work of reclustering the
+// input jet. The resulting incljets are stored in output_jets and the
+// jet definition that has been used can be deduced from their
+// associated ClusterSequence
+//
+//  - input_jet       the (input) jet that one wants to recluster
+//  - output_jets     incljets resulting from the new clustering
+//
+// returns true if the C/A optimisation has been used (this means
+// that generate_output_jet will watch out for non-explicit-ghost
+// areas that might be leftover)
+bool Recluster::get_new_jets_and_def(const PseudoJet & input_jet, 
+                                     vector<PseudoJet> & output_jets) const{
+  // generic sanity checks
+  //-------------------------------------------------------------------
+  // make sure that the jet has constituents
+  if (! input_jet.has_constituents())
+    throw Error("Recluster can only be applied on jets having constituents");
+
+  // tests particular to certain configurations
+  //-------------------------------------------------------------------
+
+  // for a whole variety of tests, we shall need the "recursive"
+  // pieces of the jet (the jet itself or recursing down to its most
+  // fundamental pieces). So we do compute these once and for all.
+  //
+  // Note that the pieces are always needed (either for C/A or for the
+  // area checks)
+  vector<PseudoJet> all_pieces;
+  if ((!_get_all_pieces(input_jet, all_pieces)) || (all_pieces.size()==0)){
+    throw Error("Recluster: failed to retrieve all the pieces composing the jet.");
+  }
+
+  // decide which jet definition to use
+  //-------------------------------------------------------------------
+  JetDefinition new_jet_def = _new_jet_def;
+  if (_acquire_recombiner){
+    _acquire_recombiner_from_pieces(all_pieces, new_jet_def);
+  }
+
+  // the vector that will ultimately hold the incljets
+  output_jets.clear();
+
+  // check if we can apply the simplification for C/A jets reclustered
+  // with C/A
+  //
+  // we apply C/A clustering iff
+  //  - the requested new_jet_def is C/A
+  //  - the jet is either directly coming from C/A or if it is a
+  //    superposition of C/A jets from the same cluster sequence
+  //  - the pieces agree with the recombination scheme of new_jet_def
+  //
+  // In this case area support will be automatically inherited so we
+  // can only worry about this later
+  // -------------------------------------------------------------------
+  if (_check_ca(all_pieces, new_jet_def)){
+    _recluster_ca(all_pieces, output_jets, new_jet_def.R());
+    output_jets = sorted_by_pt(output_jets);
+    return true;
+  }
+
+  // decide if area support has to be kept
+  //-------------------------------------------------------------------
+  bool include_area_support = input_jet.has_area();
+  if ((include_area_support) &&  (!_check_explicit_ghosts(all_pieces))){
+    _explicit_ghost_warning.warn("Recluster: the original cluster sequence is lacking explicit ghosts; area support will no longer be available after re-clustering");
+    include_area_support = false;
+  }
+
+  // extract the incljets
+  //-------------------------------------------------------------------
+  _recluster_generic(input_jet, output_jets, new_jet_def, include_area_support);
+  output_jets = sorted_by_pt(output_jets);
+
+  return false;
+}
+
+// given a set of incljets and a jet definition used, create the
+// resulting PseudoJet
+PseudoJet Recluster::generate_output_jet(std::vector<PseudoJet> & incljets, 
+                                         bool ca_optimisation_used) const{
+  // first handle the case where we only need to keep the hardest incljet
+  if (_keep == keep_only_hardest) {
+    if (incljets.size() > 0) {
+      return incljets[0];
+    } else {
+      return PseudoJet();
+    }
+  }
+
+  // now the case where all incljets have to be joined
+
+  // safekeeper
+  if (incljets.size()==0) return join(incljets);
+
+  PseudoJet reclustered = join(incljets, 
+			       *(incljets[0].associated_cluster_sequence()->jet_def().recombiner()));
+
+  // if we've used C/A optimisation, we need to get rid of the area
+  // information if it comes from a non-explicit-ghost clustering.
+  // (because in that case it can be erroneous due to the lack of
+  // information about empty areas)
+  if (ca_optimisation_used){
+    if (reclustered.has_area() &&
+        (incljets.size() > 0) &&
+        (! incljets[0].validated_csab()->has_explicit_ghosts())){
+      CompositeJetStructure *css = dynamic_cast<CompositeJetStructure *>(reclustered.structure_non_const_ptr());
+      assert(css);
+      css->discard_area();
+    }
+  }
+  
+  return reclustered;
+}
+
+//----------------------------------------------------------------------
+// the parts that really do the reclustering
+//----------------------------------------------------------------------
+
+// get the subjets in the simple case of C/A+C/A
+void Recluster::_recluster_ca(const vector<PseudoJet> & all_pieces, 
+                              vector<PseudoJet> & subjets, 
+                              double Rfilt) const{
+  subjets.clear();
+
+  // each individual piece should have a C/A cluster sequence
+  // associated to it. So a simple loop would do the job
+  for (vector<PseudoJet>::const_iterator piece_it = all_pieces.begin(); 
+       piece_it!=all_pieces.end(); piece_it++){
+    // just extract the exclusive subjets of 'jet'
+    const ClusterSequence *cs = piece_it->associated_cluster_sequence(); 
+    vector<PseudoJet> local_subjets;
+
+    double dcut = Rfilt / cs->jet_def().R();
+    if (dcut>=1.0){
+      // remember that in this case all the pairwise interpiece
+      // distances are supposed to be larger than Rfilt (this was
+      // tested in _check_ca), which means that they can never
+      // recombine with each other.
+      local_subjets.push_back(*piece_it);
+    } else {
+      local_subjets = piece_it->exclusive_subjets(dcut*dcut);
+    }
+
+    copy(local_subjets.begin(), local_subjets.end(), back_inserter(subjets));
+  }
+}
+
+
+// perform the reclustering itself for all cases where the "C/A trick"
+// does not apply
+void Recluster::_recluster_generic(const PseudoJet & jet, 
+                                   vector<PseudoJet> & incljets,
+                                   const JetDefinition & new_jet_def,
+                                   bool do_areas) const{
+  // create a new, internal, ClusterSequence from the jet constituents
+  // get the incljets directly from there
+  //
+  // If the jet has area support then we separate the ghosts from the
+  // "regular" particles so the incljets will also have area
+  // support. Note that we do this regardless of whether rho is zero
+  // or not.
+  //
+  // Note that to be able to separate the ghosts, one needs explicit
+  // ghosts!!
+  // ---------------------------------------------------------------
+  if (do_areas){
+    //vector<PseudoJet> all_constituents = jet.constituents();
+    vector<PseudoJet> regular_constituents, ghosts;  
+    SelectorIsPureGhost().sift(jet.constituents(), ghosts, regular_constituents);
+
+    // figure out the ghost area from the 1st ghost (if none, any value
+    // would probably do as the area will be 0 and subtraction will have
+    // no effect!)
+    double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
+    ClusterSequenceActiveAreaExplicitGhosts * csa
+      = new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents, 
+                                                    new_jet_def, 
+                                                    ghosts, ghost_area);
+
+    incljets = csa->inclusive_jets();
+
+    // allow the cs to be deleted when it's no longer used
+    // 
+    // Note that there is at least one constituent in the jet so there
+    // is in principle at least one incljet. But one may have used a
+    // nasty recombiner or jet def that left an empty set of incljets,
+    // so we'd rather play it safe (e.g. GridJetPlugin, with the
+    // constituents outside the range of the grid)
+    if (incljets.size())
+      csa->delete_self_when_unused();
+    else
+      delete csa;
+  } else {
+    ClusterSequence * cs = new ClusterSequence(jet.constituents(), new_jet_def);
+    incljets = cs->inclusive_jets();
+    // allow the cs to be deleted when it's no longer used (again, we
+    // add an extra safety check)
+    if (incljets.size())
+      cs->delete_self_when_unused();
+    else 
+      delete cs;
+  }
+}
+
+//----------------------------------------------------------------------
+// various checks and internal constructs
+//----------------------------------------------------------------------
+
+// fundamental info for CompositeJets
+//----------------------------------------------------------------------
+ 
+// get the pieces down to the fundamental pieces
+// 
+// Note that this just checks that there is an associated CS to the
+// fundamental pieces, not that it is still valid
+bool Recluster::_get_all_pieces(const PseudoJet &jet, vector<PseudoJet> &all_pieces) const{
+  if (jet.has_associated_cluster_sequence()){
+    all_pieces.push_back(jet);
+    return true;
+  }
+
+  if (jet.has_pieces()){
+    const vector<PseudoJet> pieces = jet.pieces();
+    for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
+      if (!_get_all_pieces(*it, all_pieces)) return false;
+    return true;
+  }
+
+  return false;
+}
+
+// construct the re-clustering jet definition using the recombiner
+// from whatever definition has been used to obtain the original jet
+//----------------------------------------------------------------------
+void Recluster::_acquire_recombiner_from_pieces(const vector<PseudoJet> &all_pieces, 
+                                                JetDefinition &new_jet_def) const{
+  // a quick safety check
+  assert(_acquire_recombiner);
+
+  // check that all the pieces have the same recombiner
+  //
+  // Note that if the jet has an associated cluster sequence that is no
+  // longer valid, an error will be thrown (needed since it could be the
+  // 1st check called after the enumeration of the pieces)
+  const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
+  for (unsigned int i=1; i<all_pieces.size(); i++){
+    if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref)){
+      throw Error("Recluster instance is configured to determine the recombination scheme (or recombiner) from the original jet, but different pieces of the jet were found to have non-equivalent recombiners.");
+    }
+  }
+
+  // get the recombiner from the original jet_def
+  new_jet_def.set_recombiner(jd_ref);
+}
+
+// area support
+//----------------------------------------------------------------------
+
+// check if the jet (or all its pieces) have explicit ghosts
+// (assuming the jet has area support).
+//
+// Note that if the jet has an associated cluster sequence that is no
+// longer valid, an error will be thrown (needed since it could be the
+// 1st check called after the enumeration of the pieces)
+bool Recluster::_check_explicit_ghosts(const vector<PseudoJet> &all_pieces) const{
+  for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
+    if (! it->validated_csab()->has_explicit_ghosts()) return false;
+  return true;
+}
+
+// C/A specific tests
+//----------------------------------------------------------------------
+
+// check if one can apply the simplification for C/A incljets
+//
+// This includes:
+//  - the incljet definition asks for C/A incljets
+//  - all the pieces share the same CS
+//  - that CS is C/A with the same recombiner as the incljet def
+//  - the re-clustering radius is not larger than any of the pairwise
+//    distance between the pieces
+//
+// Note that if the jet has an associated cluster sequence that is no
+// longer valid, an error will be thrown (needed since it could be the
+// 1st check called after the enumeration of the pieces)
+bool Recluster::_check_ca(const vector<PseudoJet> &all_pieces, 
+                          const JetDefinition &new_jet_def) const{
+  // check that optimisation is enabled
+  if (!_cambridge_optimisation_enabled) return false;
+
+  // check that we're reclustering with C/A
+  if (new_jet_def.jet_algorithm() != cambridge_algorithm) return false;
+
+  // check that the 1st of all the pieces (we're sure there is at
+  // least one) is coming from a C/A clustering. Then check that all
+  // the following pieces share the same ClusterSequence
+  const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
+  if (cs_ref->jet_def().jet_algorithm() != cambridge_algorithm) return false;
+  for (unsigned int i=1; i<all_pieces.size(); i++)
+    if (all_pieces[i].validated_cs() != cs_ref) return false;
+
+  // check that the 1st piece has the same recombiner as the one used
+  // for the incljet clustering
+  // Note that since they share the same CS, checking the 1st one is enough
+  if (!cs_ref->jet_def().has_same_recombiner(new_jet_def)) return false;
+
+  // we also have to make sure that the reclustering radius is not larger
+  // than any of the inter-piece distances
+  double Rnew2 = new_jet_def.R();
+  Rnew2 *= Rnew2;
+  for (unsigned int i=0; i<all_pieces.size()-1; i++){
+    for (unsigned int j=i+1; j<all_pieces.size(); j++){
+      if (all_pieces[i].squared_distance(all_pieces[j]) <  Rnew2) return false;
+    }
+  }
+
+  return true;
+}
+
+FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
