Changeset 273e668 in git for external/fastjet/plugins/SISCone

Timestamp:

Oct 15, 2014, 10:55:55 AM (10 years ago)

Author:

Pavel Demin <pavel.demin@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

35b9204, b25d4cf

Parents:

f14bd6a

git-author:

Pavel Demin <pavel.demin@…> (10/10/14 08:56:40)

git-committer:

Pavel Demin <pavel.demin@…> (10/15/14 10:55:55)

Message:

upgrade FastJet to version 3.1.0

Location:

external/fastjet/plugins/SISCone

Files:

• SISConePlugin.cc (modified) (6 diffs)
• config.h (modified) (2 diffs)
• fastjet/SISConeBasePlugin.hh (modified) (6 diffs)
• fastjet/SISConePlugin.hh (modified) (6 diffs)
• siscone.cc (modified) (5 diffs)
• siscone.h (modified) (4 diffs)
• split_merge.cc (modified) (6 diffs)
• split_merge.h (modified) (7 diffs)

external/fastjet/plugins/SISCone/SISConePlugin.cc

@@ -21 +21 @@
                       four_vector.E);
 }
+
+//======================================================================
+// wrap-up around siscone's user-defined scales
+namespace siscone_plugin_internal{
+  /// @ingroup internal
+  /// \class SISConeUserScale
+  /// class that makes the transition between the internal SISCone
+  /// user-defined scale choice (using SISCone's Cjet) and
+  /// user-defined scale choices in the plugn above (using FastJet's
+  /// PseudoJets)
+  class SISConeUserScale  : public siscone::Csplit_merge::Cuser_scale_base{
+  public:
+    /// ctor takes the "fastjet-style" user-defined scale as well as a
+    /// reference to the current cluster sequence (to access the
+    /// particles if needed)
+    SISConeUserScale(const SISConePlugin::UserScaleBase *user_scale,
+                     const ClusterSequence &cs)
+      : _user_scale(user_scale), _cs(cs){}
+
+    /// returns the scale associated to a given jet
+    virtual double operator()(const siscone::Cjet &jet) const{
+      return _user_scale->result(_build_PJ_from_Cjet(jet));
+    }
+
+    /// returns true id the scasle associated to jet a is larger than
+    /// the scale associated to jet b
+    virtual bool is_larger(const siscone::Cjet &a, const siscone::Cjet &b) const{
+      return _user_scale->is_larger(_build_PJ_from_Cjet(a), _build_PJ_from_Cjet(b));
+    }
+
+  private:
+    /// constructs a PseudoJet from a siscone::Cjet
+    ///
+    /// Note that it is tempting to overload the PseudoJet ctor. This
+    /// would not work because down the line we need to access the
+    /// original PseudoJet through the ClusterSequence and therefore
+    /// the PseudoJet structure needs to be aware of the
+    /// ClusterSequence.
+    PseudoJet _build_PJ_from_Cjet(const siscone::Cjet &jet) const{
+      PseudoJet j(jet.v.px, jet.v.py, jet.v.pz, jet.v.E);
+      j.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(
+                                   new SISConePlugin::UserScaleBaseStructureType<siscone::Cjet>(jet,_cs)));      
+      return j;
+    }
+
+    const SISConePlugin::UserScaleBase *_user_scale;
+    const ClusterSequence & _cs;
+  };
+}
+
+// end of the internal material
+//======================================================================
 
 
@@ -45 +97 @@
   desc << "SISCone jet algorithm with " ;
   desc << "cone_radius = "       << cone_radius        () << ", ";
-  desc << "overlap_threshold = " << overlap_threshold  () << ", ";
+  if (_progressive_removal)
+    desc << "progressive-removal mode, ";
+  else  
+    desc << "overlap_threshold = " << overlap_threshold  () << ", ";
   desc << "n_pass_max = "        << n_pass_max         () << ", ";
   desc << "protojet_ptmin = "    << protojet_ptmin()      << ", ";
-  desc <<  sm_scale_string                                << ", ";
-  desc << "caching turned "      << (caching() ? on : off);
-  desc << ", SM stop scale = "     << _split_merge_stopping_scale;
+  if (_progressive_removal && _user_scale) {
+    desc << "using a user-defined scale for ordering of stable cones";
+    string user_scale_desc = _user_scale->description();
+    if (user_scale_desc != "") desc << " (" << user_scale_desc << ")";
+  } else {
+    desc <<  sm_scale_string;
+  }
+  if (!_progressive_removal){
+    desc << ", caching turned "      << (caching() ? on : off);
+    desc << ", SM stop scale = "     << _split_merge_stopping_scale;
+  }
 
   // add a note to the description if we use the pt-weighted splitting
@@ -85 +148 @@
   bool new_siscone = true; // by default we'll be running it
 
