Changeset 273e668 in git for external/fastjet

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

Files:

• AUTHORS (modified) (1 diff)
• ClusterSequence.cc (modified) (4 diffs)
• ClusterSequence.hh (modified) (3 diffs)
• CompositeJetStructure.hh (modified) (4 diffs)
• Error.cc (modified) (5 diffs)
• Error.hh (modified) (4 diffs)
• JetDefinition.cc (modified) (4 diffs)
• JetDefinition.hh (modified) (3 diffs)
• PseudoJet.cc (modified) (4 diffs)
• RectangularGrid.hh (modified) (2 diffs)
• Selector.hh (modified) (3 diffs)
• VERSION (modified) (1 diff)
• config_auto.h (modified) (5 diffs)
• config_raw.h (modified) (5 diffs)
• config_win.h (modified) (1 diff)
• plugins/D0RunICone/fastjet/D0RunIConePlugin.hh (modified) (1 diff)
• plugins/D0RunICone/fastjet/D0RunIpre96ConePlugin.hh (modified) (1 diff)
• plugins/SISCone/SISConePlugin.cc (modified) (6 diffs)
• plugins/SISCone/config.h (modified) (2 diffs)
• plugins/SISCone/fastjet/SISConeBasePlugin.hh (modified) (6 diffs)
• plugins/SISCone/fastjet/SISConePlugin.hh (modified) (6 diffs)
• plugins/SISCone/siscone.cc (modified) (5 diffs)
• plugins/SISCone/siscone.h (modified) (4 diffs)
• plugins/SISCone/split_merge.cc (modified) (6 diffs)
• plugins/SISCone/split_merge.h (modified) (7 diffs)
• tools/Recluster.hh (modified) (2 diffs)
• tools/Subtractor.cc (modified) (6 diffs)
• tools/Subtractor.hh (modified) (5 diffs)

external/fastjet/AUTHORS

@@ -35 +35 @@
 Chris Vermilion
 Markus Wobisch 
+Christopher Young
 
 ----------------------------------------------------------------------

external/fastjet/ClusterSequence.cc

@@ -1 +1 @@
 //FJSTARTHEADER
-// $Id: ClusterSequence.cc 3619 2014-08-13 14:17:19Z salam $
+// $Id: ClusterSequence.cc 3685 2014-09-11 20:15:00Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -37 +37 @@
 #include "fastjet/internal/LazyTiling9.hh"
 #include "fastjet/internal/LazyTiling25.hh"
+#ifndef __FJCORE__
 #include "fastjet/internal/LazyTiling9SeparateGhosts.hh"
+#endif  // __FJCORE__
 #include<iostream>
 #include<sstream>
@@ -370 +372 @@
     _plugin_activated = false;
 
+#ifndef __FJCORE__
   } else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
     // attempt to use an external tiling routine -- it manipulates
@@ -377 +380 @@
     tiling.run();
     _plugin_activated = false;