-  if (caching()) {
+  if (caching() && !_progressive_removal) {
 
     // Establish if we have a cached run with the same R, npass and
@@ -138 +201 @@
     // run the jet finding
     //cout << "plg sms: " << split_merge_scale() << endl;
-    siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
-                          n_pass_max(), protojet_or_ghost_ptmin(), 
-                          Esplit_merge_scale(split_merge_scale()));
+    if (_progressive_removal){
+      // handle the optional user-defined scale choice
+      SharedPtr<siscone_plugin_internal::SISConeUserScale> internal_scale;
+      if (_user_scale){
+        internal_scale.reset(new siscone_plugin_internal::SISConeUserScale(_user_scale, clust_seq));
+        siscone->set_user_scale(internal_scale.get());
+      }
+      siscone->compute_jets_progressive_removal(siscone_momenta, cone_radius(),
+                                                n_pass_max(), protojet_or_ghost_ptmin(), 
+                                                Esplit_merge_scale(split_merge_scale()));
+    } else {
+      siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
+                            n_pass_max(), protojet_or_ghost_ptmin(), 
+                            Esplit_merge_scale(split_merge_scale()));
+    }
   } else {
     // rerun the jet finding
@@ -155 +230 @@
   SISConeExtras * extras = new SISConeExtras(n);
 
+  // the ordering in which the inclusive jets are transfered here is
+  // deliberate and ensures that when a user asks for
+  // inclusive_jets(), they are provided in the order in which SISCone
+  // created them.
   for (int ijet = njet-1; ijet >= 0; ijet--) {
     const Cjet & jet = siscone->jets[ijet]; // shorthand
@@ -220 +299 @@
 }
 
+// //======================================================================
+// // material to handle user-defined scales
+// 
+// //--------------------------------------------------
+// // SISCone structure type
+// 
+// // the textual descripotion
+// std::string SISConePlugin::UserScaleBase::StructureType::description() const{
+//   return "PseudoJet wrapping a siscone::Cjet stable cone"; 
+// }
+// 
+// // retrieve the constituents
+// //
+// // if you simply need to iterate over the constituents, it will be
+// // faster to access them via constituent(i)
+// vector<PseudoJet> SISConePlugin::UserScaleBase::StructureType::constituents(const PseudoJet &) const{ 
+//   vector<PseudoJet> constits;
+//   constits.reserve(_jet.n);
+//   for (unsigned int i=0; i<(unsigned int)_jet.n;i++)
+//     constits.push_back(constituent(i));
+//   return constits;
+// }
+// 
+// // returns the number of constituents
+// unsigned int SISConePlugin::UserScaleBase::StructureType::size() const{
+//   return _jet.n;
+// }
+// 
+// // returns the index (in the original particle list) of the ith
+// // constituent
+// int SISConePlugin::UserScaleBase::StructureType::constituent_index(unsigned int i) const{ 
+//   return _jet.contents[i];
+// }
+// 
+// // returns the ith constituent (as a PseusoJet)
+// const PseudoJet & SISConePlugin::UserScaleBase::StructureType::constituent(unsigned int i) const{
+//   return _cs.jets()[_jet.contents[i]];
+// }
+// 
+// // returns the scalar pt of this stable cone
+// double SISConePlugin::UserScaleBase::StructureType::pt_tilde() const{
+//   return _jet.pt_tilde;
+// }
+// 
+// // returns the sm_var2 (signed ordering variable squared) for this stable cone
+// double SISConePlugin::UserScaleBase::StructureType::ordering_var2() const{
+//   return _jet.sm_var2;
+// }
+
+
 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh

external/fastjet/plugins/SISCone/config.h

@@ -49 +49 @@
 
 /* Define to the full name and version of this package. */
-#define PACKAGE_STRING "SISCone 2.0.6"
+#define PACKAGE_STRING "SISCone 3.0.0"
 
 /* Define to the one symbol short name of this package. */
 #define PACKAGE_TARNAME "siscone"
 
+/* Define to the home page for this package. */
+/* #undef PACKAGE_URL */
+
 /* Define to the version of this package. */
-#define PACKAGE_VERSION "2.0.6"
+#define PACKAGE_VERSION "3.0.0"
 
 /* Define to 1 if you have the ANSI C header files. */
@@ -61 +64 @@
 
 /* Version number of package */
-#define VERSION "2.0.6"
+#define VERSION "3.0.0"

external/fastjet/plugins/SISCone/fastjet/SISConeBasePlugin.hh

@@ -43 +43 @@
   SISConeBasePlugin (){
     _use_jet_def_recombiner = false;
+    set_progressive_removal(false);
   }
 
@@ -49 +50 @@
     *this = plugin;
   }
+
+  /// set whether to use SISCone with progressive removal instead of
+  /// the default split_merge step.
+  ///
+  /// If progressive removal is enabled, the following SISCone
+  /// variables are not used:
+  ///
+  /// - overlap_threshold
+  /// - caching
+  /// - split_merge_stopping_scale
+  ///
+  /// The split_merge_scale choice is reinterpreted as the ordering
+  /// variable for progressive removal. It is also possible for the
+  /// user to supply his/her own function for the scale that orders
+  /// progressive removal, with set_user_scale(...)
+  void set_progressive_removal(bool progressive_removal_in=true){
+    _progressive_removal = progressive_removal_in;
+  }
+
+  /// returns true if progressive_removal is enabled
+  bool progressive_removal() const{ return _progressive_removal;}
 
   /// the cone radius
@@ -105 +127 @@
   }
 
+  // user-defined scale for progressive removal
+  //------------------------------------------------------------
+
+  /// \class UserScaleBase
+  /// base class for user-defined ordering of stable cones (used for
+  /// prorgessive removal)
+  ///
+  /// derived classes have to implement the () operator that returns
+  /// the scale associated with a given jet.
+  /// 
+  /// It is also highly recommended to implement the is_larger()
+  /// method whenever possible, in order to avoid rounding issues
+  /// known to lead to possible infrared unsafeties.
+  ///
+  /// The jets that are passed to this class will carry the structure
+  /// of type SISConePlugin::StructureType which allows to retreive
+  /// easily the following information:
+  ///
+  ///   vector<PseudoJet> constituents = jet.constituents();
+  ///   unsigned int n_constituents = jet.structure_of<SISConePlugin::UserScaleBase>().size();
+  ///   int index = jet.structure_of<SISConePlugin::UserScaleBase>().constituent_index(index i);
+  ///   const PseudoJet & p = jet.structure_of<SISConePlugin::UserScaleBase>().constituent(index i);
+  ///   double scalar_pt = jet.structure_of<SISConePlugin::UserScaleBase>().pt_tilde();
+  ///
+  /// see SISConePlugin::StructureType below for further details
+  class UserScaleBase : public FunctionOfPseudoJet<double>{
+  public:
+    /// empty virtual dtor
+    virtual ~UserScaleBase(){}
+
+    /// returns the scale associated with a given jet
+    ///
+    /// "progressive removal" iteratively removes the stable cone with
+    /// the largest scale
+    virtual double result(const PseudoJet & jet) const = 0;
+
+    /// returns true when the scale associated with jet a is larger than
+    /// the scale associated with jet b
+    ///
+    /// By default this does a simple direct comparison but it can be
+    /// overloaded for higher precision [recommended if possible]
+    virtual bool is_larger(const PseudoJet & a, const PseudoJet & b) const;
+
+    class StructureType; // defined below
+  };
+
+  // template class derived from UserScaleBase::StryctureType that
+  // works for both SISCone jet classes
+  // implemented below 
+  template<class Tjet> 
+  class UserScaleBaseStructureType;
+
+  /// set a user-defined scale for stable-cone ordering in
+  /// progressive removal
+  void set_user_scale(const UserScaleBase *user_scale_in){ _user_scale = user_scale_in;}
+
+  /// returns the user-defined scale in use (0 if none)
+  const UserScaleBase * user_scale() const{ return _user_scale;}
+
+
   // the things that one MUST overload required by base class
   //---------------------------------------------------------
@@ -120 +202 @@
   double _split_merge_stopping_scale;
   bool   _use_jet_def_recombiner;
+  bool   _progressive_removal;
 
   mutable double _ghost_sep_scale;
@@ -126 +209 @@
   /// call the re-clustering itself 
   virtual void reset_stored_plugin() const =0;
+
+  const UserScaleBase * _user_scale;
 
 };
@@ -185 +270 @@
 inline SISConeBaseExtras::~SISConeBaseExtras(){}
 
+//----------------------------------------------------------------------
+// implementation of the structure type associated with the UserScaleBase class
+
+/// \class SISConeBasePlugin::UserScaleBase::StructureType
+/// the structure that allows to store the information contained
+/// into a siscone::Cjet (built internally in SISCone from a stable
+/// cone) into a PseudoJet
+class SISConeBasePlugin::UserScaleBase::StructureType : public PseudoJetStructureBase {
+public:
+  /// base ctor (constructed from a ClusterSequence tin order to have
+  /// access to the initial particles
+  StructureType(const ClusterSequence &cs)
+    : _cs(cs){}
+
+  /// empty virtual dtor 
+  virtual ~StructureType(){}
+  
+  //--------------------------------------------------
+  // members inherited from the base class
+  /// the textual descripotion
+  virtual std::string description() const{
+    return "PseudoJet wrapping a siscone jet from a stable cone"; 
+  }
+
+  /// this structure has constituents
+  virtual bool has_constituents() const {return true;}
+
+  /// retrieve the constituents 
+  ///
+  /// if you simply need to iterate over the constituents, it will be
+  /// faster to access them via constituent(i)
+  virtual std::vector<PseudoJet> constituents(const PseudoJet & /*reference*/) const{ 
+    std::vector<PseudoJet> constits;
+    constits.reserve(size());
+    for (unsigned int i=0; i<size();i++)
+      constits.push_back(constituent(i));
+    return constits;
+  }
+  
+  //--------------------------------------------------
+  // additional information relevant for this structure
+
+  /// returns the number of constituents
+  virtual unsigned int size() const = 0;
+
+  /// returns the index (in the original particle list) of the ith
+  /// constituent
+  virtual int constituent_index(unsigned int i) const = 0;
+
+  /// returns the ith constituent (as a PseusoJet)
+  const PseudoJet & constituent(unsigned int i) const{
+    return _cs.jets()[constituent_index(i)];
+  }
+
+  // /// returns the scalar pt of this stable cone
+  // virtual double pt_tilde() const = 0;
+
+  /// returns the sm_var2 (signed ordering variable squared) for this stable cone
+  virtual double ordering_var2() const = 0;
+
+protected:
+  const ClusterSequence &_cs; ///< a reference to the CS (for access to the particles)
+};
+
+
+///@ingroup internal
+/// template class derived from UserScaleBase::StryctureType that
+/// works for both SISCone jet classes
+/// implemented below 
+template<class Tjet>
+class SISConeBasePlugin::UserScaleBaseStructureType : public UserScaleBase::StructureType{
+public:
+  UserScaleBaseStructureType(const Tjet &jet, const ClusterSequence &cs)
+    : UserScaleBase::StructureType(cs), _jet(jet){}
+
+  /// empty virtual dtor 
+  virtual ~UserScaleBaseStructureType(){}
+
+  //--------------------------------------------------
+  // additional information relevant for this structure
+
+  /// returns the number of constituents
+  virtual unsigned int size() const{
+    return _jet.n;
+  }
+
+  /// returns the index (in the original particle list) of the ith
+  /// constituent
+  virtual int constituent_index(unsigned int i) const{
+    return _jet.contents[i];
+  }
+
+  // /// returns the scalar pt of this stable cone
+  // virtual double pt_tilde() const{
+  //   return _jet.pt_tilde;
+  // }
+
+  /// returns the sm_var2 (signed ordering variable squared) for this stable cone
+  virtual double ordering_var2() const{
+    return _jet.sm_var2;
+  }
+
+protected:
+  const Tjet &_jet; ///< a reference to the internal jet in SISCone
+};
+
 
 FASTJET_END_NAMESPACE        // defined in fastjet/internal/base.hh