+#else 
+    throw Error("N2MHTLazy9AntiKtSeparateGhosts strategy not supported with FJCORE");
+#endif  // __FJCORE__
 
   } else if (_strategy == NlnN) {

external/fastjet/ClusterSequence.hh

@@ -3 +3 @@
 
 //FJSTARTHEADER
-// $Id: ClusterSequence.hh 3619 2014-08-13 14:17:19Z salam $
+// $Id: ClusterSequence.hh 3709 2014-09-29 13:19:11Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -641 +641 @@
   Strategy _best_strategy() const;
   
+  /// \if internal_doc
+  /// \class _Parabola
   /// returns c*(a*R**2 + b*R + 1);
   /// Written as a class in case we want to give names to different
   /// parabolas
+  /// \endif
   class _Parabola {
   public:
@@ -652 +655 @@
   };
 
+  /// \if internal_doc
+  /// \class _Line
   /// operator()(R) returns a*R+b;
+  /// \endif
   class _Line {
   public:

external/fastjet/CompositeJetStructure.hh

@@ -1 +1 @@
 //FJSTARTHEADER
-// $Id: CompositeJetStructure.hh 3453 2014-07-25 07:44:32Z soyez $
+// $Id: CompositeJetStructure.hh 3652 2014-09-03 13:31:13Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -244 +244 @@
                                     const JetDefinition::Recombiner & recombiner){
   std::vector<PseudoJet> pieces;
+  pieces.reserve(2);
   pieces.push_back(j1);
   pieces.push_back(j2);
@@ -255 +256 @@
                                     const JetDefinition::Recombiner & recombiner){
   std::vector<PseudoJet> pieces;
+  pieces.reserve(3);
   pieces.push_back(j1);
   pieces.push_back(j2);
@@ -267 +269 @@
                                     const JetDefinition::Recombiner & recombiner){
   std::vector<PseudoJet> pieces;
+  pieces.reserve(4);
   pieces.push_back(j1);
   pieces.push_back(j2);

external/fastjet/Error.cc

@@ -1 +1 @@
 //FJSTARTHEADER
-// $Id: Error.cc 3433 2014-07-23 08:17:03Z salam $
+// $Id: Error.cc 3695 2014-09-18 13:57:56Z cacciari $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -33 +33 @@
 #include <sstream>
 
+#ifndef __FJCORE__
 // printing the stack would need execinfo
 #ifdef FASTJET_HAVE_EXECINFO_H
 #include <execinfo.h>
 #include <cstdlib>
-#endif
+#ifdef FASTJET_HAVE_DEMANGLING_SUPPORT
+#include <cstdio>
+#include <cxxabi.h>
+#endif // FASTJET_HAVE_DEMANGLING_SUPPORT
+#endif // FASTJET_HAVE_EXECINFO_H
+#endif  // __FJCORE__
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
@@ -46 +52 @@
 bool Error::_print_backtrace = false;
 ostream * Error::_default_ostr = & cerr;
+#if (!defined(FASTJET_HAVE_EXECINFO_H)) || defined(__FJCORE__)
+  LimitedWarning Error::_execinfo_undefined;
+#endif
 
+//----------------------------------------------------------------------
+#ifndef __FJCORE__ 
+
+// demangling only is included, i.e. --enable-demangling is specified
+// at configure time, execinfo.h is present and the GNU C++ ABI is
+// supported
+#if defined(FASTJET_HAVE_EXECINFO_H) && defined(FASTJET_HAVE_DEMANGLING_SUPPORT)
+// demangle a given backtrace symbol
+//
+// Notes:
+//  - at the moment, only the symbol is parsed.
+//  - one can get the offset by using 
+//      "%*[^(]%*[^_]%127[^+)]%64[+x0123456789abcdef]", symbol, offset
+//    and checking if sscanf returns 0, 1 or 2
+//    (offset includes the leading +)
+//  - Similarly one could exctract the address and try to convert it
+//    into a filename+line number like addr2line does but this seems
+//    to require exteral dependencies. If we want to go down that
+//    route, one could look into the inplementation o faddr2line(.c)
+//    and/or dladdr.
+string Error::_demangle(const char* symbol) {
+  size_t size;
+  int status;
+  char temp[128];
+  char* demangled;
+  // first, try to demangle a c++ name
+  // decryption:
+  //   %*[^(]  matches any number of characters different from "("
+  //           the * tells not to store in input var 
+  //   %*[^_]  matches any number of characters different from "_"
+  //           the * tells not to store in input var 
+  //   %127[^)+]  matches at most 127 characters different from "+"
+  //              match is stored
+  if (1 == sscanf(symbol, "%*[^(]%*[^_]%127[^)+]", temp)) {
+    //cout << symbol << " -> " << temp << endl;
+    if (NULL != (demangled = abi::__cxa_demangle(temp, NULL, &size, &status))) {
+      string result(demangled);
+      free(demangled);
+      return result;
+    }
+  }
+  //if that didn't work, try to get a regular c symbol
+  if (1 == sscanf(symbol, "%127s", temp)) {
+    return temp;
+  }
+ 
+  //if all else fails, just return the symbol
+  return symbol;
+}
+#endif  // FASTJET_HAVE_DEMANGLING_SUPPORT && FASTJET_HAVE_EXECINFO_H
+#endif  // __FJCORE__
+
+
+//----------------------------------------------------------------------
 Error::Error(const std::string & message_in) {
   _message = message_in; 
+
   if (_print_errors && _default_ostr){
     ostringstream oss;
     oss << "fastjet::Error:  "<< message_in << endl;
 
+#ifndef __FJCORE__
     // only print the stack if execinfo is available and stack enabled
 #ifdef FASTJET_HAVE_EXECINFO_H
@@ -64 +129 @@
       oss << "stack:" << endl;
       for (int i = 1; i < size && messages != NULL; ++i){
+#ifdef FASTJET_HAVE_DEMANGLING_SUPPORT
+        oss << "  #" << i << ": " << _demangle(messages[i])
+            << " [" << messages[i]  << "]" << endl;
+#else
         oss << "  #" << i << ": " << messages[i] << endl;
+#endif
       }
       free(messages);
     }
-#endif
+#endif  // FASTJET_HAVE_EXECINFO_H
+#endif  // __FJCORE__
 
     *_default_ostr << oss.str();
@@ -85 +156 @@
 }
 
+//----------------------------------------------------------------------
+void Error::set_print_backtrace(bool enabled) {
+#if (!defined(FASTJET_HAVE_EXECINFO_H)) || defined(__FJCORE__)
+  if (enabled) {
+    _execinfo_undefined.warn("Error::set_print_backtrace(true) will not work with this build of FastJet");
+  }
+#endif    
+  _print_backtrace = enabled;
+}
+
 FASTJET_END_NAMESPACE
 

external/fastjet/Error.hh

@@ -3 +3 @@
 
 //FJSTARTHEADER
-// $Id: Error.hh 3433 2014-07-23 08:17:03Z salam $
+// $Id: Error.hh 3694 2014-09-18 13:21:54Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -35 +35 @@
 #include<string>
 #include "fastjet/internal/base.hh"
+#include "fastjet/config.h"
+#if (!defined(FASTJET_HAVE_EXECINFO_H)) || defined(__FJCORE__)
+#include "fastjet/LimitedWarning.hh"
+#endif
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
@@ -65 +69 @@
   /// controls whether the backtrace is printed out with the error message or not.
   /// The default is "false".
-  static void set_print_backtrace(bool enabled) {_print_backtrace = enabled;}
+  static void set_print_backtrace(bool enabled);
 
   /// sets the default output stream for all errors; by default
@@ -74 +78 @@
 
 private:
+
+#ifndef __FJCORE__
+#if defined(FASTJET_HAVE_EXECINFO_H) && defined(FASTJET_HAVE_DEMANGLING_SUPPORT)
+  /// demangle a given backtrace symbol
+  std::string _demangle(const char* symbol);
+#endif
+#endif
+
   std::string _message;                ///< error message
   static bool _print_errors;           ///< do we print anything?
   static bool _print_backtrace;        ///< do we print the backtrace?
   static std::ostream * _default_ostr; ///< the output stream (cerr if not set)
+#if (!defined(FASTJET_HAVE_EXECINFO_H)) || defined(__FJCORE__)
+  static LimitedWarning _execinfo_undefined;
+#endif
 };
 

external/fastjet/JetDefinition.cc

@@ -1 +1 @@
 //FJSTARTHEADER
-// $Id: JetDefinition.cc 3492 2014-07-30 16:58:20Z soyez $
+// $Id: JetDefinition.cc 3677 2014-09-09 22:45:25Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -270 +270 @@
   case BIpt2_scheme:
     return "boost-invariant pt2 scheme recombination";
+  case WTA_pt_scheme:
+    return "pt-ordered Winner-Takes-All recombination";
+  // Energy-ordering can lead to dangerous situations with particles at
+  // rest. We instead implement the WTA_modp_scheme
+  //
+  //   case WTA_E_scheme:
+  //     return "energy-ordered Winner-Takes-All recombination";
+  case WTA_modp_scheme:
+    return "|3-momentum|-ordered Winner-Takes-All recombination";
   default:
     ostringstream err;
@@ -309 +318 @@
     weightb = pb.perp2();
     break;
+  case WTA_pt_scheme:{
+    const PseudoJet & phard = (pa.pt2() >= pb.pt2()) ? pa : pb;
+    /// keep y,phi and m from the hardest, sum pt
+    pab.reset_PtYPhiM(pa.pt()+pb.pt(), 
+                      phard.rap(), phard.phi(), phard.m());
+    return;}
+  // Energy-ordering can lead to dangerous situations with particles at
+  // rest. We instead implement the WTA_modp_scheme
+  //
+  //   case WTA_E_scheme:{
+  //     const PseudoJet & phard = (pa.E() >= pb.E()) ? pa : pb;
+  //     /// keep 3-momentum direction and mass from the hardest, sum energies
+  //     ///
+  //     /// If the particle with the largest energy is at rest, the sum
+  //     /// remains at rest, implying that the mass of the sum is larger
+  //     /// than the mass of pa.
+  //     double Eab = pa.E() + pb.E();
+  //     double scale = (phard.modp2()==0.0)
+  //       ? 0.0
+  //       : sqrt((Eab*Eab - phard.m2())/phard.modp2());
+  //     pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale, Eab);
+  //     return;}
+  case WTA_modp_scheme:{
+    // Note: we need to compute both a and b modp. And we need pthard
+    // and its modp. If we want to avoid repeating the test and do
+    // only 2 modp calculations, we'd have to duplicate the code (or
+    // use a pair<const PJ&>). An alternative is to write modp_soft as
+    // modp_ab-modp_hard but this could suffer from larger rounding
+    // errors
+    bool a_hardest = (pa.modp2() >= pb.modp2());
+    const PseudoJet & phard = a_hardest ? pa : pb;
+    const PseudoJet & psoft = a_hardest ? pb : pa;
+    /// keep 3-momentum direction and mass from the hardest, sum modp
+    ///
+    /// If the hardest particle is at rest, the sum remains at rest
+    /// (the energy of the sum is therefore the mass of pa)
+    double modp_hard = phard.modp();
+    double modp_ab = modp_hard + psoft.modp();
+    if (phard.modp2()==0.0){
+      pab.reset(0.0, 0.0, 0.0, phard.m());
+    } else {
+      double scale = modp_ab/modp_hard;
+      pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale,
+                sqrt(modp_ab*modp_ab + phard.m2()));
+    }
+    return;}
   default:
     ostringstream err;