external/fastjet/plugins/SISCone/fastjet/SISConePlugin.hh

@@ -7 +7 @@
 namespace siscone{
   class Csiscone;
+  class Cjet;
 }
 
@@ -74 +75 @@
   /// enum for the different split-merge scale choices;
   /// Note that order _must_ be the same as in siscone
-  enum SplitMergeScale {SM_pt,     ///< transverse momentum (E-scheme), IR unsafe
-                        SM_Et,     ///< transverse energy (E-scheme), not long. boost invariant
-                                   ///< original run-II choice [may not be implemented]
-                        SM_mt,     ///< transverse mass (E-scheme), IR safe except
-                                   ///< in decays of two identical narrow heavy particles
-                        SM_pttilde ///< pt-scheme pt = \sum_{i in jet} |p_{ti}|, should
-                                   ///< be IR safe in all cases
+  enum SplitMergeScale {
+    SM_pt,     ///< transverse momentum (E-scheme), IR unsafe
+    SM_Et,     ///< transverse energy (E-scheme), not long. boost invariant
+               ///< original run-II choice [may not be implemented]
+    SM_mt,     ///< transverse mass (E-scheme), IR safe except
+               ///< in decays of two identical narrow heavy particles
+    SM_pttilde ///< pt-scheme pt = \sum_{i in jet} |p_{ti}|, should
+               ///< be IR safe in all cases
   };
 
@@ -108 +110 @@
     _split_merge_stopping_scale = split_merge_stopping_scale_in;
     _ghost_sep_scale       = 0.0;
-    _use_pt_weighted_splitting = false;}
+    _use_pt_weighted_splitting = false;
+    _user_scale = 0;}
 
 
@@ -125 +128 @@
     _split_merge_stopping_scale = 0.0;
     _split_merge_scale     = split_merge_on_transverse_mass_in ? SM_mt : SM_pttilde;
-    _ghost_sep_scale       = 0.0;}
+    _ghost_sep_scale       = 0.0;
+    _user_scale = 0;}
   
   /// backwards compatible constructor for the SISCone Plugin class
@@ -141 +145 @@
     _split_merge_stopping_scale = 0.0;
     _ghost_sep_scale       = 0.0;
-    _use_pt_weighted_splitting = false;}
+    _use_pt_weighted_splitting = false;
+    _user_scale = 0;}
 
   /// minimum pt for a protojet to be considered in the split-merge step
@@ -191 +196 @@
 };
 
+
+/////\class SISConePlugin::UserScaleBase::StructureType
+///// the structure that allows to store the information contained
+///// into a siscone::Cjet (built internally in SISCone from a stable
+///// cone) into a PseudoJet
+//class SISConePlugin::UserScaleBase::StructureType : public PseudoJetStructureBase {
+//public:
+//  StructureType(const siscone::Cjet & jet, const ClusterSequence &cs)
+//    : _jet(jet), _cs(cs){}
+//
+//  //--------------------------------------------------
+//  // members inherited from the base class
+//  /// the textual descripotion
+//  virtual std::string description() const;
+//
+//  /// this structure has constituents
+//  virtual bool has_constituents() const {return true;}
+//
+//  /// retrieve the constituents 
+//  ///
+//  /// if you simply need to iterate over the constituents, it will be
+//  /// faster to access them via constituent(i)
+//  virtual std::vector<PseudoJet> constituents(const PseudoJet & /*reference*/) const;
+//
+//  //--------------------------------------------------
+//  // additional information relevant for this structure
+//
+//  /// returns the number of constituents
+//  unsigned int size() const;
+//
+//  /// returns the index (in the original particle list) of the ith
+//  /// constituent
+//  int constituent_index(unsigned int i) const;
+//
+//  /// returns the ith constituent (as a PseusoJet)
+//  const PseudoJet & constituent(unsigned int i) const;
+//
+//  /// returns the scalar pt of this stable cone
+//  double pt_tilde() const;
+//
+//  /// returns the sm_var2 (signed ordering variable squared) for this stable cone
+//  double ordering_var2() const;
+//
+//protected:
+//  const siscone::Cjet &_jet;  ///< a dreference to the internal jet in SISCone
+//  const ClusterSequence &_cs; ///< a reference to the CS (for access to the particles)
+//};
 
 //======================================================================