@@ -344 +399 @@
   case BIpt_scheme:
   case BIpt2_scheme:
+  case WTA_pt_scheme:
+  //case WTA_E_scheme:
+  case WTA_modp_scheme:
     break;
   case pt_scheme:

external/fastjet/JetDefinition.hh

@@ -3 +3 @@
 
 //FJSTARTHEADER
-// $Id: JetDefinition.hh 3523 2014-08-02 13:15:21Z salam $
+// $Id: JetDefinition.hh 3677 2014-09-09 22:45:25Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -181 +181 @@
 
 //======================================================================
-/// the various recombination schemes
+/// The various recombination schemes
+///
+/// Note that the schemes that recombine with non-linear weighting of
+/// the directions (e.g. pt2, winner-takes-all) are collinear safe
+/// only for algorithms with a suitable ordering of the
+/// recombinations: orderings in which, for particles of comparable
+/// energies, small-angle clusterings take place before large-angle
+/// clusterings. This property is satisfied by all gen-kt algorithms.
+/// 
 enum RecombinationScheme {
   /// summing the 4-momenta
@@ -203 +211 @@
   /// no preprocessing
   BIpt2_scheme=6,
+  /// pt-based Winner-Takes-All (WTA) recombination: the
+  /// result of the recombination has the rapidity, azimuth and mass
+  /// of the the PseudoJet with the larger pt, and a pt equal to the
+  /// sum of the two pt's
+  WTA_pt_scheme=7,
+  /// mod-p-based Winner-Takes-All (WTA) recombination: the result of
+  /// the recombination gets the 3-vector direction and mass of the
+  /// PseudoJet with the larger |3-momentum| (modp), and a
+  /// |3-momentum| equal to the scalar sum of the two |3-momenta|.
+  WTA_modp_scheme=8,
+  // Energy-ordering can lead to dangerous situations with particles at
+  // rest. We instead implement the WTA_modp_scheme
+  //
+  // // energy-based Winner-Takes-All (WTA) recombination: the result of
+  // // the recombination gets the 3-vector direction and mass of the
+  // // PseudoJet with the larger energy, and an energy equal to the
+  // // to the sum of the two energies
+  // WTA_E_scheme=8,
   /// for the user's external scheme
   external_scheme = 99

external/fastjet/PseudoJet.cc

@@ -1 +1 @@
 //FJSTARTHEADER
-// $Id: PseudoJet.cc 3565 2014-08-11 15:24:39Z salam $
+// $Id: PseudoJet.cc 3652 2014-09-03 13:31:13Z salam $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -861 +861 @@
 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
   vector<PseudoJet> pieces;
+  pieces.reserve(2);
   pieces.push_back(j1);
   pieces.push_back(j2);
@@ -869 +870 @@
 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){
   vector<PseudoJet> pieces;
+  pieces.reserve(3);
   pieces.push_back(j1);
   pieces.push_back(j2);
@@ -878 +880 @@
 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){
   vector<PseudoJet> pieces;
+  pieces.reserve(4);
   pieces.push_back(j1);
   pieces.push_back(j2);

external/fastjet/RectangularGrid.hh

@@ -88 +88 @@
   ///  \param rapmax     the maximal absolute rapidity extent of the grid
   ///  \param cell_size  the grid spacing (equivalently, cell size)
-  RectangularGrid(double rapmax, double cell_size) :
-      _ymax(rapmax), _ymin(-rapmax), 
+  RectangularGrid(double rapmax_in, double cell_size) :
+      _ymax(rapmax_in), _ymin(-rapmax_in), 
       _requested_drap(cell_size), _requested_dphi(cell_size) {
     _setup_grid();
@@ -103 +103 @@
   ///                        a massless 4-vector at the center of the tile passes
   ///                        the selection
-  RectangularGrid(double rapmin, double rapmax, double drap, double dphi,
+  RectangularGrid(double rapmin_in, double rapmax_in, double drap_in, double dphi_in,
                   Selector tile_selector = Selector()) 
-    : _ymax(rapmax), _ymin(rapmin), 
-      _requested_drap(drap), _requested_dphi(dphi),
+    : _ymax(rapmax_in), _ymin(rapmin_in), 
+      _requested_drap(drap_in), _requested_dphi(dphi_in),
       _tile_selector(tile_selector)
   {

external/fastjet/Selector.hh

@@ -3 +3 @@
 
 //FJSTARTHEADER
-// $Id: Selector.hh 3504 2014-08-01 06:07:54Z soyez $
+// $Id: Selector.hh 3711 2014-09-29 13:54:51Z salam $
 //
 // Copyright (c) 2009-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -305 +305 @@
   }
 
-  /// class that gets throw when a Selector is applied despite it not
+  /// class that gets thrown when a Selector is applied despite it not
   /// having a valid underlying worker.
   class InvalidWorker : public Error {
@@ -312 +312 @@
   };
 
-  /// class that gets throw when a Selector is applied despite it not
-  /// having a valid underlying worker.
+  /// class that gets thrown when the area is requested from a Selector for which
+  /// the area is not meaningful
   class InvalidArea : public Error {
   public:

external/fastjet/VERSION

@@ -1 @@
-3.1.0-beta.1
+3.1.0

external/fastjet/config_auto.h

@@ -64 +64 @@
 #endif
 
+/* defined if demangling is enabled at configure time and is supported through
+   the GNU C++ ABI */
+/* #undef HAVE_DEMANGLING_SUPPORT */
+
 /* Define to 1 if you have the <dlfcn.h> header file. */
 #ifndef FASTJET_HAVE_DLFCN_H 
@@ -147 +151 @@
 /* Define to the full name and version of this package. */
 #ifndef FASTJET_PACKAGE_STRING 
-#define FASTJET_PACKAGE_STRING  "FastJet 3.1.0-beta.1" 
+#define FASTJET_PACKAGE_STRING  "FastJet 3.1.0" 
 #endif
 
@@ -157 +161 @@
 /* Define to the version of this package. */
 #ifndef FASTJET_PACKAGE_VERSION 
-#define FASTJET_PACKAGE_VERSION  "3.1.0-beta.1" 
+#define FASTJET_PACKAGE_VERSION  "3.1.0" 
 #endif
 
@@ -167 +171 @@
 /* Version number of package */
 #ifndef FASTJET_VERSION 
-#define FASTJET_VERSION  "3.1.0-beta.1" 
+#define FASTJET_VERSION  "3.1.0" 
 #endif
 
@@ -182 +186 @@
 /* Version of the package under the form XYYZZ (instead of X.Y.Z) */
 #ifndef FASTJET_VERSION_NUMBER 
-#define FASTJET_VERSION_NUMBER  30103 
+#define FASTJET_VERSION_NUMBER  30100 
 #endif
 
 /* Patch version of this package */
 #ifndef FASTJET_VERSION_PATCHLEVEL 
-#define FASTJET_VERSION_PATCHLEVEL  3 
+#define FASTJET_VERSION_PATCHLEVEL  0 
 #endif
 
 /* Pre-release version of this package */
-#ifndef FASTJET_VERSION_PRERELEASE 
-#define FASTJET_VERSION_PRERELEASE  ".1.0-beta.1" 
-#endif
+/* #undef VERSION_PRERELEASE */
  
 /* once: _INCLUDE_FASTJET_CONFIG_AUTO_H */

external/fastjet/config_raw.h

@@ -37 +37 @@
 /* The TrackJet plugin is enabled */
 #define ENABLE_PLUGIN_TRACKJET /**/
+
+/* defined if demangling is enabled at configure time and is supported through
+   the GNU C++ ABI */
+/* #undef HAVE_DEMANGLING_SUPPORT */
 
 /* Define to 1 if you have the <dlfcn.h> header file. */
@@ -88 +92 @@
 
 /* Define to the full name and version of this package. */
-#define PACKAGE_STRING "FastJet 3.1.0-beta.1"
+#define PACKAGE_STRING "FastJet 3.1.0"
 
 /* Define to the one symbol short name of this package. */
@@ -94 +98 @@
 
 /* Define to the version of this package. */
-#define PACKAGE_VERSION "3.1.0-beta.1"
+#define PACKAGE_VERSION "3.1.0"
 
 /* Define to 1 if you have the ANSI C header files. */
@@ -100 +104 @@
 
 /* Version number of package */
-#define VERSION "3.1.0-beta.1"
+#define VERSION "3.1.0"
 
 /* Major version of this package */
@@ -109 +113 @@
 
 /* Version of the package under the form XYYZZ (instead of X.Y.Z) */
-#define VERSION_NUMBER 30103
+#define VERSION_NUMBER 30100
 
 /* Patch version of this package */
-#define VERSION_PATCHLEVEL 3
+#define VERSION_PATCHLEVEL 0
 
 /* Pre-release version of this package */
-#define VERSION_PRERELEASE ".1.0-beta.1"
+/* #undef VERSION_PRERELEASE */

external/fastjet/config_win.h

@@ -1 @@
-#define FASTJET_PACKAGE_STRING  "FastJet 3.1.0-beta.1"
-#define FASTJET_PACKAGE_VERSION  "3.1.0-beta.1"
-#define FASTJET_VERSION  "3.1.0-beta.1"
+#define FASTJET_PACKAGE_STRING  "FastJet 3.1.0"
+#define FASTJET_PACKAGE_VERSION  "3.1.0"
+#define FASTJET_VERSION  "3.1.0"
 #define FASTJET_VERSION_MAJOR       3
 #define FASTJET_VERSION_MINOR       1
-#define FASTJET_VERSION_PATCHLEVEL  3
-#define FASTJET_VERSION_PRERELEASE  ".1.0-beta.1"
-#define FASTJET_VERSION_NUMBER      30103
+#define FASTJET_VERSION_PATCHLEVEL  0
+#define FASTJET_VERSION_NUMBER      30100
 
 /* The ATLASCone plugin is disabled by default*/

external/fastjet/plugins/D0RunICone/fastjet/D0RunIConePlugin.hh

@@ -32 +32 @@
 //FJENDHEADER
 
+#include "fastjet/internal/base.hh" // namespace macros (include explicitly to help Doxygen)
 #include "D0RunIBaseConePlugin.hh"
 

external/fastjet/plugins/D0RunICone/fastjet/D0RunIpre96ConePlugin.hh

@@ -32 +32 @@
 //FJENDHEADER
 
+#include "fastjet/internal/base.hh" // namespace macros (include explicitly to help Doxygen)
 #include "D0RunIBaseConePlugin.hh"
 

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)

external/fastjet/tools/Recluster.hh

@@ -2 +2 @@
 #define __FASTJET_TOOLS_RECLUSTER_HH__
 
-// $Id: Recluster.hh 3632 2014-08-15 10:42:53Z salam $
+// $Id: Recluster.hh 3714 2014-09-30 09:47:31Z soyez $
 //
 // Copyright (c) 2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -38 +38 @@
 
 //----------------------------------------------------------------------
+/// @ingroup tools_generic
 /// \class Recluster
 /// Recluster a jet's constituents with a new jet definition.

external/fastjet/tools/Subtractor.cc

@@ -1 +1 @@
 //FJSTARTHEADER
-// $Id: Subtractor.cc 3594 2014-08-12 15:25:11Z soyez $
+// $Id: Subtractor.cc 3670 2014-09-08 14:17:59Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -43 +43 @@
 // ctor
 Subtractor::Subtractor(double rho) : _bge(0), _rho(rho) {
-  assert(_rho>0.0);
+  if (_rho<0.0) throw Error("Subtractor(rho) was passed a negative rho value; rho should be >= 0");
   set_defaults();
 }
 
 //----------------------------------------------------------------------
+// ctor
+Subtractor::Subtractor(double rho, double rho_m) : _bge(0), _rho(rho) {
+  if (_rho<0.0) throw Error("Subtractor(rho, rho_m) was passed a negative rho value; rho should be >= 0");
+  if (rho_m<0.0) throw Error("Subtractor(rho, rho_m) was passed a negative rho_m value; rho_m should be >= 0");
+  set_defaults();
+  _rho_m = rho_m;
+  set_use_rho_m(true);
+}
+
+//----------------------------------------------------------------------
 void Subtractor::set_defaults(){
+  _rho_m = _invalid_rho;
   _use_rho_m = false; // likely to change in future releases!!
   _safe_mass = false; // likely to change in future releases!!
@@ -60 +71 @@
 PseudoJet Subtractor::result(const PseudoJet & jet) const {
   if (!jet.has_area()){
-    throw Error("Trying to subtract a jet without area support");
+    throw Error("Subtractor::result(...): Trying to subtract a jet without area support");
   }
 
@@ -152 +163 @@
     ostringstream ostr;
     ostr << "Subtractor that uses a fixed value of rho = " << _rho;
+    if (use_rho_m()) ostr << " and rho_m = " << _rho_m;
     return ostr.str();
   } else {
@@ -169 +181 @@
     rho = _rho;
   } else {
-    throw Error("default Subtractor does not have any information about the background, which is needed to perform the subtraction");
+    throw Error("Subtractor::_amount_to_subtract(...): default Subtractor does not have any information about the background, needed to perform the subtraction");
   }
 
@@ -178 +190 @@
 
   // add an optional contribution from the unknown particles masses
-  if (_use_rho_m){
-    assert(_bge != 0); // test done in "set_use_rho_m()"
-    to_subtract += _bge->rho_m(jet) * PseudoJet(0.0, 0.0, area.pz(), area.E());
+  if (_use_rho_m) {
+    double rho_m;
+    
+    if (_bge != 0) {
+      if (!_bge->has_rho_m()) throw Error("Subtractor::_amount_to_subtract(...): requested subtraction with rho_m from a background estimator, but the estimator does not have rho_m support");
+      rho_m = _bge->rho_m(jet);
+    } else if (_rho_m != _invalid_rho) {
+      rho_m = _rho_m;
+    } else {
+      throw Error("Subtractor::_amount_to_subtract(...): default Subtractor does not have any information about the background rho_m, needed to perform the rho_m subtraction");
+    }
+    to_subtract += rho_m * PseudoJet(0.0, 0.0, area.pz(), area.E());
   } else if (_bge && 
              _bge->has_rho_m() && 
              _bge->rho_m(jet) > rho_m_warning_threshold * rho) {
-    _unused_rho_m_warning.warn("Background estimator indicates non-zero rho_m, but use_rho_m()==false in subtractor; consider calling set_use_rho_m(true) to include the rho_m information");
+    _unused_rho_m_warning.warn("Subtractor::_amount_to_subtract(...): Background estimator indicates non-zero rho_m, but use_rho_m()==false in subtractor; consider calling set_use_rho_m(true) to include the rho_m information");
   }
 

external/fastjet/tools/Subtractor.hh

@@ -1 +1 @@
 //FJSTARTHEADER
-// $Id: Subtractor.hh 3584 2014-08-12 12:40:10Z soyez $
+// $Id: Subtractor.hh 3670 2014-09-08 14:17:59Z soyez $
 //
 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
@@ -32 +32 @@
 #define __FASTJET_TOOLS_SUBTRACTOR_HH__
 
+#include "fastjet/internal/base.hh"     // namespace macros (include explicitly to help Doxygen)
 #include "fastjet/tools/Transformer.hh" // to derive Subtractor from Transformer
 #include "fastjet/tools/BackgroundEstimatorBase.hh" // used as a ctor argument
@@ -69 +70 @@
   Subtractor(double rho);
 
+  /// define a subtractor that uses a fixed value of rho and rho_m;
+  /// both must be >= 0;
+  Subtractor(double rho, double rho_m);
+
   /// default constructor
   Subtractor() : _bge(0), _rho(_invalid_rho) { set_defaults(); }
@@ -96 +101 @@
   /// in a future release of FastJet
   void set_use_rho_m(bool use_rho_m_in = true){
-    if (_bge == 0) {
-      throw Error("Subtractor: rho_m support works only for Subtractors constructed with background estimator");
+    if (_bge == 0  && _rho_m < 0) {
+      throw Error("Subtractor: rho_m support works only for Subtractors constructed with a background estimator or an explicit rho_m value");
     }
     _use_rho_m=use_rho_m_in;
@@ -176 +181 @@
   /// if has to be mutable in case its underlying selector takes a reference jet
   mutable BackgroundEstimatorBase * _bge;
-  /// the fixed value of rho to use if the user has selected that option
-  double _rho;
+  /// the fixed value of rho and/or rho_m to use if the user has selected that option
+  double _rho, _rho_m;
 
   // configuration parameters/flags