external/fastjet/plugins/SISCone/siscone.cc

@@ -21 +21 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 341                                                          $//
-// $Date:: 2013-04-08 12:21:06 +0200 (Mon, 08 Apr 2013)                     $//
+// $Revision:: 371                                                          $//
+// $Date:: 2014-09-09 10:05:32 +0200 (Tue, 09 Sep 2014)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -29 +29 @@
 //#else
 //#define PACKAGE_NAME "SISCone"
-//#define VERSION "2.0.6"
+//#define VERSION "3.0.0"
 //#warning "No config.h file available, using preset values"
 //#endif
@@ -87 +87 @@
                            int _n_pass_max, double _ptmin,
                            Esplit_merge_scale _split_merge_scale){
-  // initialise random number generator
-  if (!init_done){
-    // initialise random number generator
-    ranlux_init();
-
-    // do not do this again
-    init_done=true;
-
-    // print the banner
-    if (_banner_ostr != 0){
-      (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl;
-      (*_banner_ostr) << "#                    SISCone   version " << setw(28) << left << siscone_version() << "o" << endl;
-      (*_banner_ostr) << "#              http://projects.hepforge.org/siscone                o" << endl;
-      (*_banner_ostr) << "#                                                                  o" << endl;
-      (*_banner_ostr) << "# This is SISCone: the Seedless Infrared Safe Cone Jet Algorithm   o" << endl;
-      (*_banner_ostr) << "# SISCone was written by Gavin Salam and Gregory Soyez             o" << endl;
-      (*_banner_ostr) << "# It is released under the terms of the GNU General Public License o" << endl;
-      (*_banner_ostr) << "#                                                                  o" << endl;
-      (*_banner_ostr) << "# A description of the algorithm is available in the publication   o" << endl;
-      (*_banner_ostr) << "# JHEP 05 (2007) 086 [arXiv:0704.0292 (hep-ph)].                   o" << endl;
-      (*_banner_ostr) << "# Please cite it if you use SISCone.                               o" << endl;
-      (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl;
-      (*_banner_ostr) << endl;
-
-      _banner_ostr->flush();
-    }
-  }
+  _initialise_if_needed();
 
   // run some general safety tests (NB: f will be checked in split-merge)
@@ -174 +148 @@
 }
 
+
+/*
+ * compute the jets from a given particle set doing multiple passes
+ * such pass N looks for jets among all particles not put into jets
+ * during previous passes.
+ *  - _particles   list of particles
+ *  - _radius      cone radius
+ *  - _n_pass_max  maximum number of runs
+ *  - _ptmin       minimum pT of the protojets
+ *  - _ordering_scale    the ordering scale to decide which stable
+ *                       cone is removed
+ * return the number of jets found.
+ **********************************************************************/
+int Csiscone::compute_jets_progressive_removal(vector<Cmomentum> &_particles, double _radius, 
+                                               int _n_pass_max, double _ptmin,
+                                               Esplit_merge_scale _ordering_scale){
+  _initialise_if_needed();
+
+  // run some general safety tests (NB: f will be checked in split-merge)
+  if (_radius <= 0.0 || _radius >= 0.5*M_PI) {
+    ostringstream message;
+    message << "Illegal value for cone radius, R = " << _radius 
+            << " (legal values are 0<R<pi/2)";
+    throw Csiscone_error(message.str());
+  }
+
+  ptcomparison.split_merge_scale = _ordering_scale;
+  partial_clear(); // make sure some things are initialised properly
+
+  // init the split_merge algorithm with the initial list of particles
+  // this initialises particle list p_left of remaining particles to deal with
+  //
+  // this stores the "processed" particles in p_uncol_hard
+  init_particles(_particles);
+  jets.clear();
+
+  bool unclustered_left;
+  rerun_allowed = false;
+  protocones_list.clear();
+
+  do{
+    //cout << n_left << " particle left" << endl; 
+
+    // initialise stable_cone finder
+    // here we use the list of remaining particles
+    // AFTER COLLINEAR CLUSTERING !!!!!!
+    Cstable_cones::init(p_uncol_hard);
+
+    // get stable cones (stored in 'protocones')
+    unclustered_left = get_stable_cones(_radius);
+
+    // add the hardest stable cone to the list of jets
+    if (add_hardest_protocone_to_jets(&protocones, R2, _ptmin)) break;
+  
+    _n_pass_max--;
+  } while ((unclustered_left) && (n_left>0) && (_n_pass_max!=0));
+
+  // split & merge
+  return jets.size();
+}
+
+
 /*
  * recompute the jets with a different overlap parameter.
@@ -205 +241 @@
 }  
 
+// ensure things are initialised
+void Csiscone::_initialise_if_needed(){
+  // initialise random number generator
+  if (init_done) return;
+
+  // initialise random number generator
+  ranlux_init();
+
+  // do not do this again
+  init_done=true;
+
+  // print the banner
+  if (_banner_ostr != 0){
+    (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl;
+    (*_banner_ostr) << "#                    SISCone   version " << setw(28) << left << siscone_version() << "o" << endl;
+    (*_banner_ostr) << "#              http://projects.hepforge.org/siscone                o" << endl;
+    (*_banner_ostr) << "#                                                                  o" << endl;
+    (*_banner_ostr) << "# This is SISCone: the Seedless Infrared Safe Cone Jet Algorithm   o" << endl;
+    (*_banner_ostr) << "# SISCone was written by Gavin Salam and Gregory Soyez             o" << endl;
+    (*_banner_ostr) << "# It is released under the terms of the GNU General Public License o" << endl;
+    (*_banner_ostr) << "#                                                                  o" << endl;
+    (*_banner_ostr) << "# A description of the algorithm is available in the publication   o" << endl;
+    (*_banner_ostr) << "# JHEP 05 (2007) 086 [arXiv:0704.0292 (hep-ph)].                   o" << endl;
+    (*_banner_ostr) << "# Please cite it if you use SISCone.                               o" << endl;
+    (*_banner_ostr) << "#ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" << endl;
+    (*_banner_ostr) << endl;
+
+    _banner_ostr->flush();
+  }
+}
 
 // finally, a bunch of functions to access to 

external/fastjet/plugins/SISCone/siscone.h

@@ -22 +22 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 325                                                          $//
-// $Date:: 2011-11-25 12:41:17 +0100 (Fri, 25 Nov 2011)                     $//
+// $Revision:: 369                                                          $//
+// $Date:: 2014-09-04 16:57:55 +0200 (Thu, 04 Sep 2014)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -79 +79 @@
 
   /**
+   * compute the jets from a given particle set.
+   * We are doing multiple passes such pass n_pass looks for jets among 
+   * all particles not put into jets during previous passes.
+   * By default the number of passes is infinite (0). 
+   * \param _particles   list of particles
+   * \param _radius      cone radius
+   * \param _n_pass_max  maximum number of passes (0=full search)
+   * \param _ptmin       minimum pT of the protojets
+   * \param _ordering_scale    the ordering scale to decide which stable
+   *                           cone is removed
+   *
+   * Note that the Csplit_merge::SM_var2_hardest_cut_off cut is not
+   * used in the progressive removal variant.
+   * 
+   * \return the number of jets found.
+   */
+  int compute_jets_progressive_removal(std::vector<Cmomentum> &_particles, double _radius, 
+                                       int _n_pass_max=0, double _ptmin=0.0,
+                                       Esplit_merge_scale _ordering_scale=SM_pttilde);
+
+  /**
    * recompute the jets with a different overlap parameter.
    * we use the same particles and R as in the preceeding call.
@@ -93 +114 @@
                      Esplit_merge_scale _split_merge_scale=SM_pttilde);
 
-  /// list of protocones found pass-by-pass
+  /// list of protocones found pass-by-pass (not filled by compute_jets_progressive_removal)
   std::vector<std::vector<Cmomentum> > protocones_list;
 
@@ -125 +146 @@
   bool rerun_allowed;         ///< is recompute_jets allowed ?
   static std::ostream * _banner_ostr; ///< stream to use for banners
+
+  /// ensure things are initialised
+  void _initialise_if_needed();
+
 };
 

external/fastjet/plugins/SISCone/split_merge.cc

@@ -21 +21 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 322                                                          $//
-// $Date:: 2011-11-15 10:12:36 +0100 (Tue, 15 Nov 2011)                     $//
+// $Revision:: 370                                                          $//
+// $Date:: 2014-09-04 17:03:15 +0200 (Thu, 04 Sep 2014)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -28 +28 @@
 #include "siscone_error.h"
 #include "momentum.h"
-#include <math.h>
 #include <limits>   // for max
 #include <iostream>
@@ -34 +33 @@
 #include <sstream>
 #include <cassert>
+#include <cmath>
 
 namespace siscone{
@@ -229 +229 @@
 #endif
 #endif
+  _user_scale = NULL;
   indices = NULL;
 
@@ -237 +238 @@
 
   // no hardest cut (col-unsafe)
-  SM_var2_hardest_cut_off = -1.0;
+  SM_var2_hardest_cut_off = -numeric_limits<double>::max();
 
   // no pt cutoff for the particles to put in p_uncol_hard
@@ -555 +556 @@
 }
 
+
+/*
+ * remove the hardest protocone and declare it a a jet 
+ *  - protocones  list of protocones (initial jet candidates)
+ *  - R2          cone radius (squared)
+ *  - ptmin       minimal pT allowed for jets
+ * return 0 on success, 1 on error
+ *
+ * The list of remaining particles (and the uncollinear-hard ones)
+ * is updated.
+ */
+int Csplit_merge::add_hardest_protocone_to_jets(std::vector<Cmomentum> *protocones, double R2, double ptmin){
+
+  int i;
+  Cmomentum *c;
+  Cmomentum *v;
+  double eta, phi;
+  double dx, dy;
+  double R;
+  Cjet jet, jet_candidate;
+  bool found_jet = false;
+
+  if (protocones->size()==0)
+    return 1;
+
+  pt_min2 = ptmin*ptmin;
+  R = sqrt(R2);
+
+  // browse protocones
+  // for each of them, build the list of particles in them
+  for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
+    // initialise variables
+    c = &(*p_it);
+
+    // note: cones have been tested => their (eta,phi) coordinates are computed
+    eta = c->eta;
+    phi = c->phi;
+
+    // NOTE: this is common to this method and add_protocones, so it
+    // could be moved into a 'build_jet_from_protocone' method
+    //
+    // browse particles to create cone contents
+    jet_candidate.v = Cmomentum();
+    jet_candidate.pt_tilde=0;
+    jet_candidate.contents.clear();
+    for (i=0;i<n_left;i++){
+      v = &(p_remain[i]);
+      // for security, avoid including particles with infinite rapidity)
+      // NO NEEDED ANYMORE SINCE REMOVED FROM p_remain FROM THE BEGINNING
+      //if (fabs(v->pz)!=v->E){
+      dx = eta - v->eta;
+      dy = fabs(phi - v->phi);
+      if (dy>M_PI) 
+        dy -= twopi;
+      if (dx*dx+dy*dy<R2){
+        jet_candidate.contents.push_back(v->parent_index);
+        jet_candidate.v+= *v;
+        jet_candidate.pt_tilde+= pt[v->parent_index];
+        v->index=0;
+      }
+    }
+    jet_candidate.n=jet_candidate.contents.size();
+
+    // set the momentum in protocones 
+    // (it was only known through eta and phi up to now)
+    *c = jet_candidate.v;
+    c->eta = eta; // restore exact original coords
+    c->phi = phi; // to avoid rounding error inconsistencies
+
+    // set the jet range
+    jet_candidate.range=Ceta_phi_range(eta,phi,R);
+
+    // check that the protojet has large enough pt
+    if (jet_candidate.v.perp2()<pt_min2)
+      continue;
+
+    // assign the split-merge (or progressive-removal) squared scale variable
+    if (_user_scale) {
+      // sm_var2 is the signed square of the user scale returned
+      // for the jet candidate
+      jet_candidate.sm_var2 = (*_user_scale)(jet_candidate);
+      jet_candidate.sm_var2 *= abs(jet_candidate.sm_var2);
+    } else {
+      jet_candidate.sm_var2 = get_sm_var2(jet_candidate.v, jet_candidate.pt_tilde);
+    }
+
+    // now check if it is possibly the hardest
+    if ((! found_jet) ||
+        (_user_scale ? _user_scale->is_larger(jet_candidate, jet)
+                     : ptcomparison(jet_candidate, jet))){
+      jet = jet_candidate;
+      found_jet = true;
+    }
+  }
+
+  // make sure at least one of the jets has passed the selection
+  if (!found_jet) return 1;  
+  
+  // add the jet to the list of jets
+  jets.push_back(jet);
+  jets[jets.size()-1].v.build_etaphi();
+
+#ifdef DEBUG_SPLIT_MERGE
+  cout << "PR-Jet " << jets.size() << " [size " << next_jet.contents.size() << "]:";
+#endif
+    
+  // update the list of what particles are left
+  int p_remain_index = 0;
+  int contents_index = 0;
+  //sort(next_jet.contents.begin(),next_jet.contents.end());
+  for (int index=0;index<n_left;index++){
+    if ((contents_index<(int) jet.contents.size()) &&
+        (p_remain[index].parent_index == jet.contents[contents_index])){
+      // this particle belongs to the newly found jet
+      // set pass in initial list
+      particles[p_remain[index].parent_index].index = n_pass;
+#ifdef DEBUG_SPLIT_MERGE
+      cout << " " << jet.contents[contents_index];
+#endif
+      contents_index++;
+    } else {
+      // this particle still has to be clustered
+      p_remain[p_remain_index] = p_remain[index];
+      p_remain[p_remain_index].parent_index = p_remain[index].parent_index;
+      p_remain[p_remain_index].index=1;
+      p_remain_index++;
+    }
+  }
+  p_remain.resize(n_left-jet.contents.size());
+  n_left = p_remain.size();
+  jets[jets.size()-1].pass = particles[jet.contents[0]].index;
+
+  // increase the pass index
+  n_pass++;
+
+#ifdef DEBUG_SPLIT_MERGE
+  cout << endl;
+#endif
+
+  // male sure the list of uncol_hard particles (used for the next
+  // stable cone finding) is updated [could probably be more
+  // efficient]
+  merge_collinear_and_remove_soft();
+  
+  return 0;
+}
 
 /*

external/fastjet/plugins/SISCone/split_merge.h

@@ -22 +22 @@
 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
 //                                                                           //
-// $Revision:: 268                                                          $//
-// $Date:: 2009-03-12 21:24:16 +0100 (Thu, 12 Mar 2009)                     $//
+// $Revision:: 367                                                          $//
+// $Date:: 2014-09-04 15:57:37 +0200 (Thu, 04 Sep 2014)                     $//
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -141 +141 @@
   /// the split-merge process i.e. the variable we use for 
   ///  1. ordering jet candidates;
-  ///  2. computing te overlap fraction of two candidates.
+  ///  2. computing the overlap fraction of two candidates.
   /// The default value uses pttile (p-scheme pt). Other alternatives are
   /// pt, mt=sqrt(pt^2+m^2)=sqrt(E^2-pz^2) or Et. 
@@ -151 +151 @@
   ///   the default value i.e.  to use pt only for the purpose of
   ///   investigating the IR issue
-  /// - using Et is safe but do not respect boost invariance
+  /// - using Et is safe but does not respect boost invariance
   /// - using mt solves the IR unsafety issues with the pt variable
   ///   for QCD jets but the IR unsafety remains for nack-to-back 
@@ -234 +234 @@
   int full_clear();
 
+  ///////////////////////////////////////
+  // user-defined stable-cone ordering //
+  ///////////////////////////////////////
+
+  /// \class Cuser_scale_base
+  /// base class for user-defined ordering of stable cones
+  ///
+  /// derived classes have to implement the () operator that returns
+  /// the scale associated with a given jet.
+  class Cuser_scale_base{
+  public:
+    /// empty virtual dtor
+    virtual ~Cuser_scale_base(){}
+
+    /// the scale associated with a given jet
+    ///
+    /// "progressive removal" iteratively removes the stable cone with
+    /// the largest scale
+    virtual double operator()(const Cjet & jet) const = 0;
+
+    /// returns true when the scale associated with jet a is larger than
+    /// the scale associated with jet b
+    ///
+    /// By default this does a simple direct comparison but it can be
+    /// overloaded for higher precision [recommended if possible]
+    ///
+    /// This function assumes that a.sm_var2 and b.sm_var2 have been
+    /// correctly initialised with the signed squared output of
+    /// operator(), as is by default the case when is_larger is called
+    /// from within siscone.
+    virtual bool is_larger(const Cjet & a, const Cjet & b) const{
+      return (a.sm_var2 > b.sm_var2);
+    }
+  };
+
+  /// associate a user-defined scale to order the stable cones
+  ///
+  /// Note that this is only used in "progressive-removal mode",
+  /// e.g. in add_hardest_protocone_to_jets().
+  void set_user_scale(const Cuser_scale_base * user_scale_in){
+    _user_scale = user_scale_in;
+  }
+
+  /// return the user-defined scale (NULL if none)
+  const Cuser_scale_base * user_scale() const { return _user_scale; }
+
 
   /////////////////////////////////
@@ -256 +302 @@
    */
   int add_protocones(std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0);
+
+  /**
+   * remove the hardest protocone and declare it a jet 
+   * \param protocones  list of protocones (initial jet candidates)
+   * \param R2          cone radius (squared)
+   * \param ptmin       minimal pT allowed for jets
+   * \return 0 on success, 1 on error
+   *
+   * The list of remaining particles (and the uncollinear-hard ones)
+   * is updated.
+   */
+  int add_hardest_protocone_to_jets(std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0);
 
   /**
@@ -314 +372 @@
   Csplit_merge_ptcomparison ptcomparison;
 
-  /// stop split--merge when the SM_var of the hardest protojet 
-  /// is below this cut-off. 
+  /// stop split--merge or progressive-removal when the squared SM_var
+  /// of the hardest protojet is below this cut-off. Note that this is
+  /// a signed square (ie SM_var*|SM_var|) to be able to handle
+  /// negative values.
+  ///
+  /// Note that the cut-off is set on the variable squared.
+  double SM_var2_hardest_cut_off;
+
+  /// pt cutoff for the particles to put in p_uncol_hard 
+  /// this is meant to allow removing soft particles in the
+  /// stable-cone search.
+  ///
   /// This is not collinear-safe so you should not use this
   /// variable unless you really know what you are doing
   /// Note that the cut-off is set on the variable squared.
-  double SM_var2_hardest_cut_off;
-
-  /// pt cutoff for the particles to put in p_uncol_hard 
-  /// this is meant to allow removing soft particles in the
-  /// stable-cone search.
   double stable_cone_soft_pt2_cutoff;
 
@@ -390 +453 @@
   bool use_pt_weighted_splitting;
 
+  /// use a user-defined scale to order the stable cones and jet
+  /// candidates
+  const Cuser_scale_base *_user_scale;
+
 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
   /// checkxor for the candidates (to avoid having twice the same contents)